diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h
index 5b632bbce53..a8781a0bccf 100644
--- a/include/openmc/particle_data.h
+++ b/include/openmc/particle_data.h
@@ -210,6 +210,8 @@ struct MacroXS {
   double incoherent;      //!< macroscopic incoherent xs
   double photoelectric;   //!< macroscopic photoelectric xs
   double pair_production; //!< macroscopic pair production xs
+
+  double last_E {0.0}; //!< last evaluated energy in [eV]
 };
 
 //==============================================================================
diff --git a/openmc/data/decay.py b/openmc/data/decay.py
index 7cd4bf43d41..ef616384366 100644
--- a/openmc/data/decay.py
+++ b/openmc/data/decay.py
@@ -564,7 +564,7 @@ def sources(self):
 _DECAY_PHOTON_ENERGY = {}
 
 
-def decay_photon_energy(nuclide: str) -> Univariate | None:
+def decay_photon_energy(nuclide: str, chain = None) -> Univariate | None:
     """Get photon energy distribution resulting from the decay of a nuclide
 
     This function relies on data stored in a depletion chain. Before calling it
@@ -577,7 +577,9 @@ def decay_photon_energy(nuclide: str) -> Univariate | None:
     ----------
     nuclide : str
         Name of nuclide, e.g., 'Co58'
-
+    chain : Chain, optional.
+            Chain file to get decay energy from
+    
     Returns
     -------
     openmc.stats.Univariate or None
@@ -585,6 +587,15 @@ def decay_photon_energy(nuclide: str) -> Univariate | None:
         if no photon source exists. Note that the probabilities represent
         intensities, given as [Bq/atom] (in other words, decay constants).
     """
+    from openmc.deplete.chain import Chain
+    if chain is not None:
+        cv.check_type('chain', chain, Chain)
+        for nuc in chain.nuclides:
+            if nuc.name == nuclide:
+                return nuc.sources.get('photon')
+        else:
+            return
+                                                            
     if not _DECAY_PHOTON_ENERGY:
         chain_file = openmc.config.get('chain_file')
         if chain_file is None:
@@ -593,7 +604,6 @@ def decay_photon_energy(nuclide: str) -> Univariate | None:
                 "openmc.config['chain_file'] in order to load decay data."
             )
 
-        from openmc.deplete import Chain
         chain = Chain.from_xml(chain_file)
         for nuc in chain.nuclides:
             if 'photon' in nuc.sources:
@@ -610,7 +620,7 @@ def decay_photon_energy(nuclide: str) -> Univariate | None:
 _DECAY_ENERGY = {}
 
 
-def decay_energy(nuclide: str):
+def decay_energy(nuclide: str, chain = None):
     """Get decay energy value resulting from the decay of a nuclide
 
     This function relies on data stored in a depletion chain. Before calling it
@@ -623,6 +633,8 @@ def decay_energy(nuclide: str):
     ----------
     nuclide : str
         Name of nuclide, e.g., 'H3'
+    chain : Chain, optional.
+        Chain file to get decay energy from
 
     Returns
     -------
@@ -630,6 +642,15 @@ def decay_energy(nuclide: str):
         Decay energy of nuclide in [eV]. If the nuclide is stable, a value of
         0.0 is returned.
     """
+    from openmc.deplete.chain import Chain
+    if chain is not None:
+        cv.check_type('chain', chain, Chain)
+        for nuclide in chain.nuclides:
+            if nuc.name == nuclide and nuc.decay_energy:
+                return nuc.decay_energy
+        else:
+            return 0.0
+
     if not _DECAY_ENERGY:
         chain_file = openmc.config.get('chain_file')
         if chain_file is None:
@@ -638,7 +659,6 @@ def decay_energy(nuclide: str):
                 "openmc.config['chain_file'] in order to load decay data."
             )
 
-        from openmc.deplete import Chain
         chain = Chain.from_xml(chain_file)
         for nuc in chain.nuclides:
             if nuc.decay_energy:
diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py
index 57bbe437ffb..9283a128836 100644
--- a/openmc/deplete/r2s.py
+++ b/openmc/deplete/r2s.py
@@ -230,7 +230,7 @@ def run(
         )
         self.step3_photon_transport(
             photon_time_indices, bounding_boxes, output_dir / 'photon_transport',
-            mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs
+            mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs, chain_file = chain_file
         )
 
         return output_dir
@@ -405,6 +405,7 @@ def step3_photon_transport(
         output_dir: PathLike = 'photon_transport',
         mat_vol_kwargs: dict | None = None,
         run_kwargs: dict | None = None,
+        chain_file: PathLike | None = None,
     ):
         """Run the photon transport step.
 
@@ -437,12 +438,17 @@ def step3_photon_transport(
         run_kwargs : dict, optional
             Additional keyword arguments passed to :meth:`openmc.Model.run`
             during the photon transport step. By default, output is disabled.
+        chain_file : PathLike, optional
+                    Path to the depletion chain XML file to use during activation. If
+                    not provided, the default configured chain file will be used.
         """
 
         # TODO: Automatically determine bounding box for each cell
         if bounding_boxes is None and self.method == 'cell-based':
             raise ValueError("bounding_boxes must be provided for cell-based "
                              "R2S calculations.")
+        from openmc.deplete import Chain
+        chain = Chain.from_xml(chain_file) if chain_file is not None else None
 
         # Set default run arguments if not provided
         if run_kwargs is None:
@@ -490,7 +496,7 @@ def step3_photon_transport(
             # Create decay photon source
             if self.method == 'mesh-based':
                 self.photon_model.settings.source = \
-                    self.get_decay_photon_source_mesh(time_index)
+                    self.get_decay_photon_source_mesh(time_index, chain_file=chain_file)
             else:
                 sources = []
                 results = self.results['depletion_results']
@@ -510,7 +516,7 @@ def step3_photon_transport(
 
                     # Create decay photon source source
                     space = openmc.stats.Box(*bounding_box)
-                    energy = activated_mat.get_decay_photon_energy()
+                    energy = activated_mat.get_decay_photon_energy(chain=chain)
                     strength = energy.integral() if energy is not None else 0.0
                     source = openmc.IndependentSource(
                         space=space,
@@ -539,7 +545,8 @@ def step3_photon_transport(
 
     def get_decay_photon_source_mesh(
         self,
-        time_index: int = -1
+        time_index: int = -1,
+        chain_file: PathLike | None = None,
     ) -> list[openmc.IndependentSource]:
         """Create decay photon source for a mesh-based calculation.
 
@@ -558,6 +565,9 @@ def get_decay_photon_source_mesh(
         ----------
         time_index : int, optional
             Time index for the decay photon source. Default is -1 (last time).
+        chain_file : PathLike, optional
+                    Path to the depletion chain XML file to use during activation. If
+                    not provided, the default configured chain file will be used.
 
         Returns
         -------
@@ -566,6 +576,9 @@ def get_decay_photon_source_mesh(
             each mesh element-material combination with non-zero source strength.
 
         """
+        from openmc.deplete import Chain
+        chain = Chain.from_xml(chain_file) if chain_file is not None else None
+                        
         mat_dict = self.neutron_model._get_all_materials()
 
         # List to hold all sources
@@ -607,7 +620,7 @@ def get_decay_photon_source_mesh(
                 activated_mat = results[time_index].get_material(str(original_mat.id))
 
                 # Create decay photon source
-                energy = activated_mat.get_decay_photon_energy()
+                energy = activated_mat.get_decay_photon_energy(chain=chain)
                 if energy is not None:
                     strength = energy.integral()
                     space = openmc.stats.Box(*bbox)
diff --git a/openmc/material.py b/openmc/material.py
index 239bc01f771..7b5fb8ceada 100644
--- a/openmc/material.py
+++ b/openmc/material.py
@@ -28,7 +28,6 @@
 from openmc.data.function import Tabulated1D
 from openmc.data import mass_energy_absorption_coefficient, dose_coefficients
 
-
 # Units for density supported by OpenMC
 DENSITY_UNITS = ('g/cm3', 'g/cc', 'kg/m3', 'atom/b-cm', 'atom/cm3', 'sum',
                  'macro')
@@ -338,7 +337,8 @@ def get_decay_photon_energy(
         units: str = 'Bq',
         volume: float | None = None,
         exclude_nuclides: list[str] | None = None,
-        include_nuclides: list[str] | None = None
+        include_nuclides: list[str] | None = None,
+        chain: Chain | None = None,
     ) -> Univariate | None:
         r"""Return energy distribution of decay photons from unstable nuclides.
 
@@ -359,6 +359,8 @@ def get_decay_photon_energy(
         include_nuclides : list of str, optional
             Nuclides to include in the photon source calculation. If specified,
             only these nuclides are used.
+        chain : Chain, optional.
+            Chain file to get decay energy from
 
         Returns
         -------
@@ -367,7 +369,9 @@ def get_decay_photon_energy(
             the total intensity of the photon source in the requested units.
 
         """
+        from openmc.deplete import Chain
         cv.check_value('units', units, {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3'})
+        cv.check_type('chain', chain, Chain, none_ok=True)
 
         if exclude_nuclides is not None and include_nuclides is not None:
             raise ValueError("Cannot specify both exclude_nuclides and include_nuclides")
@@ -391,7 +395,7 @@ def get_decay_photon_energy(
             if include_nuclides is not None and nuc not in include_nuclides:
                 continue
 
-            source_per_atom = openmc.data.decay_photon_energy(nuc)
+            source_per_atom = openmc.data.decay_photon_energy(nuc, chain=chain)
             if source_per_atom is not None and atoms_per_bcm > 0.0:
                 dists.append(source_per_atom)
                 probs.append(1e24 * atoms_per_bcm * multiplier)
diff --git a/src/geometry.cpp b/src/geometry.cpp
index ddb61385f18..e8e8a065112 100644
--- a/src/geometry.cpp
+++ b/src/geometry.cpp
@@ -173,12 +173,21 @@ bool find_cell_inner(
       }
 
       // Set the material, temperature and density multiplier.
-      p.material_last() = p.material();
-      p.material() = c.material(p.cell_instance());
-      p.sqrtkT_last() = p.sqrtkT();
-      p.sqrtkT() = c.sqrtkT(p.cell_instance());
-      p.density_mult_last() = p.density_mult();
-      p.density_mult() = c.density_mult(p.cell_instance());
+      if (p.sqrtkT() != -1) {
+        // if current material is set, set last material
+        p.material_last() = p.material();
+        p.material() = c.material(p.cell_instance());
+        p.sqrtkT_last() = p.sqrtkT();
+        p.sqrtkT() = c.sqrtkT(p.cell_instance());
+        p.density_mult_last() = p.density_mult();
+        p.density_mult() = c.density_mult(p.cell_instance());
+      } else {
+        // set both current data and last data to current data
+        p.material_last() = p.material() = c.material(p.cell_instance());
+        p.sqrtkT_last() = p.sqrtkT() = c.sqrtkT(p.cell_instance());
+        p.density_mult_last() = p.density_mult() =
+          c.density_mult(p.cell_instance());
+      }
 
       return true;
 
diff --git a/src/material.cpp b/src/material.cpp
index 21b11b8b9ce..16ad78cf12c 100644
--- a/src/material.cpp
+++ b/src/material.cpp
@@ -817,6 +817,7 @@ void Material::calculate_xs(Particle& p) const
   p.macro_xs().absorption = 0.0;
   p.macro_xs().fission = 0.0;
   p.macro_xs().nu_fission = 0.0;
+  p.macro_xs().last_E = p.E();
 
   if (p.type().is_neutron()) {
     this->calculate_neutron_xs(p);
diff --git a/src/mgxs.cpp b/src/mgxs.cpp
index a0fe4060ddb..640c990c634 100644
--- a/src/mgxs.cpp
+++ b/src/mgxs.cpp
@@ -615,6 +615,7 @@ void Mgxs::calculate_xs(Particle& p)
   p.macro_xs().nu_fission =
     fissionable ? xs[temperature].nu_fission(angle, p.g()) * p.density_mult()
                 : 0.;
+  p.macro_xs().last_E = p.E();
 }
 
 //==============================================================================
diff --git a/src/particle.cpp b/src/particle.cpp
index 0a063559824..953497cf376 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -198,6 +198,8 @@ void Particle::event_calculate_xs()
       cell_last(j) = coord(j).cell();
     }
     n_coord_last() = n_coord();
+
+    material_last() = material();
   }
 
   // Write particle track.
@@ -211,9 +213,10 @@ void Particle::event_calculate_xs()
   if (material() != MATERIAL_VOID) {
     if (settings::run_CE) {
       if (material() != material_last() || sqrtkT() != sqrtkT_last() ||
-          density_mult() != density_mult_last()) {
+          density_mult() != density_mult_last() || E() != macro_xs().last_E ||
+          true) {
         // If the material is the same as the last material and the
-        // temperature hasn't changed, we don't need to lookup cross
+        // energy and temperature hasn't changed, we don't need to lookup cross
         // sections again.
         model::materials[material()]->calculate_xs(*this);
       }
@@ -231,6 +234,7 @@ void Particle::event_calculate_xs()
     macro_xs().absorption = 0.0;
     macro_xs().fission = 0.0;
     macro_xs().nu_fission = 0.0;
+    macro_xs().last_E = E();
   }
 }
 
@@ -292,10 +296,7 @@ void Particle::event_advance()
 
 void Particle::event_cross_surface()
 {
-  // Saving previous cell data
-  for (int j = 0; j < n_coord(); ++j) {
-    cell_last(j) = coord(j).cell();
-  }
+
   n_coord_last() = n_coord();
 
   // Set surface that particle is on and adjust coordinate levels
@@ -309,6 +310,11 @@ void Particle::event_cross_surface()
       boundary().lattice_translation()[2] != 0) {
     // Particle crosses lattice boundary
 
+    // Saving previous cell data
+    for (int j = 0; j < n_coord(); ++j) {
+      cell_last(j) = coord(j).cell();
+    }
+
     bool verbose = settings::verbosity >= 10 || trace();
     cross_lattice(*this, boundary(), verbose);
     event() = TallyEvent::LATTICE;
@@ -397,10 +403,6 @@ void Particle::event_collide()
   // Save coordinates for tallying purposes
   r_last_current() = r();
 
-  // Set last material to none since cross sections will need to be
-  // re-evaluated
-  material_last() = C_NONE;
-
   // Set all directions to base level -- right now, after a collision, only
   // the base level directions are changed
   for (int j = 0; j < n_coord() - 1; ++j) {
@@ -419,6 +421,15 @@ void Particle::event_collide()
   if (!model::active_tallies.empty())
     score_collision_derivative(*this);
 
+  // Saving previous cell data
+  for (int j = 0; j < n_coord(); ++j) {
+    cell_last(j) = coord(j).cell();
+  }
+  n_coord_last() = n_coord();
+
+  // Saving previous material data
+  material_last() = material();
+
 #ifdef OPENMC_DAGMC_ENABLED
   history().reset();
 #endif
@@ -472,6 +483,8 @@ void Particle::event_revive_from_secondary()
           cell_last(j) = coord(j).cell();
         }
         n_coord_last() = n_coord();
+
+        material_last() = material();
       }
       pht_secondary_particles();
     }
@@ -568,7 +581,6 @@ void Particle::cross_surface(const Surface& surf)
   if (surf.geom_type() == GeometryType::CSG)
     history().reset();
 #endif
-
   // Handle any applicable boundary conditions.
   if (surf.bc_ && settings::run_mode != RunMode::PLOTTING &&
       settings::run_mode != RunMode::VOLUME) {
@@ -576,6 +588,11 @@ void Particle::cross_surface(const Surface& surf)
     return;
   }
 
+  // Saving previous cell data
+  for (int j = 0; j < n_coord(); ++j) {
+    cell_last(j) = coord(j).cell();
+  }
+
   // ==========================================================================
   // SEARCH NEIGHBOR LISTS FOR NEXT CELL
 
@@ -694,10 +711,11 @@ void Particle::cross_reflective_bc(const Surface& surf, Direction new_u)
   // Set the new particle direction
   u() = new_u;
 
-  // Reassign particle's cell and surface
+  // Reassign particle's cell, material and surface
   coord(0).cell() = cell_last(0);
-  surface() = -surface();
+  material() = material_last();
 
+  surface() = -surface();
   // If a reflective surface is coincident with a lattice or universe
   // boundary, it is necessary to redetermine the particle's coordinates in
   // the lower universes.
diff --git a/tests/regression_tests/filter_cellfrom/results_true.dat b/tests/regression_tests/filter_cellfrom/results_true.dat
index 408a4965b1d..1f9d0446104 100644
--- a/tests/regression_tests/filter_cellfrom/results_true.dat
+++ b/tests/regression_tests/filter_cellfrom/results_true.dat
@@ -1,29 +1,29 @@
 k-combined:
 9.035025E-02 2.654309E-03
 tally 1:
-5.994069E+00
-3.594398E+00
+6.586649E+00
+4.340013E+00
 tally 2:
 4.559707E-04
 2.080739E-08
 tally 3:
-5.994525E+00
-3.594945E+00
+6.587105E+00
+4.340613E+00
 tally 4:
 0.000000E+00
 0.000000E+00
 tally 5:
-1.473892E+00
-2.187509E-01
+8.813120E-01
+7.829296E-02
 tally 6:
 5.223875E-05
 2.992861E-10
 tally 7:
-1.473945E+00
-2.187663E-01
+8.813643E-01
+7.830214E-02
 tally 8:
-1.885798E+01
-3.558423E+01
+5.908648E+00
+3.492172E+00
 tally 9:
 7.467961E+00
 5.580255E+00
@@ -34,8 +34,8 @@ tally 11:
 7.468470E+00
 5.581014E+00
 tally 12:
-1.885798E+01
-3.558423E+01
+5.908648E+00
+3.492172E+00
 tally 13:
 0.000000E+00
 0.000000E+00
@@ -46,8 +46,8 @@ tally 15:
 2.739543E-04
 7.600983E-09
 tally 16:
-7.881296E+01
-6.221087E+02
+9.176229E+01
+8.428689E+02
 tally 17:
 1.051397E+02
 1.106292E+03
diff --git a/tests/unit_tests/test_energy_balance.py b/tests/unit_tests/test_energy_balance.py
new file mode 100644
index 00000000000..d1c07b417b4
--- /dev/null
+++ b/tests/unit_tests/test_energy_balance.py
@@ -0,0 +1,152 @@
+import openmc
+import pytest
+import numpy as np
+
+
+@pytest.fixture
+def two_cell_model():
+    """Simple two-cell slab model with a fixed source.
+       Cell1 occupies x in [-10, 0], cell2 x in [0, 10]. 
+    """
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+
+    m1 = openmc.Material()
+    m1.add_element('Li', 1.0)
+    m1.set_density('g/cm3', 1.0)
+
+    m2 = openmc.Material()
+    m2.add_element('Al', 1.0)
+    m2.set_density('g/cm3', 1.0)
+    
+    xmin = openmc.XPlane(-10.0, boundary_type="reflective")
+    xmid = openmc.XPlane(0.0)
+    xmax = openmc.XPlane(10.0, boundary_type="reflective")
+    ymin = openmc.YPlane(-10.0, boundary_type="reflective")
+    ymax = openmc.YPlane(10.0, boundary_type="reflective")
+    zmin = openmc.ZPlane(-10.0, boundary_type="reflective")
+    zmax = openmc.ZPlane(10.0, boundary_type="reflective")
+    cell1 = openmc.Cell(fill=m1, region=+xmin & -xmid & +ymin & -ymax & +zmin & -zmax)
+    cell2 = openmc.Cell(fill=m2, region=+xmid & -xmax & +ymin & -ymax & +zmin & -zmax)
+    model.geometry = openmc.Geometry([cell1, cell2])
+
+    src = openmc.IndependentSource()
+    src.space = openmc.stats.Point((-5.0, 0.0, 0.0))
+    src.particle = 'photon'
+    src.energy = openmc.stats.Discrete([1e6],[1.0])
+    
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 1
+    model.settings.particles = 100
+    model.settings.source = src
+
+    return model, xmid, cell1, cell2, m1, m2
+
+@pytest.fixture
+def two_sphere_model():
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+
+    mat = openmc.Material()
+    mat.add_nuclide('H2', 1.0)
+    mat.set_density('g/cm3', 1.0)
+
+    surf1 = openmc.Sphere(r=0.01)
+    surf2 = openmc.Sphere(r=1000, boundary_type='vacuum')
+    cell1 = openmc.Cell(region=-surf1)
+    cell2 = openmc.Cell(fill=mat, region=+surf1 & -surf2)
+
+    model.geometry = openmc.Geometry([cell1, cell2])
+
+    src = openmc.IndependentSource()
+    src.energy = openmc.stats.Discrete([1e6],[1.0])
+
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 1
+    model.settings.particles = 1000
+    model.settings.source = src
+
+    return model, cell1, cell2, mat
+
+def test_energy_balance_simple(two_cell_model, run_in_tmpdir):
+    model, xmid, cell1, cell2, *_  = two_cell_model
+
+    tally1 = openmc.Tally()
+    tally1.scores = ['heating']
+    tally1.filters = [openmc.CellFilter([cell1, cell2])]
+
+    tally2 = openmc.Tally()
+    tally2.scores = ['current']
+    tally2.filters = [openmc.EnergyFunctionFilter([0.0,20e6],[0.0,20e6]), openmc.SurfaceFilter([xmid])]
+
+    
+    model.tallies = [tally1, tally2]
+    model.run(apply_tally_results=True)
+
+    assert tally1.mean.sum() == pytest.approx(1e6)
+
+    assert tally2.mean[0] == pytest.approx(tally1.mean[1])
+
+def test_current_conservation(two_cell_model, run_in_tmpdir):
+    model, xmid, cell1, cell2, m1, m2 = two_cell_model
+
+    energyfunc_filter = openmc.EnergyFunctionFilter([0.0,20e6],[0.0,20e6])
+
+    tally1 = openmc.Tally()
+    tally1.scores = ['current']
+    tally1.filters = [energyfunc_filter, openmc.SurfaceFilter([xmid])]
+
+    tally2 = openmc.Tally()
+    tally2.scores = ['current']
+    tally2.filters = [energyfunc_filter,
+                      openmc.CellFromFilter([cell1]),
+                      openmc.CellFilter([cell2])]
+
+    tally3 = openmc.Tally()
+    tally3.scores = ['current']
+    tally3.filters = [energyfunc_filter,
+                      openmc.CellFromFilter([cell2]),
+                      openmc.CellFilter([cell1])]
+
+    tally4 = openmc.Tally()
+    tally4.scores = ['current']
+    tally4.filters = [energyfunc_filter,
+                      openmc.MaterialFromFilter([m1]),
+                      openmc.MaterialFilter([m2])]
+
+    tally5 = openmc.Tally()
+    tally5.scores = ['current']
+    tally5.filters = [energyfunc_filter,
+                      openmc.MaterialFromFilter([m2]),
+                      openmc.MaterialFilter([m1])]
+    
+    
+    model.tallies = [tally1, tally2, tally3, tally4, tally5]
+    model.run(apply_tally_results=True)
+
+    assert pytest.approx(tally1.mean.sum()) == tally2.mean[0] + tally3.mean[0]
+
+    assert tally2.mean[0] == pytest.approx(tally4.mean[0])
+
+    assert tally3.mean[0] == pytest.approx(tally5.mean[0])
+
+def test_cellfrom_heating(run_in_tmpdir, two_sphere_model):
+
+    model, cell1, cell2, mat = two_sphere_model
+
+    tally1 = openmc.Tally()
+    tally1.scores = ['heating']
+    tally1.filters = [openmc.CellFromFilter([cell1, cell2]), openmc.CellFilter([cell2])]
+
+    tally2 = openmc.Tally()
+    tally2.scores = ['heating']
+    tally2.filters = [openmc.MaterialFromFilter([mat]), openmc.MaterialFilter([mat])]
+
+    model.tallies = [tally1, tally2]
+
+    model.run(apply_tally_results=True)
+
+    assert np.all(tally1.mean > 0)
+
+    assert tally1.mean[1] == pytest.approx(tally2.mean[0]) 
+    
diff --git a/tests/unit_tests/test_r2s.py b/tests/unit_tests/test_r2s.py
index a94f85c8c08..f0b325573eb 100644
--- a/tests/unit_tests/test_r2s.py
+++ b/tests/unit_tests/test_r2s.py
@@ -53,7 +53,7 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path):
     r2s = R2SManager(model, mesh)
 
     # Use custom reduced chain file for Ni
-    chain = Chain.from_xml(Path(__file__).parents[1] / "chain_ni.xml")
+    chain_file = Path(__file__).parents[1] / "chain_ni.xml"
 
     # Run R2S calculation
     outdir = r2s.run(
@@ -61,7 +61,7 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path):
         source_rates=[1.0],
         photon_time_indices=[1],
         output_dir=tmp_path,
-        chain_file=chain,
+        chain_file=chain_file,
     )
 
     # Check directories and files exist
@@ -105,7 +105,7 @@ def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path):
     r2s = R2SManager(model, [c1, c2])
 
     # Use custom reduced chain file for Ni
-    chain = Chain.from_xml(Path(__file__).parents[1] / "chain_ni.xml")
+    chain_file = Path(__file__).parents[1] / "chain_ni.xml"
 
     # Run R2S calculation
     bounding_boxes = {c1.id: c1.bounding_box, c2.id: c2.bounding_box}
@@ -115,7 +115,7 @@ def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path):
         photon_time_indices=[1],
         output_dir=tmp_path,
         bounding_boxes=bounding_boxes,
-        chain_file=chain
+        chain_file=chain_file
     )
 
     # Check directories and files exist