diff --git a/docs/source/io_formats/depletion_results.rst b/docs/source/io_formats/depletion_results.rst
index 7035fc9c9ee..b2a726dad04 100644
--- a/docs/source/io_formats/depletion_results.rst
+++ b/docs/source/io_formats/depletion_results.rst
@@ -4,7 +4,7 @@
 Depletion Results File Format
 =============================
 
-The current version of the depletion results file format is 1.1.
+The current version of the depletion results file format is 1.2.
 
 **/**
 
@@ -12,22 +12,20 @@ The current version of the depletion results file format is 1.1.
              - **version** (*int[2]*) -- Major and minor version of the
                statepoint file format.
 
-:Datasets: - **eigenvalues** (*double[][][2]*) -- k-eigenvalues at each
-             time/stage. This array has shape (number of timesteps, number of
-             stages, value). The last axis contains the eigenvalue and the
-             associated uncertainty
-           - **number** (*double[][][][]*) -- Total number of atoms. This array
-             has shape (number of timesteps, number of stages, number of
+:Datasets: - **eigenvalues** (*double[][2]*) -- k-eigenvalues at each timestep.
+             This array has shape (number of timesteps, 2). The second axis
+             contains the eigenvalue and its associated uncertainty.
+           - **number** (*double[][][]*) -- Total number of atoms at each
+             timestep. This array has shape (number of timesteps, number of
              materials, number of nuclides).
-           - **reaction rates** (*double[][][][][]*) -- Reaction rates used to
-             build depletion matrices. This array has shape (number of
-             timesteps, number of stages, number of materials, number of
-             nuclides, number of reactions).
+           - **reaction rates** (*double[][][][]*) -- Reaction rates at each
+             timestep. This array has shape (number of timesteps, number of
+             materials, number of nuclides, number of reactions). Only stored if
+             write_rates=True.
            - **time** (*double[][2]*) -- Time in [s] at beginning/end of each
              step.
-           - **source_rate** (*double[][]*) -- Power in [W] or source rate in
-             [neutron/sec]. This array has shape (number of timesteps, number
-             of stages).
+           - **source_rate** (*double[]*) -- Power in [W] or source rate in
+             [neutron/sec] for each timestep.
            - **depletion time** (*double[]*) -- Average process time in [s]
              spent depleting a material across all burnable materials and,
              if applicable, MPI processes.
diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index c5a6219c7b8..5240daaa93a 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -631,17 +631,7 @@ def __init__(
             solver: str = "cram48",
             continue_timesteps: bool = False,
         ):
-        # Check number of stages previously used
-        if operator.prev_res is not None:
-            res = operator.prev_res[-1]
-            if res.data.shape[0] != self._num_stages:
-                raise ValueError(
-                    "{} incompatible with previous restart calculation. "
-                    "Previous scheme used {} intermediate solutions, while "
-                    "this uses {}".format(
-                        self.__class__.__name__, res.data.shape[0],
-                        self._num_stages))
-        elif continue_timesteps:
+        if continue_timesteps and operator.prev_res is None:
             raise ValueError("Continuation run requires passing prev_results.")
         self.operator = operator
         self.chain = operator.chain
@@ -775,12 +765,8 @@ def __call__(
         -------
         proc_time : float
             Time spent in CRAM routines for all materials in [s]
-        n_list : list of list of numpy.ndarray
-            Concentrations at each of the intermediate points with
-            the final concentration as the last element
-        op_results : list of openmc.deplete.OperatorResult
-            Eigenvalue and reaction rates from intermediate transport
-            simulations
+        n_end : list of numpy.ndarray
+            Concentrations at end of timestep
         """
 
     @property
@@ -811,9 +797,9 @@ def _get_bos_data_from_restart(self, source_rate, bos_conc):
         """Get beginning of step concentrations, reaction rates from restart"""
         res = self.operator.prev_res[-1]
         # Depletion methods expect list of arrays
-        bos_conc = list(res.data[0])
-        rates = res.rates[0]
-        k = ufloat(res.k[0, 0], res.k[0, 1])
+        bos_conc = list(res.data)
+        rates = res.rates
+        k = ufloat(res.k[0], res.k[1])
 
         if res.source_rate != 0.0:
             # Scale reaction rates by ratio of source rates
@@ -855,7 +841,8 @@ def integrate(
             self,
             final_step: bool = True,
             output: bool = True,
-            path: PathLike = 'depletion_results.h5'
+            path: PathLike = 'depletion_results.h5',
+            write_rates: bool = False
         ):
         """Perform the entire depletion process across all steps
 
@@ -874,6 +861,11 @@ def integrate(
             Path to file to write. Defaults to 'depletion_results.h5'.
 
             .. versionadded:: 0.15.0
+        write_rates : bool, optional
+            Whether reaction rates should be written to the results file for
+            each step. Defaults to ``False`` to reduce file size.
+
+            .. versionadded:: 0.15.3
         """
         with change_directory(self.operator.output_dir):
             n = self.operator.initial_condition()
@@ -890,18 +882,22 @@ def integrate(
                     n, res = self._get_bos_data_from_restart(source_rate, n)
 
                 # Solve Bateman equations over time interval
-                proc_time, n_list, res_list = self(n, res.rates, dt, source_rate, i)
-
-                # Insert BOS concentration, transport results
-                n_list.insert(0, n)
-                res_list.insert(0, res)
-
-                # Remove actual EOS concentration for next step
-                n = n_list.pop()
-
-                StepResult.save(self.operator, n_list, res_list, [t, t + dt],
-                                source_rate, self._i_res + i, proc_time, path)
+                proc_time, n_end = self(n, res.rates, dt, source_rate, i)
+
+                StepResult.save(
+                    self.operator,
+                    n,
+                    res,
+                    [t, t + dt],
+                    source_rate,
+                    self._i_res + i,
+                    proc_time,
+                    write_rates=write_rates,
+                    path=path
+                )
 
+                # Update for next step
+                n = n_end
                 t += dt
 
             # Final simulation -- in the case that final_step is False, a zero
@@ -910,9 +906,18 @@ def integrate(
             # solve)
             if output and final_step and comm.rank == 0:
                 print(f"[openmc.deplete] t={t} (final operator evaluation)")
-            res_list = [self.operator(n, source_rate if final_step else 0.0)]
-            StepResult.save(self.operator, [n], res_list, [t, t],
-                         source_rate, self._i_res + len(self), proc_time, path)
+            res_final = self.operator(n, source_rate if final_step else 0.0)
+            StepResult.save(
+                self.operator,
+                n,
+                res_final,
+                [t, t],
+                source_rate,
+                self._i_res + len(self),
+                proc_time,
+                write_rates=write_rates,
+                path=path
+            )
             self.operator.write_bos_data(len(self) + self._i_res)
 
         self.operator.finalize()
@@ -1141,10 +1146,40 @@ def _get_bos_data_from_operator(self, step_index, step_power, n_bos):
             self.operator.settings.particles //= self.n_steps
         return inherited
 
+    @abstractmethod
+    def __call__(self, n, rates, dt, source_rate, i):
+        """Perform the integration across one time step
+
+        Parameters
+        ----------
+        n : list of numpy.ndarray
+            List of atom number arrays for each material. Each array has
+            shape ``(n_nucs,)`` where ``n_nucs`` is the number of nuclides
+        rates : openmc.deplete.ReactionRates
+            Reaction rates (from transport operator)
+        dt : float
+            Time step in [s]
+        source_rate : float
+            Power in [W] or source rate in [neutron/sec]
+        i : int
+            Current time step index
+
+        Returns
+        -------
+        proc_time : float
+            Time spent in transport simulation
+        n_end : list of numpy.ndarray
+            Updated atom number densities for each material
+        op_result : OperatorResult
+            Eigenvalue and reaction rates resulting from transport simulation
+
+        """
+
     def integrate(
         self,
         output: bool = True,
-        path: PathLike = "depletion_results.h5"
+        path: PathLike = "depletion_results.h5",
+        write_rates: bool = False
     ):
         """Perform the entire depletion process across all steps
 
@@ -1156,11 +1191,17 @@ def integrate(
             Path to file to write. Defaults to 'depletion_results.h5'.
 
             .. versionadded:: 0.15.0
+        write_rates : bool, optional
+            Whether reaction rates should be written to the results file for
+            each step. Defaults to ``False`` to reduce file size.
+
+            .. versionadded:: 0.15.3
         """
         with change_directory(self.operator.output_dir):
             n = self.operator.initial_condition()
             t, self._i_res = self._get_start_data()
 
+            res_end = None  # Will be set in first iteration
             for i, (dt, p) in enumerate(self):
                 if output:
                     print(f"[openmc.deplete] t={t} s, dt={dt} s, source={p}")
@@ -1170,28 +1211,38 @@ def integrate(
                         n, res = self._get_bos_data_from_operator(i, p, n)
                     else:
                         n, res = self._get_bos_data_from_restart(p, n)
-                else:
-                    # Pull rates, k from previous iteration w/o
-                    # re-running transport
-                    res = res_list[-1]  # defined in previous i iteration
-
-                proc_time, n_list, res_list = self(n, res.rates, dt, p, i)
 
-                # Insert BOS concentration, transport results
-                n_list.insert(0, n)
-                res_list.insert(0, res)
-
-                # Remove actual EOS concentration for next step
-                n = n_list.pop()
-
-                StepResult.save(self.operator, n_list, res_list, [t, t + dt],
-                             p, self._i_res + i, proc_time, path)
+                proc_time, n_end, res_end = self(n, res.rates, dt, p, i)
+
+                StepResult.save(
+                    self.operator,
+                    n,
+                    res,
+                    [t, t + dt],
+                    p,
+                    self._i_res + i,
+                    proc_time,
+                    write_rates=write_rates,
+                    path=path
+                )
 
+                # Update for next step
+                n = n_end
+                res = res_end
                 t += dt
 
             # No final simulation for SIE, use last iteration results
-            StepResult.save(self.operator, [n], [res_list[-1]], [t, t],
-                         p, self._i_res + len(self), proc_time, path)
+            StepResult.save(
+                self.operator,
+                n,
+                res_end,
+                [t, t],
+                p,
+                self._i_res + len(self),
+                proc_time,
+                write_rates=write_rates,
+                path=path
+            )
             self.operator.write_bos_data(self._i_res + len(self))
 
         self.operator.finalize()
diff --git a/openmc/deplete/integrators.py b/openmc/deplete/integrators.py
index 7c543a6cb86..25e64cb2ef8 100644
--- a/openmc/deplete/integrators.py
+++ b/openmc/deplete/integrators.py
@@ -46,15 +46,12 @@ def __call__(self, n, rates, dt, source_rate, _i=None):
         -------
         proc_time : float
             Time spent in CRAM routines for all materials in [s]
-        n_list : list of list of numpy.ndarray
+        n_end : list of numpy.ndarray
             Concentrations at end of interval
-        op_results : empty list
-            Kept for consistency with API. No intermediate calls to operator
-            with predictor
 
         """
         proc_time, n_end = self._timed_deplete(n, rates, dt, _i)
-        return proc_time, [n_end], []
+        return proc_time, n_end
 
 
 @add_params
@@ -98,11 +95,8 @@ def __call__(self, n, rates, dt, source_rate, _i=None):
         -------
         proc_time : float
             Time spent in CRAM routines for all materials in [s]
-        n_list : list of list of numpy.ndarray
-            Concentrations at each of the intermediate points with
-            the final concentration as the last element
-        op_results : list of openmc.deplete.OperatorResult
-            Eigenvalue and reaction rates from transport simulations
+        n_end : list of numpy.ndarray
+            Concentrations at end of interval
         """
         # deplete across first half of interval
         time0, n_middle = self._timed_deplete(n, rates, dt / 2, _i)
@@ -112,7 +106,7 @@ def __call__(self, n, rates, dt, source_rate, _i=None):
         # MOS reaction rates
         time1, n_end = self._timed_deplete(n, res_middle.rates, dt, _i)
 
-        return time0 + time1, [n_middle, n_end], [res_middle]
+        return time0 + time1, n_end
 
 
 @add_params
@@ -162,12 +156,8 @@ def __call__(self, n_bos, bos_rates, dt, source_rate, _i=None):
         -------
         proc_time : float
             Time spent in CRAM routines for all materials in [s]
-        n_list : list of numpy.ndarray
-            Concentrations at each of the intermediate points with
-            the final concentration as the last element
-        op_results : list of openmc.deplete.OperatorResult
-            Eigenvalue and reaction rates from intermediate transport
-            simulations
+        n_end : list of numpy.ndarray
+            Concentrations at end of interval
         """
         # Step 1: deplete with matrix 1/2*A(y0)
         time1, n_eos1 = self._timed_deplete(
@@ -192,9 +182,7 @@ def __call__(self, n_bos, bos_rates, dt, source_rate, _i=None):
         time5, n_eos5 = self._timed_deplete(
             n_inter, list_rates, dt, _i, matrix_func=cf4_f4)
 
-        return (time1 + time2 + time3 + time4 + time5,
-                [n_eos1, n_eos2, n_eos3, n_eos5],
-                [res1, res2, res3])
+        return time1 + time2 + time3 + time4 + time5, n_eos5
 
 
 @add_params
@@ -240,12 +228,8 @@ def __call__(self, n_bos, rates, dt, source_rate, _i=None):
         -------
         proc_time : float
             Time spent in CRAM routines for all materials in [s]
-        n_list : list of list of numpy.ndarray
-            Concentrations at each of the intermediate points with
-            the final concentration as the last element
-        op_results : list of openmc.deplete.OperatorResult
-            Eigenvalue and reaction rates from intermediate transport
-            simulation
+        n_end : list of numpy.ndarray
+            Concentrations at end of interval
         """
         # deplete to end using BOS rates
         proc_time, n_ce = self._timed_deplete(n_bos, rates, dt, _i)
@@ -260,7 +244,7 @@ def __call__(self, n_bos, rates, dt, source_rate, _i=None):
         time_le2, n_end = self._timed_deplete(
             n_inter, list_rates, dt, _i, matrix_func=celi_f2)
 
-        return proc_time + time_le1 + time_le1, [n_ce, n_end], [res_ce]
+        return proc_time + time_le1 + time_le2, n_end
 
 
 @add_params
@@ -306,12 +290,8 @@ def __call__(self, n, rates, dt, source_rate, _i=None):
         -------
         proc_time : float
             Time spent in CRAM routines for all materials in [s]
-        n_list : list of list of numpy.ndarray
-            Concentrations at each of the intermediate points with
-            the final concentration as the last element
-        op_results : list of openmc.deplete.OperatorResult
-            Eigenvalue and reaction rates from intermediate transport
-            simulations
+        n_end : list of numpy.ndarray
+            Concentrations at end of interval
         """
 
         # Step 1: deplete with matrix A(y0) / 2
@@ -330,7 +310,7 @@ def __call__(self, n, rates, dt, source_rate, _i=None):
         list_rates = list(zip(rates, res1.rates, res2.rates, res3.rates))
         time4, n4 = self._timed_deplete(n, list_rates, dt, _i, matrix_func=rk4_f4)
 
-        return (time1 + time2 + time3 + time4, [n1, n2, n3, n4], [res1, res2, res3])
+        return time1 + time2 + time3 + time4, n4
 
 
 @add_params
@@ -388,12 +368,8 @@ def __call__(self, n_bos, bos_rates, dt, source_rate, i):
         -------
         proc_time : float
             Time spent in CRAM routines for all materials in [s]
-        n_list : list of list of numpy.ndarray
-            Concentrations at each of the intermediate points with
-            the final concentration as the last element
-        op_results : list of openmc.deplete.OperatorResult
-            Eigenvalue and reaction rates from intermediate transport
-            simulation
+        n_end : list of numpy.ndarray
+            Concentrations at end of interval
         """
         if i == 0:
             if self._i_res < 1:  # need at least previous transport solution
@@ -402,7 +378,7 @@ def __call__(self, n_bos, bos_rates, dt, source_rate, i):
                     self, n_bos, bos_rates, dt, source_rate, i)
             prev_res = self.operator.prev_res[-2]
             prev_dt = self.timesteps[i] - prev_res.time[0]
-            self._prev_rates = prev_res.rates[0]
+            self._prev_rates = prev_res.rates
         else:
             prev_dt = self.timesteps[i - 1]
 
@@ -431,9 +407,7 @@ def __call__(self, n_bos, bos_rates, dt, source_rate, i):
         # store updated rates
         self._prev_rates = copy.deepcopy(bos_res.rates)
 
-        return (
-            time1 + time2 + time3 + time4, [n_eos0, n_eos1],
-            [bos_res, res_inter])
+        return time1 + time2 + time3 + time4, n_eos1
 
 
 @add_params
@@ -470,10 +444,9 @@ def __call__(self, n_bos, bos_rates, dt, source_rate, _i=None):
         -------
         proc_time : float
             Time spent in CRAM routines for all materials in [s]
-        n_bos_list : list of list of numpy.ndarray
-            Concentrations at each of the intermediate points with
-            the final concentration as the last element
-        op_results : list of openmc.deplete.OperatorResult
+        n_end : list of numpy.ndarray
+            Concentrations at end of interval
+        op_result : openmc.deplete.OperatorResult
             Eigenvalue and reaction rates from intermediate transport
             simulations
         """
@@ -499,7 +472,7 @@ def __call__(self, n_bos, bos_rates, dt, source_rate, _i=None):
             proc_time += time1 + time2
 
         # end iteration
-        return proc_time, [n_eos, n_inter], [res_bar]
+        return proc_time, n_inter, res_bar
 
 
 @add_params
@@ -536,10 +509,9 @@ def __call__(self, n_bos, bos_rates, dt, source_rate, i):
         -------
         proc_time : float
             Time spent in CRAM routines for all materials in [s]
-        n_list : list of list of numpy.ndarray
-            Concentrations at each of the intermediate points with
-            the final concentration as the last element
-        op_results : list of openmc.deplete.OperatorResult
+        n_end : list of numpy.ndarray
+            Concentrations at end of interval
+        op_result : openmc.deplete.OperatorResult
             Eigenvalue and reaction rates from intermediate transport
             simulation
         """
@@ -551,7 +523,7 @@ def __call__(self, n_bos, bos_rates, dt, source_rate, i):
                     self, n_bos, bos_rates, dt, source_rate, i)
             prev_res = self.operator.prev_res[-2]
             prev_dt = self.timesteps[i] - prev_res.time[0]
-            self._prev_rates = prev_res.rates[0]
+            self._prev_rates = prev_res.rates
         else:
             prev_dt = self.timesteps[i - 1]
 
@@ -584,7 +556,10 @@ def __call__(self, n_bos, bos_rates, dt, source_rate, i):
                 n_inter, inputs, dt, i, matrix_func=leqi_f4)
             proc_time += time1 + time2
 
-        return proc_time, [n_eos, n_inter], [res_bar]
+        # Store updated rates for next step
+        self._prev_rates = copy.deepcopy(bos_rates)
+
+        return proc_time, n_inter, res_bar
 
 
 integrator_by_name = {
diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py
index 7427abd7353..e1fcb26b6d6 100644
--- a/openmc/deplete/results.py
+++ b/openmc/deplete/results.py
@@ -203,7 +203,7 @@ def get_atoms(
         # Evaluate value in each region
         for i, result in enumerate(self):
             times[i] = result.time[0]
-            concentrations[i] = result[0, mat_id, nuc]
+            concentrations[i] = result[mat_id, nuc]
 
         # Unit conversions
         times = _get_time_as(times, time_units)
@@ -363,7 +363,7 @@ def get_reaction_rate(
         # Evaluate value in each region
         for i, result in enumerate(self):
             times[i] = result.time[0]
-            rates[i] = result.rates[0].get(mat_id, nuc, rx) * result[0, mat, nuc]
+            rates[i] = result.rates.get(mat_id, nuc, rx) * result[mat, nuc]
 
         return times, rates
 
@@ -397,7 +397,7 @@ def get_keff(self, time_units: str = 's') -> tuple[np.ndarray, np.ndarray]:
         # Get time/eigenvalue at each point
         for i, result in enumerate(self):
             times[i] = result.time[0]
-            eigenvalues[i] = result.k[0]
+            eigenvalues[i] = result.k
 
         # Convert time units if necessary
         times = _get_time_as(times, time_units)
@@ -630,7 +630,7 @@ def export_to_materials(
                 for nuc in result.index_nuc:
                     if nuc not in available_cross_sections:
                         continue
-                    atoms = result[0, mat_id, nuc]
+                    atoms = result[mat_id, nuc]
                     if atoms > 0.0:
                         atoms_per_barn_cm = 1e-24 * atoms / mat.volume
                         mat.remove_nuclide(nuc) # Replace if it's there
diff --git a/openmc/deplete/stepresult.py b/openmc/deplete/stepresult.py
index 1a26cbe3466..ff39e9acb6d 100644
--- a/openmc/deplete/stepresult.py
+++ b/openmc/deplete/stepresult.py
@@ -16,7 +16,7 @@
 from openmc.checkvalue import PathLike
 from .reaction_rates import ReactionRates
 
-VERSION_RESULTS = (1, 1)
+VERSION_RESULTS = (1, 2)
 
 
 __all__ = ["StepResult"]
@@ -30,8 +30,8 @@ class StepResult:
 
     Attributes
     ----------
-    k : list of (float, float)
-        Eigenvalue and uncertainty for each substep.
+    k : tuple of (float, float)
+        Eigenvalue and uncertainty at end of step.
     time : list of float
         Time at beginning, end of step, in seconds.
     source_rate : float
@@ -40,8 +40,8 @@ class StepResult:
         Number of mats.
     n_nuc : int
         Number of nuclides.
-    rates : list of ReactionRates
-        The reaction rates for each substep.
+    rates : ReactionRates
+        The reaction rates at end of step.
     volume : dict of str to float
         Dictionary mapping mat id to volume.
     index_mat : dict of str to int
@@ -52,10 +52,8 @@ class StepResult:
         A dictionary mapping mat ID as string to global index.
     n_hdf5_mats : int
         Number of materials in entire geometry.
-    n_stages : int
-        Number of stages in simulation.
     data : numpy.ndarray
-        Atom quantity, stored by stage, mat, then by nuclide.
+        Atom quantity, stored by mat, then by nuclide.
     proc_time : int
         Average time spent depleting a material across all
         materials and processes
@@ -86,17 +84,17 @@ def __getitem__(self, pos):
         Parameters
         ----------
         pos : tuple
-            A three-length tuple containing a stage index, mat index and a nuc
-            index.  All can be integers or slices.  The second two can be
+            A two-length tuple containing a mat index and a nuc
+            index.  Both can be integers or slices, or can be
             strings corresponding to their respective dictionary.
 
         Returns
         -------
         float
-            The atoms for stage, mat, nuc
+            The atoms for mat, nuc
 
         """
-        stage, mat, nuc = pos
+        mat, nuc = pos
         if isinstance(mat, openmc.Material):
             mat = str(mat.id)
         if isinstance(mat, str):
@@ -104,7 +102,7 @@ def __getitem__(self, pos):
         if isinstance(nuc, str):
             nuc = self.index_nuc[nuc]
 
-        return self.data[stage, mat, nuc]
+        return self.data[mat, nuc]
 
     def __setitem__(self, pos, val):
         """Sets an item from results.
@@ -112,21 +110,21 @@ def __setitem__(self, pos, val):
         Parameters
         ----------
         pos : tuple
-            A three-length tuple containing a stage index, mat index and a nuc
-            index.  All can be integers or slices.  The second two can be
+            A two-length tuple containing a mat index and a nuc
+            index.  Both can be integers or slices, or can be
             strings corresponding to their respective dictionary.
 
         val : float
             The value to set data to.
 
         """
-        stage, mat, nuc = pos
+        mat, nuc = pos
         if isinstance(mat, str):
             mat = self.index_mat[mat]
         if isinstance(nuc, str):
             nuc = self.index_nuc[nuc]
 
-        self.data[stage, mat, nuc] = val
+        self.data[mat, nuc] = val
 
     @property
     def n_mat(self):
@@ -140,11 +138,7 @@ def n_nuc(self):
     def n_hdf5_mats(self):
         return len(self.mat_to_hdf5_ind)
 
-    @property
-    def n_stages(self):
-        return self.data.shape[0]
-
-    def allocate(self, volume, nuc_list, burn_list, full_burn_list, stages):
+    def allocate(self, volume, nuc_list, burn_list, full_burn_list):
         """Allocate memory for depletion step data
 
         Parameters
@@ -157,8 +151,6 @@ def allocate(self, volume, nuc_list, burn_list, full_burn_list, stages):
             A list of all mat IDs to be burned.  Used for sorting the simulation.
         full_burn_list : list of str
             List of all burnable material IDs
-        stages : int
-            Number of stages in simulation.
 
         """
         self.volume = copy.deepcopy(volume)
@@ -167,7 +159,7 @@ def allocate(self, volume, nuc_list, burn_list, full_burn_list, stages):
         self.mat_to_hdf5_ind = {mat: i for i, mat in enumerate(full_burn_list)}
 
         # Create storage array
-        self.data = np.zeros((stages, self.n_mat, self.n_nuc))
+        self.data = np.zeros((self.n_mat, self.n_nuc))
 
     def distribute(self, local_materials, ranges):
         """Create a new object containing data for distributed materials
@@ -196,8 +188,8 @@ def distribute(self, local_materials, ranges):
         for attr in direct_attrs:
             setattr(new, attr, getattr(self, attr))
         # Get applicable slice of data
-        new.data = self.data[:, ranges]
-        new.rates = [r[ranges] for r in self.rates]
+        new.data = self.data[ranges]
+        new.rates = self.rates[ranges]
         return new
 
     def get_material(self, mat_id):
@@ -232,7 +224,7 @@ def get_material(self, mat_id):
                 f'values are {list(self.volume.keys())}'
             ) from e
         for nuc, _ in sorted(self.index_nuc.items(), key=lambda x: x[1]):
-            atoms = self[0, mat_id, nuc]
+            atoms = self[mat_id, nuc]
             if atoms <= 0.0:
                 continue
             atom_per_bcm = atoms / vol * 1e-24
@@ -240,7 +232,7 @@ def get_material(self, mat_id):
         material.volume = vol
         return material
 
-    def export_to_hdf5(self, filename, step):
+    def export_to_hdf5(self, filename, step, write_rates: bool = False):
         """Export results to an HDF5 file
 
         Parameters
@@ -249,6 +241,8 @@ def export_to_hdf5(self, filename, step):
             The filename to write to
         step : int
             What step is this?
+        write_rates : bool, optional
+            Whether to include reaction rate datasets in the results file.
 
         """
         # Write new file if first time step, else add to existing file
@@ -259,7 +253,8 @@ def export_to_hdf5(self, filename, step):
             kwargs['driver'] = 'mpio'
             kwargs['comm'] = comm
             with h5py.File(filename, **kwargs) as handle:
-                self._to_hdf5(handle, step, parallel=True)
+                self._to_hdf5(handle, step, parallel=True,
+                              write_rates=write_rates)
         else:
             # Gather results at root process
             all_results = comm.gather(self)
@@ -268,15 +263,18 @@ def export_to_hdf5(self, filename, step):
             if comm.rank == 0:
                 with h5py.File(filename, **kwargs) as handle:
                     for res in all_results:
-                        res._to_hdf5(handle, step, parallel=False)
+                        res._to_hdf5(handle, step, parallel=False,
+                                     write_rates=write_rates)
 
-    def _write_hdf5_metadata(self, handle):
+    def _write_hdf5_metadata(self, handle, write_rates):
         """Writes result metadata in HDF5 file
 
         Parameters
         ----------
         handle : h5py.File or h5py.Group
             An hdf5 file or group type to store this in.
+        write_rates : bool
+            Whether reaction rate datasets are being written.
 
         """
         # Create and save the 5 dictionaries:
@@ -284,8 +282,8 @@ def _write_hdf5_metadata(self, handle):
         #   self.index_mat -> self.volume (TODO: support for changing volumes)
         #   self.index_nuc
         # reactions
-        #   self.rates[0].index_nuc (can be different from above, above is superset)
-        #   self.rates[0].index_rx
+        #   self.rates.index_nuc (can be different from above, above is superset)
+        #   self.rates.index_rx
         # these are shared by every step of the simulation, and should be deduplicated.
 
         # Store concentration mat and nuclide dictionaries (along with volumes)
@@ -295,13 +293,19 @@ def _write_hdf5_metadata(self, handle):
 
         mat_list = sorted(self.mat_to_hdf5_ind, key=int)
         nuc_list = sorted(self.index_nuc)
-        rxn_list = sorted(self.rates[0].index_rx)
+
+        include_rates = (
+            write_rates
+            and self.rates is not None
+            and bool(self.rates.index_nuc)
+            and bool(self.rates.index_rx)
+        )
+        rxn_list = sorted(self.rates.index_rx) if include_rates else []
 
         n_mats = self.n_hdf5_mats
         n_nuc_number = len(nuc_list)
-        n_nuc_rxn = len(self.rates[0].index_nuc)
+        n_nuc_rxn = len(self.rates.index_nuc) if include_rates else 0
         n_rxn = len(rxn_list)
-        n_stages = self.n_stages
 
         mat_group = handle.create_group("materials")
 
@@ -315,41 +319,44 @@ def _write_hdf5_metadata(self, handle):
         for nuc in nuc_list:
             nuc_single_group = nuc_group.create_group(nuc)
             nuc_single_group.attrs["atom number index"] = self.index_nuc[nuc]
-            if nuc in self.rates[0].index_nuc:
-                nuc_single_group.attrs["reaction rate index"] = self.rates[0].index_nuc[nuc]
+            if include_rates and nuc in self.rates.index_nuc:
+                nuc_single_group.attrs["reaction rate index"] = (
+                    self.rates.index_nuc[nuc])
 
-        rxn_group = handle.create_group("reactions")
+        if include_rates:
+            rxn_group = handle.create_group("reactions")
 
-        for rxn in rxn_list:
-            rxn_single_group = rxn_group.create_group(rxn)
-            rxn_single_group.attrs["index"] = self.rates[0].index_rx[rxn]
+            for rxn in rxn_list:
+                rxn_single_group = rxn_group.create_group(rxn)
+                rxn_single_group.attrs["index"] = (
+                    self.rates.index_rx[rxn])
 
         # Construct array storage
 
-        handle.create_dataset("number", (1, n_stages, n_mats, n_nuc_number),
-                              maxshape=(None, n_stages, n_mats, n_nuc_number),
+        handle.create_dataset("number", (1, n_mats, n_nuc_number),
+                              maxshape=(None, n_mats, n_nuc_number),
                               chunks=True,
                               dtype='float64')
 
-        if n_nuc_rxn > 0 and n_rxn > 0:
-            handle.create_dataset("reaction rates", (1, n_stages, n_mats, n_nuc_rxn, n_rxn),
-                                maxshape=(None, n_stages, n_mats, n_nuc_rxn, n_rxn),
-                                chunks=True,
-                                dtype='float64')
+        if include_rates and n_nuc_rxn > 0 and n_rxn > 0:
+            handle.create_dataset(
+                "reaction rates", (1, n_mats, n_nuc_rxn, n_rxn),
+                maxshape=(None, n_mats, n_nuc_rxn, n_rxn),
+                chunks=True, dtype='float64')
 
-        handle.create_dataset("eigenvalues", (1, n_stages, 2),
-                              maxshape=(None, n_stages, 2), dtype='float64')
+        handle.create_dataset("eigenvalues", (1, 2),
+                              maxshape=(None, 2), dtype='float64')
 
         handle.create_dataset("time", (1, 2), maxshape=(None, 2), dtype='float64')
 
-        handle.create_dataset("source_rate", (1, n_stages), maxshape=(None, n_stages),
+        handle.create_dataset("source_rate", (1,), maxshape=(None,),
                               dtype='float64')
 
         handle.create_dataset(
             "depletion time", (1,), maxshape=(None,),
             dtype="float64")
 
-    def _to_hdf5(self, handle, index, parallel=False):
+    def _to_hdf5(self, handle, index, parallel=False, write_rates: bool = False):
         """Converts results object into an hdf5 object.
 
         Parameters
@@ -360,12 +367,14 @@ def _to_hdf5(self, handle, index, parallel=False):
             What step is this?
         parallel : bool
             Being called with parallel HDF5?
+        write_rates : bool, optional
+            Whether reaction rate datasets are being written.
 
         """
         if "/number" not in handle:
             if parallel:
                 comm.barrier()
-            self._write_hdf5_metadata(handle)
+            self._write_hdf5_metadata(handle, write_rates)
 
         if parallel:
             comm.barrier()
@@ -417,18 +426,14 @@ def _to_hdf5(self, handle, index, parallel=False):
             return
 
         # Add data
-        # Note, for the last step, self.n_stages = 1, even if n_stages != 1.
-        n_stages = self.n_stages
         inds = [self.mat_to_hdf5_ind[mat] for mat in self.index_mat]
         low = min(inds)
         high = max(inds)
-        for i in range(n_stages):
-            number_dset[index, i, low:high+1] = self.data[i]
-            if has_reactions:
-                rxn_dset[index, i, low:high+1] = self.rates[i]
-            if comm.rank == 0:
-                eigenvalues_dset[index, i] = self.k[i]
+        number_dset[index, low:high+1] = self.data
+        if has_reactions:
+            rxn_dset[index, low:high+1] = self.rates
         if comm.rank == 0:
+            eigenvalues_dset[index] = self.k
             time_dset[index] = self.time
             source_rate_dset[index] = self.source_rate
             if self.proc_time is not None:
@@ -459,10 +464,24 @@ def from_hdf5(cls, handle, step):
             # Older versions used "power" instead of "source_rate"
             source_rate_dset = handle["/power"]
 
-        results.data = number_dset[step, :, :, :]
-        results.k = eigenvalues_dset[step, :]
+        # Check if this is an old format file (with stages dimension) or new format
+        # Old format: number has shape (n_steps, n_stages, n_mats, n_nucs)
+        # New format: number has shape (n_steps, n_mats, n_nucs)
+        has_stages = len(number_dset.shape) == 4
+
+        if has_stages:
+            # Old format - extract data from first stage (index 0)
+            results.data = number_dset[step, 0, :, :]
+            results.k = eigenvalues_dset[step, 0, :]
+            # source_rate had shape (n_steps, n_stages) in old format
+            results.source_rate = source_rate_dset[step, 0]
+        else:
+            # New format - no stages dimension
+            results.data = number_dset[step, :, :]
+            results.k = eigenvalues_dset[step, :]
+            results.source_rate = source_rate_dset[step]
+
         results.time = time_dset[step, :]
-        results.source_rate = source_rate_dset[step, 0]
 
         if "depletion time" in handle:
             proc_time_dset = handle["/depletion time"]
@@ -493,33 +512,45 @@ def from_hdf5(cls, handle, step):
             if "reaction rate index" in nuc_handle.attrs:
                 rxn_nuc_to_ind[nuc] = nuc_handle.attrs["reaction rate index"]
 
-        for rxn, rxn_handle in handle["/reactions"].items():
-            rxn_to_ind[rxn] = rxn_handle.attrs["index"]
-
-        results.rates = []
-        # Reconstruct reactions
-        for i in range(results.n_stages):
-            rate = ReactionRates(results.index_mat, rxn_nuc_to_ind, rxn_to_ind, True)
+        if "reactions" in handle:
+            for rxn, rxn_handle in handle["/reactions"].items():
+                rxn_to_ind[rxn] = rxn_handle.attrs["index"]
 
-            if "reaction rates" in handle:
-                rate[:] = handle["/reaction rates"][step, i, :, :, :]
-            results.rates.append(rate)
+        # Reconstruct reaction rates
+        rate = ReactionRates(results.index_mat, rxn_nuc_to_ind, rxn_to_ind, True)
+        if "reaction rates" in handle:
+            if has_stages:
+                # Old format: (n_steps, n_stages, n_mats, n_nucs, n_rxns)
+                rate[:] = handle["/reaction rates"][step, 0, :, :, :]
+            else:
+                # New format: (n_steps, n_mats, n_nucs, n_rxns)
+                rate[:] = handle["/reaction rates"][step, :, :, :]
+        results.rates = rate
 
         return results
 
     @staticmethod
-    def save(op, x, op_results, t, source_rate, step_ind, proc_time=None,
-             path: PathLike = "depletion_results.h5"):
+    def save(
+        op,
+        x,
+        op_results,
+        t,
+        source_rate,
+        step_ind,
+        proc_time=None,
+        write_rates: bool = False,
+        path: PathLike = "depletion_results.h5"
+    ):
         """Creates and writes depletion results to disk
 
         Parameters
         ----------
         op : openmc.deplete.abc.TransportOperator
             The operator used to generate these results.
-        x : list of list of numpy.array
-            The prior x vectors.  Indexed [i][cell] using the above equation.
-        op_results : list of openmc.deplete.OperatorResult
-            Results of applying transport operator
+        x : numpy.array
+            End-of-step concentrations for each material
+        op_results : openmc.deplete.OperatorResult
+            Result of applying transport operator at end of step
         t : list of float
             Time indices.
         source_rate : float
@@ -530,7 +561,8 @@ def save(op, x, op_results, t, source_rate, step_ind, proc_time=None,
             Total process time spent depleting materials. This may
             be process-dependent and will be reduced across MPI
             processes.
-
+        write_rates : bool, optional
+            Whether reaction rates should be written to the results file.
         path : PathLike
             Path to file to write. Defaults to 'depletion_results.h5'.
 
@@ -539,26 +571,20 @@ def save(op, x, op_results, t, source_rate, step_ind, proc_time=None,
         # Get indexing terms
         vol_dict, nuc_list, burn_list, full_burn_list = op.get_results_info()
 
-        stages = len(x)
-
         # Create results
         results = StepResult()
-        results.allocate(vol_dict, nuc_list, burn_list, full_burn_list, stages)
+        results.allocate(vol_dict, nuc_list, burn_list, full_burn_list)
 
         n_mat = len(burn_list)
 
-        for i in range(stages):
-            for mat_i in range(n_mat):
-                results[i, mat_i, :] = x[i][mat_i]
+        for mat_i in range(n_mat):
+            results[mat_i, :] = x[mat_i]
 
-        ks = []
-        for r in op_results:
-            if isinstance(r.k, type(None)):
-                ks += [(None, None)]
-            else:
-                ks += [(r.k.nominal_value, r.k.std_dev)]
-        results.k = ks
-        results.rates = [r.rates for r in op_results]
+        if isinstance(op_results.k, type(None)):
+            results.k = (None, None)
+        else:
+            results.k = (op_results.k.nominal_value, op_results.k.std_dev)
+        results.rates = op_results.rates
         results.time = t
         results.source_rate = source_rate
         results.proc_time = proc_time
@@ -567,7 +593,7 @@ def save(op, x, op_results, t, source_rate, step_ind, proc_time=None,
 
         if not Path(path).is_file():
             Path(path).parent.mkdir(parents=True, exist_ok=True)
-        results.export_to_hdf5(path, step_ind)
+        results.export_to_hdf5(path, step_ind, write_rates)
 
     def transfer_volumes(self, model):
         """Transfers volumes from depletion results to geometry
diff --git a/tests/unit_tests/test_deplete_continue.py b/tests/unit_tests/test_deplete_continue.py
index 637c9d5e440..1b6eac2384f 100644
--- a/tests/unit_tests/test_deplete_continue.py
+++ b/tests/unit_tests/test_deplete_continue.py
@@ -17,7 +17,7 @@ def test_continue(run_in_tmpdir):
     operator = dummy_operator.DummyOperator()
 
     # initial depletion
-    bundle.solver(operator, [1.0, 2.0], [1.0, 2.0]).integrate()
+    bundle.solver(operator, [1.0, 2.0], [1.0, 2.0]).integrate(write_rates=True)
 
     # set up continue run
     prev_res = openmc.deplete.Results(operator.output_dir / "depletion_results.h5")
@@ -25,7 +25,7 @@ def test_continue(run_in_tmpdir):
 
     # if continue run happens, test passes
     bundle.solver(operator, [1.0, 2.0, 3.0, 4.0], [1.0, 2.0, 3.0, 4.0],
-                  continue_timesteps=True).integrate()
+                  continue_timesteps=True).integrate(write_rates=True)
 
     final_res = openmc.deplete.Results(operator.output_dir / "depletion_results.h5")
 
@@ -42,7 +42,7 @@ def test_continue_continue(run_in_tmpdir):
     operator = dummy_operator.DummyOperator()
 
     # initial depletion
-    bundle.solver(operator, [1.0, 2.0], [1.0, 2.0]).integrate()
+    bundle.solver(operator, [1.0, 2.0], [1.0, 2.0]).integrate(write_rates=True)
 
     # set up continue run
     prev_res = openmc.deplete.Results(operator.output_dir / "depletion_results.h5")
@@ -50,7 +50,7 @@ def test_continue_continue(run_in_tmpdir):
 
     # first continue run
     bundle.solver(operator, [1.0, 2.0, 3.0, 4.0], [1.0, 2.0, 3.0, 4.0],
-                  continue_timesteps=True).integrate()
+                  continue_timesteps=True).integrate(write_rates=True)
 
     prev_res = openmc.deplete.Results(operator.output_dir / "depletion_results.h5")
     # second continue run
diff --git a/tests/unit_tests/test_deplete_integrator.py b/tests/unit_tests/test_deplete_integrator.py
index b1d2cb950eb..6463eaa2084 100644
--- a/tests/unit_tests/test_deplete_integrator.py
+++ b/tests/unit_tests/test_deplete_integrator.py
@@ -10,6 +10,7 @@
 from random import uniform
 from unittest.mock import MagicMock
 
+import h5py
 import numpy as np
 from uncertainties import ufloat
 import pytest
@@ -38,8 +39,6 @@
 def test_results_save(run_in_tmpdir):
     """Test data save module"""
 
-    stages = 3
-
     rng = np.random.RandomState(comm.rank)
 
     # Mock geometry
@@ -63,31 +62,22 @@ def test_results_save(run_in_tmpdir):
     op.get_results_info.return_value = (
         vol_dict, nuc_list, burn_list, full_burn_list)
 
-    # Construct x
-    x1 = []
-    x2 = []
-
-    for i in range(stages):
-        x1.append([rng.random(2), rng.random(2)])
-        x2.append([rng.random(2), rng.random(2)])
+    # Construct end-of-step concentrations
+    x1 = [rng.random(2), rng.random(2)]
+    x2 = [rng.random(2), rng.random(2)]
 
-    # Construct r
+    # Construct reaction rates
     r1 = ReactionRates(burn_list, ["na", "nb"], ["ra", "rb"])
     r1[:] = rng.random((2, 2, 2))
+    rate1 = copy.deepcopy(r1)
 
-    rate1 = []
-    rate2 = []
-
-    for i in range(stages):
-        rate1.append(copy.deepcopy(r1))
-        r1[:] = rng.random((2, 2, 2))
-        rate2.append(copy.deepcopy(r1))
-        r1[:] = rng.random((2, 2, 2))
+    r2 = ReactionRates(burn_list, ["na", "nb"], ["ra", "rb"])
+    r2[:] = rng.random((2, 2, 2))
+    rate2 = copy.deepcopy(r2)
 
-    # Create global terms
-    # Col 0: eig, Col 1: uncertainty
-    eigvl1 = rng.random((stages, 2))
-    eigvl2 = rng.random((stages, 2))
+    # Create global terms (eigenvalue and uncertainty)
+    eigvl1 = rng.random(2)
+    eigvl2 = rng.random(2)
 
     eigvl1 = comm.bcast(eigvl1, root=0)
     eigvl2 = comm.bcast(eigvl2, root=0)
@@ -95,29 +85,35 @@ def test_results_save(run_in_tmpdir):
     t1 = [0.0, 1.0]
     t2 = [1.0, 2.0]
 
-    op_result1 = [OperatorResult(ufloat(*k), rates)
-                  for k, rates in zip(eigvl1, rate1)]
-    op_result2 = [OperatorResult(ufloat(*k), rates)
-                  for k, rates in zip(eigvl2, rate2)]
+    op_result1 = OperatorResult(ufloat(*eigvl1), rate1)
+    op_result2 = OperatorResult(ufloat(*eigvl2), rate2)
 
     # saves within a subdirectory
-    StepResult.save(op, x1, op_result1, t1, 0, 0, path='out/put/depletion.h5')
+    StepResult.save(
+        op,
+        x1,
+        op_result1,
+        t1,
+        0,
+        0,
+        write_rates=True,
+        path='out/put/depletion.h5'
+    )
     res = Results('out/put/depletion.h5')
 
     # saves with default filename
-    StepResult.save(op, x1, op_result1, t1, 0, 0)
-    StepResult.save(op, x2, op_result2, t2, 0, 1)
+    StepResult.save(op, x1, op_result1, t1, 0, 0, write_rates=True)
+    StepResult.save(op, x2, op_result2, t2, 0, 1, write_rates=True)
 
     # Load the files
     res = Results("depletion_results.h5")
 
-    for i in range(stages):
-        for mat_i, mat in enumerate(burn_list):
-            for nuc_i, nuc in enumerate(nuc_list):
-                assert res[0][i, mat, nuc] == x1[i][mat_i][nuc_i]
-                assert res[1][i, mat, nuc] == x2[i][mat_i][nuc_i]
-        np.testing.assert_array_equal(res[0].rates[i], rate1[i])
-        np.testing.assert_array_equal(res[1].rates[i], rate2[i])
+    for mat_i, mat in enumerate(burn_list):
+        for nuc_i, nuc in enumerate(nuc_list):
+            assert res[0][mat, nuc] == x1[mat_i][nuc_i]
+            assert res[1][mat, nuc] == x2[mat_i][nuc_i]
+    np.testing.assert_array_equal(res[0].rates, rate1)
+    np.testing.assert_array_equal(res[1].rates, rate2)
 
     np.testing.assert_array_equal(res[0].k, eigvl1)
     np.testing.assert_array_equal(res[0].time, t1)
@@ -126,6 +122,31 @@ def test_results_save(run_in_tmpdir):
     np.testing.assert_array_equal(res[1].time, t2)
 
 
+def test_results_save_without_rates(run_in_tmpdir):
+    """StepResult.save skips reaction-rate datasets by default"""
+
+    op = MagicMock()
+    op.prev_res = None
+    vol_dict = {"0": 1.0}
+    nuc_list = ["na"]
+    burn_list = ["0"]
+    op.get_results_info.return_value = (vol_dict, nuc_list, burn_list, burn_list)
+
+    x = [np.array([1.0])]
+    rates = ReactionRates(burn_list, nuc_list, ["ra"])
+    rates[:] = np.array([[[2.0]]])
+    op_result = OperatorResult(ufloat(1.0, 0.1), rates)
+
+    StepResult.save(op, x, op_result, [0.0, 1.0], 0.0, 0)
+
+    with h5py.File('depletion_results.h5', 'r') as handle:
+        assert 'reaction rates' not in handle
+        assert 'reactions' not in handle
+
+    res = Results('depletion_results.h5')
+    assert res[0].rates.size == 0
+
+
 def test_bad_integrator_inputs():
     """Test failure modes for Integrator inputs"""
 
diff --git a/tests/unit_tests/test_deplete_restart.py b/tests/unit_tests/test_deplete_restart.py
index e8bfc062a08..1cbff30a5fb 100644
--- a/tests/unit_tests/test_deplete_restart.py
+++ b/tests/unit_tests/test_deplete_restart.py
@@ -21,7 +21,7 @@ def test_restart_predictor_cecm(run_in_tmpdir):
     # Perform simulation using the predictor algorithm
     dt = [0.75]
     power = 1.0
-    openmc.deplete.PredictorIntegrator(op, dt, power).integrate()
+    openmc.deplete.PredictorIntegrator(op, dt, power).integrate(write_rates=True)
 
     # Load the files
     prev_res = openmc.deplete.Results(op.output_dir / "depletion_results.h5")
@@ -30,10 +30,6 @@ def test_restart_predictor_cecm(run_in_tmpdir):
     op = dummy_operator.DummyOperator(prev_res)
     op.output_dir = output_dir
 
-    # check ValueError is raised, indicating previous and current stages
-    with pytest.raises(ValueError, match="incompatible.* 1.*2"):
-        openmc.deplete.CECMIntegrator(op, dt, power)
-
 
 def test_restart_cecm_predictor(run_in_tmpdir):
     """Integral regression test of integrator algorithm using CE/CM for the
@@ -47,7 +43,7 @@ def test_restart_cecm_predictor(run_in_tmpdir):
     dt = [0.75]
     power = 1.0
     cecm = openmc.deplete.CECMIntegrator(op, dt, power)
-    cecm.integrate()
+    cecm.integrate(write_rates=True)
 
     # Load the files
     prev_res = openmc.deplete.Results(op.output_dir / "depletion_results.h5")
@@ -56,10 +52,6 @@ def test_restart_cecm_predictor(run_in_tmpdir):
     op = dummy_operator.DummyOperator(prev_res)
     op.output_dir = output_dir
 
-    # check ValueError is raised, indicating previous and current stages
-    with pytest.raises(ValueError, match="incompatible.* 2.*1"):
-        openmc.deplete.PredictorIntegrator(op, dt, power)
-
 
 @pytest.mark.parametrize("scheme", dummy_operator.SCHEMES)
 def test_restart(run_in_tmpdir, scheme):
@@ -70,7 +62,7 @@ def test_restart(run_in_tmpdir, scheme):
     operator = dummy_operator.DummyOperator()
 
     # take first step
-    bundle.solver(operator, [0.75], 1.0).integrate()
+    bundle.solver(operator, [0.75], 1.0).integrate(write_rates=True)
 
     # restart
     prev_res = openmc.deplete.Results(