diff --git a/docs/source/methods/neutron_physics.rst b/docs/source/methods/neutron_physics.rst
index 70ace4a3532..9515055692f 100644
--- a/docs/source/methods/neutron_physics.rst
+++ b/docs/source/methods/neutron_physics.rst
@@ -290,7 +290,10 @@ create and store fission sites for the following generation. First, the average
 number of prompt and delayed neutrons must be determined to decide whether the
 secondary neutrons will be prompt or delayed. This is important because delayed
 neutrons have a markedly different spectrum from prompt neutrons, one that has a
-lower average energy of emission. The total number of neutrons emitted
+lower average energy of emission. Furthermore, in simulations where tracking
+time of neutrons is important, we need to consider the emission time delay of
+the secondary neutrons, which is dependent on the decay constant of the
+delayed neutron precursor. The total number of neutrons emitted
 :math:`\nu_t` is given as a function of incident energy in the ENDF format. Two
 representations exist for :math:`\nu_t`. The first is a polynomial of order
 :math:`N` with coefficients :math:`c_0,c_1,\dots,c_N`. If :math:`\nu_t` has this
@@ -306,8 +309,8 @@ interpolation law. The number of prompt neutrons released per fission event
 :math:`\nu_p` is also given as a function of incident energy and can be
 specified in a polynomial or tabular format. The number of delayed neutrons
 released per fission event :math:`\nu_d` can only be specified in a tabular
-format. In practice, we only need to determine :math:`nu_t` and
-:math:`nu_d`. Once these have been determined, we can calculated the delayed
+format. In practice, we only need to determine :math:`\nu_t` and
+:math:`\nu_d`. Once these have been determined, we can calculate the delayed
 neutron fraction
 
 .. math::
@@ -335,8 +338,14 @@ neutrons. Otherwise, we produce :math:`\lfloor \nu \rfloor + 1` neutrons. Then,
 for each fission site produced, we sample the outgoing angle and energy
 according to the algorithms given in :ref:`sample-angle` and
 :ref:`sample-energy` respectively. If the neutron is to be born delayed, then
-there is an extra step of sampling a delayed neutron precursor group since they
-each have an associated secondary energy distribution.
+there is an extra step of sampling a delayed neutron precursor group to get the
+associated secondary energy distribution and the decay constant
+:math:`\lambda`, which is needed to sample the emission delay time :math:`t_d`:
+
+.. math::
+    :label: sample-delay-time
+
+    t_d = -\frac{\ln \xi}{\lambda}.
 
 The sampled outgoing angle and energy of fission neutrons along with the
 position of the collision site are stored in an array called the fission
diff --git a/examples/pincell_pulsed/build_xml.py b/examples/pincell_pulsed/build_xml.py
new file mode 100644
index 00000000000..7cd1617dfeb
--- /dev/null
+++ b/examples/pincell_pulsed/build_xml.py
@@ -0,0 +1,86 @@
+import numpy as np
+import openmc
+
+###############################################################################
+# Create materials for the problem
+
+uo2 = openmc.Material(name="UO2 fuel at 2.4% wt enrichment")
+uo2.set_density("g/cm3", 10.29769)
+uo2.add_element("U", 1.0, enrichment=2.4)
+uo2.add_element("O", 2.0)
+
+helium = openmc.Material(name="Helium for gap")
+helium.set_density("g/cm3", 0.001598)
+helium.add_element("He", 2.4044e-4)
+
+zircaloy = openmc.Material(name="Zircaloy 4")
+zircaloy.set_density("g/cm3", 6.55)
+zircaloy.add_element("Sn", 0.014, "wo")
+zircaloy.add_element("Fe", 0.00165, "wo")
+zircaloy.add_element("Cr", 0.001, "wo")
+zircaloy.add_element("Zr", 0.98335, "wo")
+
+borated_water = openmc.Material(name="Borated water")
+borated_water.set_density("g/cm3", 0.740582)
+borated_water.add_element("B", 2.0e-4)  # 3x the original pincell
+borated_water.add_element("H", 5.0e-2)
+borated_water.add_element("O", 2.4e-2)
+borated_water.add_s_alpha_beta("c_H_in_H2O")
+
+###############################################################################
+# Define problem geometry
+
+# Create cylindrical surfaces
+fuel_or = openmc.ZCylinder(r=0.39218, name="Fuel OR")
+clad_ir = openmc.ZCylinder(r=0.40005, name="Clad IR")
+clad_or = openmc.ZCylinder(r=0.45720, name="Clad OR")
+
+# Create a region represented as the inside of a rectangular prism
+pitch = 1.25984
+box = openmc.model.RectangularPrism(pitch, pitch, boundary_type="reflective")
+
+# Create cells, mapping materials to regions
+fuel = openmc.Cell(fill=uo2, region=-fuel_or)
+gap = openmc.Cell(fill=helium, region=+fuel_or & -clad_ir)
+clad = openmc.Cell(fill=zircaloy, region=+clad_ir & -clad_or)
+water = openmc.Cell(fill=borated_water, region=+clad_or & -box)
+
+# Create a model and assign geometry
+model = openmc.Model()
+model.geometry = openmc.Geometry([fuel, gap, clad, water])
+
+###############################################################################
+# Define problem settings
+
+# Set the mode
+model.settings.run_mode = "fixed source"
+
+# Indicate how many batches and particles to run
+model.settings.batches = 10
+model.settings.particles = 10000
+
+# Set time cutoff
+# (because we only care about solutions in t < 100 seconds, see tally below)
+model.settings.cutoff = {"time_neutron": 100}
+
+# Create the neutron pulse source
+# (by default, isotropic direction, at t = 0 s)
+space = openmc.stats.Point()  # At the origin (0, 0, 0)
+energy = openmc.stats.Discrete([1.41e7], [1.0])  # At 14.1 MeV
+model.settings.source = openmc.IndependentSource(space=space, energy=energy)
+
+###############################################################################
+# Define tallies
+
+# Create time filter
+t_grid = np.insert(np.logspace(-6, 2, 100), 0, 0.0)
+time_filter = openmc.TimeFilter(t_grid)
+
+# Tally for total neutron density in time
+density_tally = openmc.Tally(name="Density")
+density_tally.filters = [time_filter]
+density_tally.scores = ["inverse-velocity"]
+
+# Add tallies to model
+model.tallies = openmc.Tallies([density_tally])
+model.export_to_model_xml()
diff --git a/examples/pincell_pulsed/plot_result.py b/examples/pincell_pulsed/plot_result.py
new file mode 100644
index 00000000000..29dae470c5f
--- /dev/null
+++ b/examples/pincell_pulsed/plot_result.py
@@ -0,0 +1,24 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import openmc
+
+
+# Get results from statepoint
+with openmc.StatePoint("statepoint.10.h5") as sp:
+    tally = sp.get_tally(name="Density")
+
+# Get the time grid
+t = tally.filters[0].values
+dt = np.diff(t)
+t_mid = 0.5 * (t[:-1] + t[1:])
+
+# Bin-averaged result
+density_mean = tally.mean.ravel() / dt
+
+# Plot flux spectrum
+fig, ax = plt.subplots()
+ax.loglog(t_mid, density_mean, "ok", fillstyle="none")
+ax.set_xlabel("Time [s]")
+ax.set_ylabel("Total density")
+ax.grid()
+plt.show()
diff --git a/include/openmc/settings.h b/include/openmc/settings.h
index 1e44c08801d..f2e306e059e 100644
--- a/include/openmc/settings.h
+++ b/include/openmc/settings.h
@@ -132,6 +132,7 @@ extern std::unordered_set<int>
 
 extern int
   max_history_splits; //!< maximum number of particle splits for weight windows
+extern int max_secondaries; //!< maximum number of secondaries in the bank
 extern int64_t ssw_max_particles; //!< maximum number of particles to be
                                   //!< banked on surfaces per process
 extern int64_t ssw_max_files;     //!< maximum number of surface source files
diff --git a/openmc/settings.py b/openmc/settings.py
index 0a78fb564f8..6c7186ac273 100644
--- a/openmc/settings.py
+++ b/openmc/settings.py
@@ -120,6 +120,10 @@ class Settings:
         Maximum number of times a particle can split during a history
 
         .. versionadded:: 0.13
+    max_secondaries : int
+        Maximum secondary bank size
+
+        .. versionadded:: 0.14.1
     max_tracks : int
         Maximum number of tracks written to a track file (per MPI process).
 
@@ -1056,6 +1060,16 @@ def max_history_splits(self, value: int):
         cv.check_greater_than('max particle splits', value, 0)
         self._max_history_splits = value
 
+    @property
+    def max_secondaries(self) -> int:
+        return self._max_secondaries
+
+    @max_secondaries.setter
+    def max_secondaries(self, value: int):
+        cv.check_type('maximum secondary bank size', value, Integral)
+        cv.check_greater_than('max secondary bank size', value, 0)
+        self._max_secondaries = value
+
     @property
     def max_tracks(self) -> int:
         return self._max_tracks
@@ -1528,6 +1542,11 @@ def _create_max_history_splits_subelement(self, root):
             elem = ET.SubElement(root, "max_history_splits")
             elem.text = str(self._max_history_splits)
 
+    def _create_max_secondaries_subelement(self, root):
+        if self._max_secondaries is not None:
+            elem = ET.SubElement(root, "max_secondaries")
+            elem.text = str(self._max_secondaries)
+
     def _create_max_tracks_subelement(self, root):
         if self._max_tracks is not None:
             elem = ET.SubElement(root, "max_tracks")
@@ -1893,6 +1912,11 @@ def _max_history_splits_from_xml_element(self, root):
         if text is not None:
             self.max_history_splits = int(text)
 
+    def _max_secondaries_from_xml_element(self, root):
+        text = get_text(root, 'max_secondaries')
+        if text is not None:
+            self.max_secondaries = int(text)
+
     def _max_tracks_from_xml_element(self, root):
         text = get_text(root, 'max_tracks')
         if text is not None:
diff --git a/src/finalize.cpp b/src/finalize.cpp
index 08c2fced308..0564c0454cd 100644
--- a/src/finalize.cpp
+++ b/src/finalize.cpp
@@ -91,6 +91,7 @@ int openmc_finalize()
   settings::max_lost_particles = 10;
   settings::max_order = 0;
   settings::max_particles_in_flight = 100000;
+  settings::max_secondaries = 10000;
   settings::max_particle_events = 1'000'000;
   settings::max_history_splits = 10'000'000;
   settings::max_tracks = 1000;
diff --git a/src/particle.cpp b/src/particle.cpp
index 64c50c9438f..74638be806d 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -44,6 +44,15 @@ namespace openmc {
 
 double Particle::speed() const
 {
+  // Multigroup speed?
+  if (!settings::run_CE) {
+    auto& macro_xs = data::mg.macro_xs_[this->material()];
+    int macro_t = this->mg_xs_cache().t;
+    int macro_a = macro_xs.get_angle_index(this->u());
+    return 1.0 / macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, this->g(), nullptr,
+                   nullptr, nullptr, macro_t, macro_a);
+  }
+
   // Determine mass in eV/c^2
   double mass;
   switch (this->type()) {
@@ -244,6 +253,9 @@ void Particle::event_advance()
     double push_back_distance = speed() * dt;
     this->move_distance(-push_back_distance);
     hit_time_boundary = true;
+
+    // Reduce the distance traveled for tallying
+    distance -= push_back_distance;
   }
 
   // Score track-length tallies
diff --git a/src/physics.cpp b/src/physics.cpp
index 69e74f2eafd..683293bfb49 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -113,7 +113,7 @@ void sample_neutron_reaction(Particle& p)
 
       // Make sure particle population doesn't grow out of control for
       // subcritical multiplication problems.
-      if (p.secondary_bank().size() >= 10000) {
+      if (p.secondary_bank().size() >= settings::max_secondaries) {
         fatal_error(
           "The secondary particle bank appears to be growing without "
           "bound. You are likely running a subcritical multiplication problem "
@@ -202,13 +202,23 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
     site.particle = ParticleType::neutron;
     site.time = p.time();
     site.wgt = 1. / weight;
-    site.parent_id = p.id();
-    site.progeny_id = p.n_progeny()++;
     site.surf_id = 0;
 
     // Sample delayed group and angle/energy for fission reaction
     sample_fission_neutron(i_nuclide, rx, &site, p);
 
+    // If neutron is delayed, sample time of emission
+    if (site.delayed_group > 0) {
+      double decay_rate = rx.products_[site.delayed_group].decay_rate_;
+      site.time -= std::log(prn(p.current_seed())) / decay_rate;
+
+      // Reject site if it exceeds time cutoff
+      double t_cutoff = settings::time_cutoff[static_cast<int>(site.particle)];
+      if (site.time > t_cutoff) {
+        continue;
+      }
+    }
+
     // Store fission site in bank
     if (use_fission_bank) {
       int64_t idx = simulation::fission_bank.thread_safe_append(site);
@@ -218,11 +228,6 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
           "in this generation will not be banked. Results may be "
           "non-deterministic.");
 
-        // Decrement number of particle progeny as storage was unsuccessful.
-        // This step is needed so that the sum of all progeny is equal to the
-        // size of the shared fission bank.
-        p.n_progeny()--;
-
         // Break out of loop as no more sites can be added to fission bank
         break;
       }
@@ -230,6 +235,10 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
       p.secondary_bank().push_back(site);
     }
 
+    // Set parent and progeny ID
+    site.parent_id = p.id();
+    site.progeny_id = p.n_progeny()++;
+
     // Set the delayed group on the particle as well
     p.delayed_group() = site.delayed_group;
 
diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp
index 361cf5affcf..cbe954cc9a4 100644
--- a/src/physics_mg.cpp
+++ b/src/physics_mg.cpp
@@ -127,9 +127,8 @@ void create_fission_sites(Particle& p)
     SourceSite site;
     site.r = p.r();
     site.particle = ParticleType::neutron;
+    site.time = p.time();
     site.wgt = 1. / weight;
-    site.parent_id = p.id();
-    site.progeny_id = p.n_progeny()++;
 
     // Sample the cosine of the angle, assuming fission neutrons are emitted
     // isotropically
@@ -154,6 +153,20 @@ void create_fission_sites(Particle& p)
     // of the code, 0 is prompt.
     site.delayed_group = dg + 1;
 
+    // If delayed product production, sample time of emission
+    if (dg != -1) {
+      auto& macro_xs = data::mg.macro_xs_[p.material()];
+      double decay_rate =
+        macro_xs.get_xs(MgxsType::DECAY_RATE, 0, nullptr, nullptr, &dg, 0, 0);
+      site.time -= std::log(prn(p.current_seed())) / decay_rate;
+
+      // Reject site if it exceeds time cutoff
+      double t_cutoff = settings::time_cutoff[static_cast<int>(site.particle)];
+      if (site.time > t_cutoff) {
+        continue;
+      }
+    }
+
     // Store fission site in bank
     if (use_fission_bank) {
       int64_t idx = simulation::fission_bank.thread_safe_append(site);
@@ -163,11 +176,6 @@ void create_fission_sites(Particle& p)
           "in this generation will not be banked. Results may be "
           "non-deterministic.");
 
-        // Decrement number of particle progeny as storage was unsuccessful.
-        // This step is needed so that the sum of all progeny is equal to the
-        // size of the shared fission bank.
-        p.n_progeny()--;
-
         // Break out of loop as no more sites can be added to fission bank
         break;
       }
@@ -175,6 +183,10 @@ void create_fission_sites(Particle& p)
       p.secondary_bank().push_back(site);
     }
 
+    // Set parent and progeny ID
+    site.parent_id = p.id();
+    site.progeny_id = p.n_progeny()++;
+
     // Set the delayed group on the particle as well
     p.delayed_group() = dg + 1;
 
diff --git a/src/settings.cpp b/src/settings.cpp
index 3a905bdf961..f60a8f3bb9f 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -107,6 +107,7 @@ int max_order {0};
 int n_log_bins {8000};
 int n_batches;
 int n_max_batches;
+int max_secondaries {10000};
 int max_history_splits {10'000'000};
 int max_tracks {1000};
 ResScatMethod res_scat_method {ResScatMethod::rvs};
@@ -1056,6 +1057,11 @@ void read_settings_xml(pugi::xml_node root)
     weight_windows_on = get_node_value_bool(root, "weight_windows_on");
   }
 
+  if (check_for_node(root, "max_secondaries")) {
+    settings::max_secondaries =
+      std::stoi(get_node_value(root, "max_secondaries"));
+  }
+
   if (check_for_node(root, "max_history_splits")) {
     settings::max_history_splits =
       std::stoi(get_node_value(root, "max_history_splits"));
diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index e798161ec2f..fc3236b8aee 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -1558,8 +1558,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
       tally.estimator_ == TallyEstimator::COLLISION) {
     if (settings::survival_biasing) {
       // Determine weight that was absorbed
-      wgt_absorb = p.wgt_last() * p.neutron_xs(p.event_nuclide()).absorption /
-                   p.neutron_xs(p.event_nuclide()).total;
+      wgt_absorb = p.wgt_last() * p.macro_xs().absorption / p.macro_xs().total;
 
       // Then we either are alive and had a scatter (and so g changed),
       // or are dead and g did not change
diff --git a/tests/unit_tests/test_settings.py b/tests/unit_tests/test_settings.py
index 02a47625162..1179c2bbaf5 100644
--- a/tests/unit_tests/test_settings.py
+++ b/tests/unit_tests/test_settings.py
@@ -65,8 +65,8 @@ def test_export_to_xml(run_in_tmpdir):
             space=openmc.stats.Box((-1., -1., -1.), (1., 1., 1.))
         )
     }
-
     s.max_particle_events = 100
+    s.max_secondaries = 1_000_000
 
     # Make sure exporting XML works
     s.export_to_xml()
@@ -142,3 +142,4 @@ def test_export_to_xml(run_in_tmpdir):
     assert s.random_ray['distance_active'] == 100.0
     assert s.random_ray['ray_source'].space.lower_left == [-1., -1., -1.]
     assert s.random_ray['ray_source'].space.upper_right == [1., 1., 1.]
+    assert s.max_secondaries == 1_000_000