diff --git a/doc/source/userguide/discrete_ordinates_problems.rst b/doc/source/userguide/discrete_ordinates_problems.rst
index 7f9b1bc99..87d6983b0 100644
--- a/doc/source/userguide/discrete_ordinates_problems.rst
+++ b/doc/source/userguide/discrete_ordinates_problems.rst
@@ -208,10 +208,18 @@ The most important options for everyday use are:
 * ``ags_convergence_check``
 * ``field_function_prefix_option``
 * ``field_function_prefix``
+* ``use_precursors``
 
 There are also restart, precursor, and message-size options for more
 specialized workflows.
 
+``use_precursors`` controls precursor storage and transient delayed-neutron
+evolution. Steady-state and k-eigenvalue solves include delayed neutron
+production whenever the active cross-section data contains it, independent of
+this option. In transient mode, ``use_precursors=False`` selects a prompt-only
+transient, while ``use_precursors=True`` evolves precursor concentrations and
+includes delayed production through those precursor equations.
+
 Restart-related options include:
 
 * ``restart_writes_enabled``: enable restart dump writes.
diff --git a/doc/source/userguide/materials_xs.rst b/doc/source/userguide/materials_xs.rst
index 394dbb5ab..edf0c9038 100644
--- a/doc/source/userguide/materials_xs.rst
+++ b/doc/source/userguide/materials_xs.rst
@@ -267,11 +267,88 @@ fission production. This is useful for advanced data preparation, but most user
 inputs are clearer when written in terms of ``SIGMA_F``, ``NU`` or
 ``NU_PROMPT``/``NU_DELAYED``, and the corresponding spectra.
 
+Delayed-Neutron Data and Solver Behavior
+----------------------------------------
+
+Delayed-neutron data has two related but distinct roles in OpenSn:
+
+* it contributes delayed fission production to steady-state and k-eigenvalue
+  fission sources, and
+* it can define precursor concentrations for time-dependent delayed-neutron
+  evolution.
+
+These roles are intentionally not controlled by the same switch. If a
+steady-state or k-eigenvalue problem uses cross sections with delayed-neutron
+data, OpenSn includes delayed neutron production in the steady fission source
+and in fission-production normalization. The steady eigenvalue therefore
+corresponds to total fission production, prompt plus delayed. This is true even
+when the problem option ``use_precursors`` is ``False``.
+
+For transient problems, ``use_precursors`` controls whether delayed-neutron
+precursor state is stored and evolved:
+
+* ``use_precursors=True`` uses prompt production and precursor-mediated delayed
+  production.
+* ``use_precursors=False`` is a prompt-only transient, even if the cross-section
+  file contains delayed-neutron data.
+
+This distinction is useful when a steady k-eigenvalue solve is used only to
+prepare an initial condition. The steady solve still computes the physical
+total-production eigenvalue, while a later transient can deliberately choose
+whether to evolve delayed precursors.
+
+.. important::
+
+   If delayed-neutron production is active and any fissionable material in the
+   active cross-section map contains precursor data, all fissionable materials
+   in that map must contain compatible precursor data. Non-fissionable materials
+   may have zero precursors. This avoids solving a mixed prompt-only and
+   prompt-plus-delayed fission model by accident.
+
+Native OpenSn delayed data is interpreted as follows:
+
+* ``NU_PROMPT`` and ``CHI_PROMPT`` define the prompt fission production matrix.
+* ``NU_DELAYED`` defines the delayed fission production by source group.
+* ``PRECURSOR_FRACTIONAL_YIELDS`` partitions delayed production among precursor
+  families. The yields are normalized when read.
+* ``CHI_DELAYED`` defines each precursor family's delayed emission spectrum by
+  destination group.
+* ``PRECURSOR_DECAY_CONSTANTS`` are used for precursor concentrations and
+  transient evolution; they do not change the steady-state k-eigenvalue.
+
+For steady-state source assembly and k-eigenvalue normalization, OpenSn adds the
+delayed term
+
+.. math::
+
+   \sum_j \chi_{d,j,g}\,\gamma_j\,\nu_{d,g'}\Sigma_{f,g'}\phi_{g'}
+
+to the prompt production matrix contribution. Here ``j`` is the precursor family
+index, ``g`` is the destination group, and ``g'`` is the source group.
+
+When the native ``PRODUCTION_MATRIX`` form is used, delayed-neutron precursor
+data is currently not allowed in the same material. This avoids ambiguity over
+whether the matrix already contains delayed production.
+
+Practical guidance:
+
+* For production steady-state or k-eigenvalue calculations, leave delayed data in
+  the XS when it exists. Do not disable ``use_precursors`` to try to remove
+  delayed neutron production from a steady eigenvalue solve.
+* For production transients with delayed-neutron physics, use
+  ``options={"use_precursors": True, "save_angular_flux": True}``.
+* Use ``use_precursors=False`` with delayed XS data only when a prompt-only
+  transient is intentional.
+* If a k-eigenvalue changes when delayed data is removed from the XS, that is
+  expected. If it changes only because ``use_precursors`` was toggled in a
+  steady-state calculation, that is not expected.
+
 Loading OpenMC MGXS Files
 =========================
 
 Use :py:meth:`pyopensn.xs.MultiGroupXS.LoadFromOpenMC` to load cross sections
-from an OpenMC multi-group HDF5 library.
+from an OpenMC multi-group HDF5 library. Delayed-neutron data is imported when
+it is present in the OpenMC file.
 
 .. code-block:: python
 
@@ -309,6 +386,13 @@ Example:
    absorption dataset directly. This makes the imported data less sensitive to
    small inconsistencies or statistical noise.
 
+   When delayed-neutron data is present, OpenSn also imports precursor-group
+   information, prompt and delayed fission production, precursor decay
+   constants, and delayed emission spectra. If ``chi-prompt`` is not present,
+   OpenSn falls back to the steady fission spectrum ``chi`` as the prompt
+   spectrum. Once imported, the same steady-state and transient delayed-neutron
+   behavior described above applies.
+
 What OpenMC Import Expects
 --------------------------
 
@@ -326,6 +410,14 @@ The imported object may include:
 * fission cross section
 * fission production data
 * fission spectrum
+* delayed-neutron data, including:
+
+  - ``prompt-nu-fission``
+  - ``delayed-nu-fission``
+  - ``decay-rate``
+  - ``chi-delayed``
+  - optional ``chi-prompt``
+
 * any requested named custom one-dimensional XS in ``extra_xs_names``
 
 If the requested dataset or temperature is not present, the load fails with an
diff --git a/doc/source/userguide/solvers.rst b/doc/source/userguide/solvers.rst
index 9bddc660d..05920bf29 100644
--- a/doc/source/userguide/solvers.rst
+++ b/doc/source/userguide/solvers.rst
@@ -178,10 +178,18 @@ Problem options available in Python include:
 * ``field_function_prefix_option``
 * ``field_function_prefix``
 
-``use_precursors`` controls whether delayed-neutron precursor treatment is kept
-active for the problem. The default is ``True``. This should usually stay
-enabled for transient and k-eigen workflows unless you explicitly want a
-prompt-only model.
+``use_precursors`` controls delayed-neutron precursor storage and transient
+precursor evolution. It does not remove delayed neutron production from
+steady-state or k-eigenvalue fission sources. If the active cross sections
+contain delayed-neutron data, steady-state and k-eigenvalue solves use total
+fission production, prompt plus delayed, regardless of ``use_precursors``.
+
+For transient solves, ``use_precursors=True`` evolves precursor concentrations
+and includes delayed production through the transient precursor equations.
+``use_precursors=False`` is a prompt-only transient, even if the cross sections
+contain delayed-neutron data. Keep the default ``True`` for production
+delayed-neutron transients; set it to ``False`` only when a prompt-only transient
+is intentional.
 
 ``read_restart_path`` reads a full restart for continuing a compatible solve.
 ``read_initial_condition_path`` reads restart data as an initial condition; this
@@ -205,9 +213,11 @@ new families start at zero and removed families are discarded. If a cell passes
 through a material with zero precursors, its precursor history is dropped and
 any later reintroduced precursor families restart from zero.
 
-If any fissionable material in the active map contains precursor data and
-``use_precursors=True``, then all fissionable materials in that map must
-contain precursor data. Non-fissionable materials may have zero precursors.
+If delayed-neutron production is active and any fissionable material in the
+active map contains precursor data, then all fissionable materials in that map
+must contain precursor data. Non-fissionable materials may have zero precursors.
+In steady-state mode, delayed-neutron production is active whenever delayed data
+exists. In transient mode, it is active only when ``use_precursors=True``.
 
 .. note::
 
@@ -516,6 +526,15 @@ Meaning of ``initial_state``:
    ``initial_state="existing"`` is the usual choice when a transient begins
    from a previously computed steady-state or eigenvalue solution.
 
+.. note::
+
+   A steady-state or k-eigenvalue initial condition computed with delayed
+   cross-section data uses total fission production. A following transient uses
+   delayed precursors only if the problem option ``use_precursors`` is ``True``.
+   If ``use_precursors`` is ``False``, the transient is prompt-only even though
+   the steady initial condition may have included delayed production in its
+   eigenvalue or source balance.
+
 .. note::
 
    ``initial_state="zero"`` is mostly useful for controlled studies where the
@@ -730,6 +749,12 @@ At a high level, the solver:
 This solver uses the groupset inner solver and AGS settings configured on the
 problem.
 
+If the cross sections contain delayed-neutron data, the fission source and
+``k_eff`` update use total fission production, prompt plus delayed. This
+steady-state behavior is independent of ``use_precursors``. The
+``use_precursors`` option only controls whether precursor concentrations are
+stored and computed after the solve.
+
 .. note::
 
    Power iteration is usually the first k-eigen method to try because it uses
diff --git a/framework/materials/multi_group_xs/openmc_import.cc b/framework/materials/multi_group_xs/openmc_import.cc
index afedd7dda..5862221c4 100644
--- a/framework/materials/multi_group_xs/openmc_import.cc
+++ b/framework/materials/multi_group_xs/openmc_import.cc
@@ -8,6 +8,7 @@
 #include "framework/utils/utils.h"
 #include <numeric>
 #include <algorithm>
+#include <cstdint>
 
 namespace opensn
 {
@@ -71,7 +72,11 @@ MultiGroupXS::LoadFromOpenMC(const std::string& file_name,
     throw std::runtime_error(file_name + " has an unsupported scatter shape");
 
   // Number of precursors
-  mgxs.num_precursors_ = 0;
+  unsigned int delayed_groups = 0;
+  if (not H5ReadAttribute<unsigned int>(file.Id(), "delayed_groups", delayed_groups))
+    delayed_groups = 0;
+  OpenSnLogicalErrorIf(delayed_groups < 0, "OpenMC number of delayed_groups must be non-negative.");
+  mgxs.num_precursors_ = delayed_groups;
 
   // Inverse velocity
   H5ReadDataset1D<double>(file.Id(), path + "inverse-velocity", mgxs.inv_velocity_);
@@ -162,6 +167,96 @@ MultiGroupXS::LoadFromOpenMC(const std::string& file_name,
       if (mgxs.sigma_f_[g] > 0.0)
         nu[g] /= mgxs.sigma_f_[g];
 
+    if (mgxs.num_precursors_ > 0)
+    {
+      std::vector<double> nu_prompt_sigma_f;
+      if (H5Has(file.Id(), path + "prompt-nu-fission"))
+        H5ReadDataset1D<double>(file.Id(), path + "prompt-nu-fission", nu_prompt_sigma_f);
+
+      std::vector<std::vector<double>> delayed_nu_sigma_f_by_prec;
+      if (not H5ReadDataset2D<double>(
+            file.Id(), path + "delayed-nu-fission", delayed_nu_sigma_f_by_prec))
+      {
+        std::vector<double> delayed_nu_sigma_f_1d;
+        OpenSnLogicalErrorIf(not H5ReadDataset1D<double>(
+                               file.Id(), path + "delayed-nu-fission", delayed_nu_sigma_f_1d),
+                             "Failed to read dataset \"" + path + "delayed-nu-fission\".");
+        OpenSnLogicalErrorIf(
+          mgxs.num_precursors_ != 1 or delayed_nu_sigma_f_1d.size() != mgxs.num_groups_,
+          "Dataset \"" + path + "delayed-nu-fission\" has incorrect dimensions.");
+        delayed_nu_sigma_f_by_prec.push_back(std::move(delayed_nu_sigma_f_1d));
+      }
+      std::vector<double> nu_delayed_sigma_f(mgxs.num_groups_, 0.0);
+      for (const auto& delayed_family : delayed_nu_sigma_f_by_prec)
+        for (unsigned int g = 0; g < mgxs.num_groups_; ++g)
+          nu_delayed_sigma_f[g] += delayed_family[g];
+
+      if (nu_prompt_sigma_f.empty())
+      {
+        nu_prompt_sigma_f = mgxs.nu_sigma_f_;
+        for (unsigned int g = 0; g < mgxs.num_groups_; ++g)
+          nu_prompt_sigma_f[g] -= nu_delayed_sigma_f[g];
+      }
+
+      OpenSnLogicalErrorIf(not IsNonNegative(nu_prompt_sigma_f),
+                           "OpenMC delayed fission data imply a negative prompt nu-sigma-f.");
+      OpenSnLogicalErrorIf(not IsNonNegative(nu_delayed_sigma_f),
+                           "OpenMC delayed nu-sigma-f must be non-negative.");
+
+      std::vector<double> precursor_decay_constants;
+      H5ReadDataset1D<double>(file.Id(), path + "decay-rate", precursor_decay_constants);
+      OpenSnLogicalErrorIf(precursor_decay_constants.size() != mgxs.num_precursors_,
+                           "OpenMC decay-rate data has incorrect size.");
+
+      std::vector<std::vector<double>> chi_delayed_by_prec;
+      if (not H5ReadDataset2D<double>(file.Id(), path + "chi-delayed", chi_delayed_by_prec))
+      {
+        std::vector<double> chi_delayed_1d;
+        OpenSnLogicalErrorIf(
+          not H5ReadDataset1D<double>(file.Id(), path + "chi-delayed", chi_delayed_1d),
+          "Failed to read dataset \"" + path + "chi-delayed\".");
+        OpenSnLogicalErrorIf(mgxs.num_precursors_ != 1 or chi_delayed_1d.size() != mgxs.num_groups_,
+                             "Dataset \"" + path + "chi-delayed\" has incorrect dimensions.");
+        chi_delayed_by_prec.push_back(std::move(chi_delayed_1d));
+      }
+
+      std::vector<double> chi_prompt = mgxs.chi_;
+      if (H5Has(file.Id(), path + "chi-prompt"))
+      {
+        H5ReadDataset1D<double>(file.Id(), path + "chi-prompt", chi_prompt);
+        OpenSnLogicalErrorIf(chi_prompt.size() != mgxs.num_groups_,
+                             "OpenMC chi-prompt data has incorrect size.");
+      }
+      else
+      {
+        log.Log0Warning() << "OpenMC delayed-neutron data found in \"" << file_name
+                          << "\" without chi-prompt. Using chi as the prompt spectrum.";
+      }
+
+      mgxs.nu_prompt_sigma_f_ = std::move(nu_prompt_sigma_f);
+      mgxs.nu_delayed_sigma_f_ = std::move(nu_delayed_sigma_f);
+      mgxs.precursors_.resize(mgxs.num_precursors_);
+
+      double delayed_total = 0.0;
+      for (const auto& delayed_family : delayed_nu_sigma_f_by_prec)
+        delayed_total += std::accumulate(delayed_family.begin(), delayed_family.end(), 0.0);
+
+      for (unsigned int j = 0; j < mgxs.num_precursors_; ++j)
+      {
+        auto& precursor = mgxs.precursors_[j];
+        precursor.decay_constant = precursor_decay_constants[j];
+        precursor.emission_spectrum = chi_delayed_by_prec[j];
+        const double family_total = std::accumulate(
+          delayed_nu_sigma_f_by_prec[j].begin(), delayed_nu_sigma_f_by_prec[j].end(), 0.0);
+        precursor.fractional_yield = delayed_total > 0.0 ? family_total / delayed_total : 0.0;
+      }
+
+      mgxs.production_matrix_.assign(mgxs.num_groups_, std::vector<double>(mgxs.num_groups_, 0.0));
+      for (unsigned int gp = 0; gp < mgxs.num_groups_; ++gp)
+        for (unsigned int g = 0; g < mgxs.num_groups_; ++g)
+          mgxs.production_matrix_[g][gp] = chi_prompt[g] * mgxs.nu_prompt_sigma_f_[gp];
+    }
+
     // Production matrix (computed)
     // TODO: This path uses chi * nu_sigma_f (total production). If delayed-neutron data is
     // introduced here in the future, ensure LBS source/fission-production routines remain
@@ -189,6 +284,7 @@ MultiGroupXS::LoadFromOpenMC(const std::string& file_name,
   else
   {
     // Clear fission data if not fissionable
+    mgxs.num_precursors_ = 0;
     mgxs.sigma_f_.clear();
     mgxs.nu_sigma_f_.clear();
     mgxs.nu_prompt_sigma_f_.clear();
diff --git a/framework/utils/hdf_utils.h b/framework/utils/hdf_utils.h
index e95639498..1cda8ff17 100644
--- a/framework/utils/hdf_utils.h
+++ b/framework/utils/hdf_utils.h
@@ -312,6 +312,52 @@ H5ReadDataset1D(hid_t id, const std::string& name, std::vector<T>& data)
   return success;
 }
 
+template <typename T>
+bool
+H5ReadDataset2D(hid_t id, const std::string& name, std::vector<std::vector<T>>& data)
+{
+  static_assert(
+    !std::is_same_v<T, bool>,
+    "H5ReadDataset2D does not support std::vector<std::vector<bool>>. Use unsigned char.");
+
+  bool success = false;
+  data.clear();
+
+  auto dataset = H5Dopen2(id, name.c_str(), H5P_DEFAULT);
+  if (dataset == H5I_INVALID_HID)
+    return false;
+
+  auto dataspace = H5Dget_space(dataset);
+  if (dataspace != H5I_INVALID_HID)
+  {
+    const int rank = H5Sget_simple_extent_ndims(dataspace);
+    if (rank == 2)
+    {
+      hsize_t dims[2] = {0, 0};
+      if (H5Sget_simple_extent_dims(dataspace, dims, nullptr) >= 0)
+      {
+        auto n0 = static_cast<size_t>(dims[0]);
+        auto n1 = static_cast<size_t>(dims[1]);
+        std::vector<T> flat_data(n0 * n1);
+        if (flat_data.empty() or
+            H5Dread(dataset, get_datatype<T>(), H5S_ALL, H5S_ALL, H5P_DEFAULT, flat_data.data()) >=
+              0)
+        {
+          data.assign(n0, std::vector<T>(n1));
+          for (size_t i = 0; i < n0; ++i)
+            for (size_t j = 0; j < n1; ++j)
+              data[i][j] = flat_data[i * n1 + j];
+          success = true;
+        }
+      }
+    }
+    H5Sclose(dataspace);
+  }
+
+  H5Dclose(dataset);
+  return success;
+}
+
 template <typename T>
 bool
 H5ReadAttribute(hid_t id, const std::string& name, T& value)
diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/acceleration/cmfd_acceleration.cc b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/acceleration/cmfd_acceleration.cc
index 100444c9e..82a29c2bb 100644
--- a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/acceleration/cmfd_acceleration.cc
+++ b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/acceleration/cmfd_acceleration.cc
@@ -943,6 +943,7 @@ CMFDAcceleration::CondensedProductionXS(const MultiGroupXS& xs,
                                         const unsigned int col_coarse_group) const
 {
   const auto& F = xs.GetProductionMatrix();
+  const bool include_delayed_production = UseDelayedNeutronProduction(do_problem_);
   const bool has_flux = not fine_coarse_phi.empty();
   const auto fine_offset = coarse_cell.local_id * num_groups_;
   double numerator = 0.0;
@@ -953,7 +954,11 @@ CMFDAcceleration::CondensedProductionXS(const MultiGroupXS& xs,
     const double phi = has_flux ? fine_coarse_phi[fine_offset + gp] : 1.0;
     denominator += phi;
     for (unsigned int g = FineGroupBegin(row_coarse_group); g < FineGroupEnd(row_coarse_group); ++g)
+    {
       numerator += F[g][gp] * phi;
+      if (include_delayed_production)
+        numerator += ComputeDelayedFissionProduction(xs, g, gp) * phi;
+    }
   }
   if (std::fabs(denominator) <= 1.0e-30)
   {
@@ -965,7 +970,11 @@ CMFDAcceleration::CondensedProductionXS(const MultiGroupXS& xs,
       denominator += 1.0;
       for (unsigned int g = FineGroupBegin(row_coarse_group); g < FineGroupEnd(row_coarse_group);
            ++g)
+      {
         numerator += F[g][gp];
+        if (include_delayed_production)
+          numerator += ComputeDelayedFissionProduction(xs, g, gp);
+      }
     }
   }
   return numerator / denominator;
diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/source_functions/source_function.cc b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/source_functions/source_function.cc
index 72b407b80..b1b2369ee 100644
--- a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/source_functions/source_function.cc
+++ b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/source_functions/source_function.cc
@@ -3,6 +3,7 @@
 
 #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/source_functions/source_function.h"
 #include "modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.h"
+#include "modules/linear_boltzmann_solvers/lbs_problem/compute/lbs_compute.h"
 #include "framework/mesh/mesh_continuum/mesh_continuum.h"
 #include "framework/runtime.h"
 #include "framework/logging/log.h"
@@ -123,7 +124,7 @@ SourceFunction::operator()(const LBSGroupset& groupset,
               for (size_t gp = gs_i_; gp <= gs_f_; ++gp)
                 rhs += F_g[gp] * phi_im[gp];
 
-            if (lbs_problem_.GetOptions().use_precursors)
+            if (UseDelayedNeutronProduction(lbs_problem_))
               rhs += DelayedFission(precursors, nu_delayed_sigma_f, &phi[uk_map]);
           }
 
diff --git a/modules/linear_boltzmann_solvers/lbs_problem/compute/lbs_compute.cc b/modules/linear_boltzmann_solvers/lbs_problem/compute/lbs_compute.cc
index 72b3dead7..9b88a40a5 100644
--- a/modules/linear_boltzmann_solvers/lbs_problem/compute/lbs_compute.cc
+++ b/modules/linear_boltzmann_solvers/lbs_problem/compute/lbs_compute.cc
@@ -3,12 +3,40 @@
 
 #include "modules/linear_boltzmann_solvers/lbs_problem/compute/lbs_compute.h"
 #include "modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.h"
+#include "framework/materials/multi_group_xs/multi_group_xs.h"
 #include "framework/runtime.h"
 #include <cassert>
 
 namespace opensn
 {
 
+bool
+UseDelayedNeutronProduction(const LBSProblem& lbs_problem)
+{
+  return lbs_problem.GetNumPrecursors() > 0 and
+         (not lbs_problem.IsTimeDependent() or lbs_problem.GetOptions().use_precursors);
+}
+
+double
+ComputeDelayedFissionProduction(const MultiGroupXS& xs,
+                                const unsigned int to_group,
+                                const unsigned int from_group)
+{
+  if (xs.GetNumPrecursors() == 0)
+    return 0.0;
+
+  const auto& nu_delayed_sigma_f = xs.GetNuDelayedSigmaF();
+  if (nu_delayed_sigma_f.empty())
+    return 0.0;
+
+  double delayed_production = 0.0;
+  for (const auto& precursor : xs.GetPrecursors())
+    delayed_production += precursor.emission_spectrum[to_group] * precursor.fractional_yield *
+                          nu_delayed_sigma_f[from_group];
+
+  return delayed_production;
+}
+
 double
 ComputeFissionProduction(const LBSProblem& lbs_problem, const std::vector<double>& phi)
 {
@@ -16,9 +44,9 @@ ComputeFissionProduction(const LBSProblem& lbs_problem, const std::vector<double
   const auto& grid = lbs_problem.GetGrid();
   const auto& cell_transport_views = lbs_problem.GetCellTransportViews();
   const auto& unit_cell_matrices = lbs_problem.GetUnitCellMatrices();
-  const auto& options = lbs_problem.GetOptions();
   assert(phi.size() == lbs_problem.GetPhiNewLocal().size() &&
          "ComputeFissionProduction size mismatch.");
+  const bool include_delayed_production = UseDelayedNeutronProduction(lbs_problem);
 
   const auto first_grp = 0;
   const auto last_grp = lbs_problem.GetNumGroups() - 1;
@@ -33,7 +61,6 @@ ComputeFissionProduction(const LBSProblem& lbs_problem, const std::vector<double
     // Obtain xs
     const auto& xs = transport_view.GetXS();
     const auto& F = xs.GetProductionMatrix();
-    const auto& nu_delayed_sigma_f = xs.GetNuDelayedSigmaF();
 
     if (not xs.IsFissionable())
       continue;
@@ -52,12 +79,13 @@ ComputeFissionProduction(const LBSProblem& lbs_problem, const std::vector<double
         for (unsigned int gp = 0; gp <= last_grp; ++gp)
           local_production += prod[gp] * phi[uk_map + gp] * IntV_ShapeI;
 
-        // When delayed-neutrons are enabled, this utility assumes production matrix
-        // contributions are prompt-only and adds delayed production separately via
-        // nu_delayed_sigma_f. If the production matrix changes to include delayed
-        // production, adjust this sum to prevent double counting.
-        if (options.use_precursors and xs.GetNumPrecursors() > 0)
-          local_production += nu_delayed_sigma_f[g] * phi[uk_map + g] * IntV_ShapeI;
+        // When delayed-neutron data is included, the production matrix is prompt-only
+        // and delayed production is added separately. If the production matrix changes
+        // to include delayed production, adjust this sum to prevent double counting.
+        if (include_delayed_production)
+          for (unsigned int gp = 0; gp <= last_grp; ++gp)
+            local_production +=
+              ComputeDelayedFissionProduction(xs, g, gp) * phi[uk_map + gp] * IntV_ShapeI;
       }
     } // for node
   } // for cell
diff --git a/modules/linear_boltzmann_solvers/lbs_problem/compute/lbs_compute.h b/modules/linear_boltzmann_solvers/lbs_problem/compute/lbs_compute.h
index f0ed9740f..650eb33a2 100644
--- a/modules/linear_boltzmann_solvers/lbs_problem/compute/lbs_compute.h
+++ b/modules/linear_boltzmann_solvers/lbs_problem/compute/lbs_compute.h
@@ -9,6 +9,15 @@ namespace opensn
 {
 
 class LBSProblem;
+class MultiGroupXS;
+
+/// Returns true when delayed neutron production contributes to transport sources.
+bool UseDelayedNeutronProduction(const LBSProblem& lbs_problem);
+
+/// Computes the delayed fission production from a source group to a destination group.
+double ComputeDelayedFissionProduction(const MultiGroupXS& xs,
+                                       unsigned int to_group,
+                                       unsigned int from_group);
 
 /**
  * Compute the total fission production in the problem.
diff --git a/modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.cc b/modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.cc
index 3e1db9361..92eb1bb4f 100644
--- a/modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.cc
+++ b/modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.cc
@@ -2,6 +2,7 @@
 // SPDX-License-Identifier: MIT
 
 #include "modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.h"
+#include "modules/linear_boltzmann_solvers/lbs_problem/compute/lbs_compute.h"
 #include "modules/linear_boltzmann_solvers/lbs_problem/point_source/point_source.h"
 #include "modules/linear_boltzmann_solvers/lbs_problem/groupset/lbs_groupset.h"
 #include "framework/field_functions/field_function_grid_based.h"
@@ -535,7 +536,11 @@ LBSProblem::GetOptionsBlock()
   params.AddOptionalParameter("write_restart_time_interval",
                               0,
                               "Time interval in seconds at which restart data is to be written.");
-  params.AddOptionalParameter("use_precursors", true, "Flag for using delayed neutron precursors.");
+  params.AddOptionalParameter(
+    "use_precursors",
+    true,
+    "Flag for storing and evolving delayed neutron precursors. Steady-state fission sources "
+    "include delayed neutron production whenever delayed data exists.");
   params.AddOptionalParameter("use_source_moments",
                               false,
                               "Flag for ignoring fixed sources and selectively using source "
@@ -986,17 +991,18 @@ LBSProblem::InitializeMaterials()
                          "precursor coupling.";
   }
 
-  // check compatibility when at least one fissionable material has delayed-neutron data
-  if (options_.use_precursors and has_fissionable_precursors)
+  // Check compatibility when delayed-neutron production is active and at least one
+  // fissionable material has delayed-neutron data.
+  if (UseDelayedNeutronProduction(*this) and has_fissionable_precursors)
   {
     for (const auto& [mat_id, xs] : block_id_to_xs_map_)
     {
       OpenSnInvalidArgumentIf(xs->IsFissionable() and xs->GetNumPrecursors() == 0,
                               GetName() + ": incompatible cross-section data for material id " +
                                 std::to_string(mat_id) +
-                                ". When options.use_precursors=true and "
-                                "delayed-neutron precursor data is present for one fissionable "
-                                "material, it must be present for all fissionable materials.");
+                                ". When delayed-neutron production is active and delayed-neutron "
+                                "precursor data is present for one fissionable material, it must "
+                                "be present for all fissionable materials.");
     }
   }
 
diff --git a/python/lib/solver.cc b/python/lib/solver.cc
index 6c7fc3ced..2515399d9 100644
--- a/python/lib/solver.cc
+++ b/python/lib/solver.cc
@@ -797,13 +797,15 @@ WrapLBS(py::module& slv)
           - write_restart_time_interval: int, default=0
             (must be 0 or >=30)
           - use_precursors: bool, default=True
-            Enable delayed-neutron precursor treatment. This is treated as user intent and remains
-            active across later ``SetXSMap`` calls, even if the current XS map temporarily contains
-            no precursor-bearing material. When XS maps are swapped, existing precursor
+            Enable delayed-neutron precursor storage and transient precursor evolution. Steady-state
+            fission sources include delayed neutron production whenever the active XS map contains
+            delayed data, independent of this option. This option is treated as user intent and
+            remains active across later ``SetXSMap`` calls, even if the current XS map temporarily
+            contains no precursor-bearing material. When XS maps are swapped, existing precursor
             concentrations are remapped by cell and precursor-family index; new families start at
-            zero and removed families are discarded. If any fissionable material in the active map
-            has precursor data and ``use_precursors=True``, all fissionable materials in that map
-            must have precursor data. Non-fissionable materials may have zero precursors.
+            zero and removed families are discarded. If delayed-neutron production is active and any
+            fissionable material in the active map has precursor data, all fissionable materials in
+            that map must have precursor data. Non-fissionable materials may have zero precursors.
           - use_source_moments: bool, default=False
           - save_angular_flux: bool, default=False
             Store angular flux state (`psi`) for transient mode, angular-flux
@@ -1282,11 +1284,12 @@ WrapLBS(py::module& slv)
               ``.restart.h5`` to this stem.
           - write_restart_time_interval: int, default=0
           - use_precursors: bool, default=True
-            Enable delayed-neutron precursor treatment. This remains active across later
-            ``SetXSMap`` calls, even if the current map temporarily has no precursor-bearing
-            material. When XS maps change, precursor concentrations are remapped by cell and
-            precursor-family index; newly introduced families start at zero and removed families
-            are discarded.
+            Enable delayed-neutron precursor storage and transient precursor evolution. Steady-state
+            fission sources include delayed neutron production whenever the active XS map contains
+            delayed data, independent of this option. This remains active across later ``SetXSMap``
+            calls, even if the current map temporarily has no precursor-bearing material. When XS
+            maps change, precursor concentrations are remapped by cell and precursor-family index;
+            newly introduced families start at zero and removed families are discarded.
           - use_source_moments: bool, default=False
           - save_angular_flux: bool, default=False
             Store angular flux state (`psi`) for transient mode, angular-flux
diff --git a/python/lib/xs.cc b/python/lib/xs.cc
index c3ae3c362..dcf4cf6c7 100644
--- a/python/lib/xs.cc
+++ b/python/lib/xs.cc
@@ -182,7 +182,16 @@ WrapMultiGroupXS(py::module& xs)
 
     Notes
     -----
-    This method mutates ``self`` by replacing its current contents.
+    When delayed-neutron data is present in the OpenMC MGXS library, this
+    method also imports:
+
+    - the number of precursor groups
+    - prompt and delayed fission production data
+    - precursor decay constants
+    - delayed emission spectra
+
+    If delayed data is present but ``chi-prompt`` is not available, the steady
+    fission spectrum ``chi`` is used as the prompt spectrum.
     )",
     py::arg("file_name"),
     py::arg("dataset_name"),
diff --git a/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_1d_1g_cmfd.py b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_1d_1g_cmfd.py
index 104a5252b..910e89fdc 100644
--- a/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_1d_1g_cmfd.py
+++ b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_1d_1g_cmfd.py
@@ -4,7 +4,7 @@
 """
 1D 1G k-eigenvalue smoke test using power iteration with CMFD acceleration.
 
-Test: Final k-eigenvalue: 0.9945456
+Test: Final k-eigenvalue: 0.9995433
 """
 
 import os
@@ -148,7 +148,7 @@ def run_cmfd():
         print(f"Python CMFD sweep speedup: {cmfd_speedup}")
 
     if not get_option("cmfd_limiter_stress", False):
-        expected_k = 0.9945456
+        expected_k = 0.9995433
         if abs(pi_k - expected_k) > 5.0e-6:
             raise RuntimeError(f"Expected PI k={expected_k}, got {pi_k}")
         if abs(cmfd_k - expected_k) > 5.0e-6:
diff --git a/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_1d_1g_precursor_toggle.py b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_1d_1g_precursor_toggle.py
new file mode 100644
index 000000000..c032d5a75
--- /dev/null
+++ b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_1d_1g_precursor_toggle.py
@@ -0,0 +1,96 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+"""
+1D 1G delayed k-eigenvalue test.
+
+Steady-state k should include delayed neutron production from the XS data
+whether or not precursor storage is enabled.
+"""
+
+import os
+import sys
+
+if "opensn_console" not in globals():
+    from mpi4py import MPI
+    size = MPI.COMM_WORLD.size
+    rank = MPI.COMM_WORLD.rank
+    sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), "../../../../../")))
+    from pyopensn.mesh import OrthogonalMeshGenerator
+    from pyopensn.xs import MultiGroupXS
+    from pyopensn.aquad import GLProductQuadrature1DSlab
+    from pyopensn.solver import DiscreteOrdinatesProblem, PowerIterationKEigenSolver
+
+
+def run_case(use_precursors):
+    length = 100.0
+    num_cells = 50
+    dx = length / num_cells
+    nodes = [i * dx for i in range(num_cells + 1)]
+    meshgen = OrthogonalMeshGenerator(node_sets=[nodes])
+    grid = meshgen.Execute()
+    grid.SetUniformBlockID(0)
+
+    xs = MultiGroupXS()
+    xs.LoadFromOpenSn("../../../../assets/xs/simple_fissile.xs")
+
+    problem = DiscreteOrdinatesProblem(
+        mesh=grid,
+        num_groups=1,
+        groupsets=[
+            {
+                "groups_from_to": (0, 0),
+                "angular_quadrature": GLProductQuadrature1DSlab(
+                    n_polar=32,
+                    scattering_order=0,
+                ),
+                "inner_linear_method": "petsc_gmres",
+                "l_max_its": 200,
+                "l_abs_tol": 1.0e-10,
+            }
+        ],
+        xs_map=[
+            {
+                "block_ids": [0],
+                "xs": xs,
+            }
+        ],
+        options={
+            "use_precursors": use_precursors,
+            "verbose_inner_iterations": False,
+            "verbose_outer_iterations": False,
+        },
+    )
+
+    solver = PowerIterationKEigenSolver(
+        problem=problem,
+        max_iters=2000,
+        k_tol=1.0e-8,
+    )
+    solver.Initialize()
+    solver.Execute()
+    return solver.GetEigenvalue()
+
+
+if __name__ == "__main__":
+    expected_num_procs = 4
+    if size != expected_num_procs:
+        sys.exit(
+            f"Incorrect number of processors. Expected {expected_num_procs} processors "
+            f"but got {size}."
+        )
+
+    k_without_precursors = run_case(False)
+    k_with_precursors = run_case(True)
+    k_difference = abs(k_with_precursors - k_without_precursors)
+
+    if rank == 0:
+        print(f"Python use_precursors=False k-eigenvalue: {k_without_precursors}")
+        print(f"Python use_precursors=True k-eigenvalue: {k_with_precursors}")
+        print(f"Python precursor toggle k-difference: {k_difference}")
+
+    if k_difference > 1.0e-10:
+        raise RuntimeError(
+            "Expected steady-state delayed k-eigenvalue to be independent of "
+            f"use_precursors, difference is {k_difference}."
+        )
diff --git a/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_1d_openmcxs_delayed_slab.py b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_1d_openmcxs_delayed_slab.py
new file mode 100644
index 000000000..7bb8cc32b
--- /dev/null
+++ b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_1d_openmcxs_delayed_slab.py
@@ -0,0 +1,145 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+"""
+1D slab k-eigenvalue solve using OpenMC-generated delayed-neutron MGXS data.
+
+Side-by-side OpenMC/OpenSn comparison on a simple reflected-slab geometry used
+to generate ``slab_delayed_mgxs.h5``.
+
+OpenMC model:
+- fuel slab on [-2, 2]
+- water reflector on [2, 10]
+- vacuum boundaries on the slab ends
+- reflective transverse directions on the OpenMC side, so the problem is
+  effectively 1D
+- MG k-effective = 0.3587+/-0.0004
+
+The script prints:
+- imported delayed-data summaries for each material
+- OpenSn k_eff
+- a global balance table
+- per-group scalar-flux maxima in the fuel and water regions
+"""
+
+import os
+import sys
+
+if "opensn_console" not in globals():
+    from mpi4py import MPI
+    size = MPI.COMM_WORLD.size
+    rank = MPI.COMM_WORLD.rank
+    sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), "../../../../../")))
+    from pyopensn.mesh import OrthogonalMeshGenerator
+    from pyopensn.logvol import RPPLogicalVolume
+    from pyopensn.xs import MultiGroupXS
+    from pyopensn.aquad import GLProductQuadrature1DSlab
+    from pyopensn.fieldfunc import FieldFunctionInterpolationVolume
+    from pyopensn.solver import DiscreteOrdinatesProblem, PowerIterationKEigenSolver
+
+
+def xs_summary(name, xs):
+    nu_sigma_f = list(xs.nu_sigma_f) if xs.is_fissionable else []
+    nu_prompt_sigma_f = list(xs.nu_prompt_sigma_f) if xs.is_fissionable else []
+    nu_delayed_sigma_f = list(xs.nu_delayed_sigma_f) if xs.is_fissionable else []
+    return {
+        "dataset": name,
+        "is_fissionable": xs.is_fissionable,
+        "num_groups": xs.num_groups,
+        "num_precursors": xs.num_precursors,
+        "sum_nu_sigma_f": sum(nu_sigma_f),
+        "sum_nu_prompt_sigma_f": sum(nu_prompt_sigma_f),
+        "sum_nu_delayed_sigma_f": sum(nu_delayed_sigma_f),
+    }
+
+
+def region_flux_max(field_function, logical_volume):
+    ffi = FieldFunctionInterpolationVolume()
+    ffi.SetOperationType("max")
+    ffi.SetLogicalVolume(logical_volume)
+    ffi.AddFieldFunction(field_function)
+    ffi.Execute()
+    return ffi.GetValue()
+
+
+if __name__ == "__main__":
+    xs_fuel = MultiGroupXS()
+    xs_fuel.LoadFromOpenMC("slab_delayed_mgxs.h5", "fuel", 294.0)
+
+    xs_water = MultiGroupXS()
+    xs_water.LoadFromOpenMC("slab_delayed_mgxs.h5", "water", 294.0)
+
+    if xs_fuel.num_groups != xs_water.num_groups:
+        raise RuntimeError("Fuel and water MGXS datasets must have the same number of groups.")
+
+    slab_min = -2.0
+    fuel_max = 2.0
+    slab_max = 10.0
+    n_cells = 120
+    dz = (slab_max - slab_min) / n_cells
+    nodes = [slab_min + i * dz for i in range(n_cells + 1)]
+
+    meshgen = OrthogonalMeshGenerator(node_sets=[nodes])
+    grid = meshgen.Execute()
+    grid.SetUniformBlockID(1)
+
+    fuel_region = RPPLogicalVolume(infx=True, infy=True, zmin=slab_min, zmax=fuel_max)
+    water_region = RPPLogicalVolume(infx=True, infy=True, zmin=fuel_max, zmax=slab_max)
+    grid.SetBlockIDFromLogicalVolume(fuel_region, 0, True)
+
+    num_groups = xs_fuel.num_groups
+    pquad = GLProductQuadrature1DSlab(
+        n_polar=8,
+        scattering_order=max(xs_fuel.scattering_order, xs_water.scattering_order),
+    )
+
+    problem = DiscreteOrdinatesProblem(
+        mesh=grid,
+        num_groups=num_groups,
+        groupsets=[
+            {
+                "groups_from_to": (0, num_groups - 1),
+                "angular_quadrature": pquad,
+                "inner_linear_method": "petsc_gmres",
+                "l_abs_tol": 1.0e-8,
+                "l_max_its": 500,
+            },
+        ],
+        xs_map=[
+            {"block_ids": [0], "xs": xs_fuel},
+            {"block_ids": [1], "xs": xs_water},
+        ],
+        boundary_conditions=[
+            {"name": "zmin", "type": "vacuum"},
+            {"name": "zmax", "type": "vacuum"},
+        ],
+        options={
+            "use_precursors": True,
+            "verbose_inner_iterations": True,
+            "verbose_outer_iterations": True,
+        },
+    )
+
+    solver = PowerIterationKEigenSolver(
+        problem=problem,
+        max_iters=200,
+        k_tol=1.0e-8,
+    )
+    solver.Initialize()
+    solver.Execute()
+
+    balance = solver.ComputeBalanceTable()
+    scalar_flux = problem.GetScalarFluxFieldFunction()
+
+    if rank == 0:
+        print("FUEL_XS_SUMMARY", xs_summary("fuel", xs_fuel))
+        print("WATER_XS_SUMMARY", xs_summary("water", xs_water))
+        print(f"OPENSN_KEFF {solver.GetEigenvalue():.8f}")
+        print(f"OPENSN_BALANCE {dict(balance)}")
+
+    for g in range(num_groups):
+        fuel_phi_max = region_flux_max(scalar_flux[g], fuel_region)
+        water_phi_max = region_flux_max(scalar_flux[g], water_region)
+        if rank == 0:
+            print(f"GROUP_{g}_FUEL_PHI_MAX {fuel_phi_max:.10e}")
+            print(f"GROUP_{g}_WATER_PHI_MAX {water_phi_max:.10e}")
diff --git a/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_1g_reflecting_cmfd_agg.py b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_1g_reflecting_cmfd_agg.py
index 5fd44e6d5..b51f34ea0 100644
--- a/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_1g_reflecting_cmfd_agg.py
+++ b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_1g_reflecting_cmfd_agg.py
@@ -4,7 +4,7 @@
 """
 3D 1G reflecting-boundary k-eigenvalue test using aggregated CMFD acceleration.
 
-Test: Final k-eigenvalue: 0.995
+Test: Final k-eigenvalue: 1.0
 """
 
 import os
@@ -101,7 +101,7 @@
         print(f"Python sweeps: {sweeps}")
         print(f"Balance={balance['balance']:.6e}")
 
-    if abs(k - 0.995) > 1.0e-4:
-        raise RuntimeError(f"Expected k near 0.995, got {k}")
+    if abs(k - 1.0) > 1.0e-4:
+        raise RuntimeError(f"Expected k near 1.0, got {k}")
     if abs(balance["balance"]) > 5.0e-4:
         raise RuntimeError(f"Expected balance closure, got {balance['balance']}")
diff --git a/test/python/modules/linear_boltzmann_solvers/transport_keigen/slab_delayed_mgxs.h5 b/test/python/modules/linear_boltzmann_solvers/transport_keigen/slab_delayed_mgxs.h5
new file mode 100644
index 000000000..ed250972c
Binary files /dev/null and b/test/python/modules/linear_boltzmann_solvers/transport_keigen/slab_delayed_mgxs.h5 differ
diff --git a/test/python/modules/linear_boltzmann_solvers/transport_keigen/tests.json b/test/python/modules/linear_boltzmann_solvers/transport_keigen/tests.json
index 0b8a8ec98..01f80bd55 100644
--- a/test/python/modules/linear_boltzmann_solvers/transport_keigen/tests.json
+++ b/test/python/modules/linear_boltzmann_solvers/transport_keigen/tests.json
@@ -8,11 +8,36 @@
         "type": "FloatCompare",
         "key": "NLKE final",
         "wordnum": 13,
-        "gold": 0.99454,
+        "gold": 0.99954,
         "abs_tol": 1e-05
       }
     ]
   },
+  {
+    "file": "keigenvalue_transport_1d_1g_precursor_toggle.py",
+    "comment": "1D delayed k-eigenvalue is independent of precursor storage in steady state",
+    "num_procs": 4,
+    "checks": [
+      {
+        "type": "KeyValuePair",
+        "key": "Python use_precursors=False k-eigenvalue:",
+        "goldvalue": 0.9995433,
+        "abs_tol": 1e-05
+      },
+      {
+        "type": "KeyValuePair",
+        "key": "Python use_precursors=True k-eigenvalue:",
+        "goldvalue": 0.9995433,
+        "abs_tol": 1e-05
+      },
+      {
+        "type": "KeyValuePair",
+        "key": "Python precursor toggle k-difference:",
+        "goldvalue": 0.0,
+        "abs_tol": 1e-10
+      }
+    ]
+  },
   {
     "file": "keigenvalue_transport_1d_1g_cmfd.py",
     "comment": "1D keigenvalue test using power iteration with aggregated CMFD acceleration",
@@ -21,13 +46,13 @@
       {
         "type": "KeyValuePair",
         "key": "Python PI k-eigenvalue:",
-        "goldvalue": 0.9945456,
+        "goldvalue": 0.9995433,
         "abs_tol": 5e-06
       },
       {
         "type": "KeyValuePair",
         "key": "Python CMFD k-eigenvalue:",
-        "goldvalue": 0.9945456,
+        "goldvalue": 0.9995433,
         "abs_tol": 5e-06
       },
       {
@@ -49,13 +74,13 @@
       {
         "type": "KeyValuePair",
         "key": "Python PI k-eigenvalue:",
-        "goldvalue": 0.9945456,
+        "goldvalue": 0.9995433,
         "abs_tol": 5e-06
       },
       {
         "type": "KeyValuePair",
         "key": "Python CMFD k-eigenvalue:",
-        "goldvalue": 0.9945456,
+        "goldvalue": 0.9995433,
         "abs_tol": 5e-06
       },
       {
@@ -306,7 +331,7 @@
         "type": "FloatCompare",
         "key": "NLKE final",
         "wordnum": 13,
-        "gold": 0.99454,
+        "gold": 0.99954,
         "abs_tol": 1e-05
       }
     ]
@@ -386,6 +411,19 @@
       }
     ]
   },
+  {
+    "file": "keigenvalue_transport_1d_openmcxs_delayed_slab.py",
+    "comment": "1D 2G delayed-neutron slab keigenvalue test using OpenMC MGXS library",
+    "num_procs": 4,
+    "checks": [
+      {
+        "type": "KeyValuePair",
+        "key": "OPENSN_KEFF",
+        "goldvalue": 0.35819387,
+        "abs_tol": 1e-05
+      }
+    ]
+  },
   {
     "file": "c5g7_restart.py",
     "outfileprefix": "c5g7_restart",
@@ -538,14 +576,14 @@
         "type": "FloatCompare",
         "key": "PI final",
         "wordnum": 13,
-        "gold": 0.995,
+        "gold": 1.0,
         "abs_tol": 0.0001
       },
       {
         "type": "KeyValuePair",
         "key": "Balance=",
         "goldvalue": 0.0,
-        "abs_tol": 1e-08
+        "abs_tol": 1e-07
       }
     ]
   },
@@ -557,7 +595,7 @@
       {
         "type": "KeyValuePair",
         "key": "Python k-eigenvalue:",
-        "goldvalue": 0.995,
+        "goldvalue": 1.0,
         "abs_tol": 0.0001
       },
       {
@@ -579,7 +617,7 @@
       {
         "type": "KeyValuePair",
         "key": "Python k-eigenvalue:",
-        "goldvalue": 0.995,
+        "goldvalue": 1.0,
         "abs_tol": 0.0001
       },
       {
diff --git a/test/python/modules/linear_boltzmann_solvers/transport_transient/tests.json b/test/python/modules/linear_boltzmann_solvers/transport_transient/tests.json
index 4b1fef5de..7e04989bb 100644
--- a/test/python/modules/linear_boltzmann_solvers/transport_transient/tests.json
+++ b/test/python/modules/linear_boltzmann_solvers/transport_transient/tests.json
@@ -528,6 +528,20 @@
       }
     ]
   },
+  {
+    "file": "transient_keigen_1d_delayed_prke_vs_stk_openmc.py",
+    "comment": "1D delayed homogeneous step: PRKE vs space-time kinetics with OpenMC cross sections",
+    "num_procs": 4,
+    "checks": [
+      {
+        "type": "FloatCompare",
+        "key": "PRKE_STK_PASS",
+        "wordnum": 1,
+        "gold": 1,
+        "abs_tol": 1e-12
+      }
+    ]
+  },
   {
     "file": "transient_keigen_1d_delayed_fission_prod_count.py",
     "comment": "1D delayed fission production invariant to precursor count (1p vs 2p)",
@@ -598,6 +612,20 @@
       }
     ]
   },
+  {
+    "file": "transient_keigen_3d_delayed_prke_vs_stk_2p_openmc.py",
+    "comment": "3D delayed homogeneous step: PRKE vs space-time kinetics with OpenMC cross sections",
+    "num_procs": 4,
+    "checks": [
+      {
+        "type": "FloatCompare",
+        "key": "PRKE_STK_PASS",
+        "wordnum": 1,
+        "gold": 1,
+        "abs_tol": 1e-12
+      }
+    ]
+  },
   {
     "file": "transient_keigen_3d_delayed_prke_vs_stk_2p_callbacks.py",
     "comment": "3D delayed homogeneous step: PRKE vs space-time kinetics",
diff --git a/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_keigen_1d_delayed_prke_vs_stk_openmc.py b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_keigen_1d_delayed_prke_vs_stk_openmc.py
new file mode 100644
index 000000000..c135ed9c0
--- /dev/null
+++ b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_keigen_1d_delayed_prke_vs_stk_openmc.py
@@ -0,0 +1,170 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+"""
+1D delayed transient: OpenMC XS homogeneous step vs PRKE.
+
+Validate space-time kinetics against point-reactor kinetics (PRKE) for a
+homogeneous perturbation in a homogeneous system.
+
+1-group, 1 precursor, reflecting boundaries (infinite-medium). A step
+to a supercritical xs is applied at t=0. Space-time kinetics should follow
+PRKE for this homogeneous case.
+
+Transport cross sections are loaded from OpenMC MGXS HDF5 files. The PRKE
+reference values use the same benchmark precursor constants as the original
+OpenSn XS source.
+
+PRKE_STK_PASS is 1 if the space-time fission-production ratio matches PRKE within 2%
+for t<=0.2.
+"""
+
+import math
+import os
+import sys
+
+ASSET_XS_DIR = os.path.abspath(os.path.join(os.path.dirname(__file__), "../../../../assets/xs"))
+
+if "opensn_console" not in globals():
+    from mpi4py import MPI
+    size = MPI.COMM_WORLD.size
+    rank = MPI.COMM_WORLD.rank
+    sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), "../../../../../")))
+    from pyopensn.solver import (
+        DiscreteOrdinatesProblem,
+        PowerIterationKEigenSolver,
+        TransientSolver,
+    )
+    from pyopensn.aquad import GLProductQuadrature1DSlab
+    from pyopensn.xs import MultiGroupXS
+    from pyopensn.mesh import OrthogonalMeshGenerator
+
+
+def read_precursor_value(path, block_name):
+    begin = f"{block_name}_BEGIN"
+    end = f"{block_name}_END"
+    in_block = False
+    with open(path, "r", encoding="utf-8") as handle:
+        for line in handle:
+            line = line.strip()
+            if not line or line.startswith("#"):
+                continue
+            if line == begin:
+                in_block = True
+                continue
+            if line == end:
+                in_block = False
+                continue
+            if in_block:
+                parts = line.split()
+                if len(parts) >= 2:
+                    return float(parts[1])
+    raise RuntimeError(f"Failed to find {block_name} in {path}")
+
+
+def prke_phi_ratio(t, beta, lam, rho, Lambda):
+    a = (rho - beta) / Lambda - lam
+    b = math.sqrt(((rho - beta) / Lambda + lam) ** 2 + 4.0 * beta * lam / Lambda)
+    w1 = 0.5 * (a + b)
+    w2 = 0.5 * (a - b)
+
+    k1 = (beta / Lambda) / (w1 + lam)
+    k2 = (beta / Lambda) / (w2 + lam)
+
+    c2 = (beta / (Lambda * lam) - k1) / (k2 - k1)
+    c1 = 1.0 - c2
+
+    return c1 * math.exp(w1 * t) + c2 * math.exp(w2 * t)
+
+
+if __name__ == "__main__":
+    dx = 8.0 / 40
+    nodes = [i * dx for i in range(40 + 1)]
+    meshgen = OrthogonalMeshGenerator(node_sets=[nodes])
+    grid = meshgen.Execute()
+    grid.SetUniformBlockID(0)
+
+    xs_crit = MultiGroupXS()
+    xs_crit.LoadFromOpenMC("xs1g_delayed_crit_1p.h5", "set1", 294.0)
+
+    xs_super = MultiGroupXS()
+    xs_super.LoadFromOpenMC("xs1g_delayed_super_1p.h5", "set1", 294.0)
+
+    pquad = GLProductQuadrature1DSlab(n_polar=4, scattering_order=0)
+
+    num_groups = 1
+    phys = DiscreteOrdinatesProblem(
+        mesh=grid,
+        num_groups=num_groups,
+        groupsets=[
+            {
+                "groups_from_to": (0, num_groups - 1),
+                "angular_quadrature": pquad,
+                "inner_linear_method": "petsc_richardson",
+                "l_abs_tol": 1.0e-8,
+                "l_max_its": 200,
+            },
+        ],
+        xs_map=[{"block_ids": [0], "xs": xs_crit}],
+        boundary_conditions=[
+            {"name": "zmin", "type": "reflecting"},
+            {"name": "zmax", "type": "reflecting"},
+        ],
+        options={
+            "save_angular_flux": True,
+            "use_precursors": True,
+            "verbose_inner_iterations": False,
+            "verbose_outer_iterations": False,
+        },
+    )
+
+    keigen = PowerIterationKEigenSolver(problem=phys, max_iters=200, k_tol=1.0e-10)
+    keigen.Initialize()
+    keigen.Execute()
+
+    phys.SetTimeDependentMode()
+
+    solver = TransientSolver(problem=phys, initial_state="existing")
+    solver.Initialize()
+
+    phys.SetXSMap(xs_map=[{"block_ids": [0], "xs": xs_super}])
+
+    sigma_a = xs_super.sigma_a[0]
+    nu_sigma_f = xs_super.nu_sigma_f[0]
+    inv_vel = xs_super.inv_velocity[0]
+    v = 1.0 / inv_vel
+
+    k_eff = nu_sigma_f / sigma_a
+    rho = (k_eff - 1.0) / k_eff
+
+    beta = read_precursor_value(
+        os.path.join(ASSET_XS_DIR, "xs1g_delayed_super_1p.cxs"),
+        "PRECURSOR_FRACTIONAL_YIELDS",
+    )
+    lam = read_precursor_value(
+        os.path.join(ASSET_XS_DIR, "xs1g_delayed_super_1p.cxs"),
+        "PRECURSOR_DECAY_CONSTANTS",
+    )
+    Lambda = 1.0 / (v * nu_sigma_f)
+
+    dt = 1.0e-2
+    solver.SetTimeStep(dt)
+    solver.SetTheta(1.0)
+
+    fp0 = phys.ComputeFissionProduction("new")
+
+    t_end = 0.2
+    rel_tol = 2.0e-2
+    ok = True
+    while phys.GetTime() < t_end:
+        solver.Advance()
+        fp_new = phys.ComputeFissionProduction("new")
+        t_to = phys.GetTime()
+
+        ratio_num = fp_new / fp0
+        ratio_prke = prke_phi_ratio(t_to, beta, lam, rho, Lambda)
+        if abs(ratio_num - ratio_prke) > rel_tol * ratio_prke:
+            ok = False
+
+    if rank == 0:
+        print(f"PRKE_STK_PASS {1 if ok else 0}")
diff --git a/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_keigen_3d_delayed_prke_vs_stk_2p_openmc.py b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_keigen_3d_delayed_prke_vs_stk_2p_openmc.py
new file mode 100644
index 000000000..5a30b244f
--- /dev/null
+++ b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_keigen_3d_delayed_prke_vs_stk_2p_openmc.py
@@ -0,0 +1,217 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+"""
+3D delayed transient: OpenMC XS homogeneous step vs PRKE (2 precursors).
+
+Validate space-time kinetics against point-reactor kinetics (PRKE) for a
+homogeneous perturbation in a homogeneous system with two precursors.
+
+1-group, 2 precursors, reflecting boundaries (infinite-medium). A step
+to a supercritical xs is applied at t=0. Space-time kinetics should follow
+PRKE for this homogeneous case.
+
+Transport cross sections are loaded from OpenMC MGXS HDF5 files. The PRKE
+reference values use the same benchmark precursor constants as the original
+OpenSn XS source.
+
+PRKE_STK_PASS is 1 if the space-time fission-production ratio matches PRKE within 2%
+for t<=0.15.
+"""
+
+import os
+import sys
+
+ASSET_XS_DIR = os.path.abspath(os.path.join(os.path.dirname(__file__), "../../../../assets/xs"))
+
+if "opensn_console" not in globals():
+    from mpi4py import MPI
+    size = MPI.COMM_WORLD.size
+    rank = MPI.COMM_WORLD.rank
+    sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), "../../../../../")))
+    from pyopensn.solver import (
+        DiscreteOrdinatesProblem,
+        PowerIterationKEigenSolver,
+        TransientSolver,
+    )
+    from pyopensn.aquad import GLCProductQuadrature3DXYZ
+    from pyopensn.xs import MultiGroupXS
+    from pyopensn.mesh import OrthogonalMeshGenerator
+
+
+def read_block_values(path, block_name):
+    begin = f"{block_name}_BEGIN"
+    end = f"{block_name}_END"
+    values = []
+    in_block = False
+    with open(path, "r", encoding="utf-8") as handle:
+        for line in handle:
+            line = line.strip()
+            if not line or line.startswith("#"):
+                continue
+            if line == begin:
+                in_block = True
+                continue
+            if line == end:
+                in_block = False
+                continue
+            if in_block:
+                parts = line.split()
+                if len(parts) >= 2:
+                    values.append(float(parts[1]))
+    if not values:
+        raise RuntimeError(f"Failed to find {block_name} in {path}")
+    return values
+
+
+def solve_linear(A, b):
+    n = len(b)
+    a = [row[:] for row in A]
+    x = b[:]
+    for i in range(n):
+        pivot = i
+        for r in range(i + 1, n):
+            if abs(a[r][i]) > abs(a[pivot][i]):
+                pivot = r
+        if abs(a[pivot][i]) < 1.0e-14:
+            raise RuntimeError("Singular system in PRKE solve")
+        if pivot != i:
+            a[i], a[pivot] = a[pivot], a[i]
+            x[i], x[pivot] = x[pivot], x[i]
+        piv = a[i][i]
+        for j in range(i, n):
+            a[i][j] /= piv
+        x[i] /= piv
+        for r in range(n):
+            if r == i:
+                continue
+            factor = a[r][i]
+            if factor == 0.0:
+                continue
+            for j in range(i, n):
+                a[r][j] -= factor * a[i][j]
+            x[r] -= factor * x[i]
+    return x
+
+
+def prke_step(phi, C, dt, beta, lambdas, rho, Lambda):
+    m = len(lambdas)
+    beta_total = sum(beta)
+    size = 1 + m
+    A = [[0.0 for _ in range(size)] for _ in range(size)]
+    b = [0.0 for _ in range(size)]
+
+    A[0][0] = 1.0 - dt * (rho - beta_total) / Lambda
+    for i in range(m):
+        A[0][1 + i] = -dt * lambdas[i]
+
+    b[0] = phi
+
+    for i in range(m):
+        A[1 + i][0] = -dt * (beta[i] / Lambda)
+        A[1 + i][1 + i] = 1.0 + dt * lambdas[i]
+        b[1 + i] = C[i]
+
+    x = solve_linear(A, b)
+    return x[0], x[1:]
+
+
+if __name__ == "__main__":
+    dx = 8.0 / 4
+    nodes = [i * dx for i in range(4 + 1)]
+    meshgen = OrthogonalMeshGenerator(node_sets=[nodes, nodes, nodes])
+    grid = meshgen.Execute()
+    grid.SetUniformBlockID(0)
+
+    xs_crit = MultiGroupXS()
+    xs_crit.LoadFromOpenMC("xs1g_delayed_crit_2p.h5", "set1", 294.0)
+
+    xs_super = MultiGroupXS()
+    xs_super.LoadFromOpenMC("xs1g_delayed_super_2p.h5", "set1", 294.0)
+
+    pquad = GLCProductQuadrature3DXYZ(n_polar=2, n_azimuthal=4, scattering_order=0)
+
+    num_groups = 1
+    phys = DiscreteOrdinatesProblem(
+        mesh=grid,
+        num_groups=num_groups,
+        groupsets=[
+            {
+                "groups_from_to": (0, num_groups - 1),
+                "angular_quadrature": pquad,
+                "inner_linear_method": "petsc_gmres",
+                "l_abs_tol": 1.0e-8,
+                "l_max_its": 200,
+                "gmres_restart_interval": 10,
+            },
+        ],
+        xs_map=[{"block_ids": [0], "xs": xs_crit}],
+        boundary_conditions=[
+            {"name": "xmin", "type": "reflecting"},
+            {"name": "xmax", "type": "reflecting"},
+            {"name": "ymin", "type": "reflecting"},
+            {"name": "ymax", "type": "reflecting"},
+            {"name": "zmin", "type": "reflecting"},
+            {"name": "zmax", "type": "reflecting"},
+        ],
+        options={
+            "save_angular_flux": True,
+            "use_precursors": True,
+            "verbose_inner_iterations": False,
+            "verbose_outer_iterations": False,
+        },
+    )
+
+    keigen = PowerIterationKEigenSolver(problem=phys, max_iters=200, k_tol=1.0e-10)
+    keigen.Initialize()
+    keigen.Execute()
+
+    phys.SetTimeDependentMode()
+
+    solver = TransientSolver(problem=phys, initial_state="existing")
+    solver.Initialize()
+
+    phys.SetXSMap(xs_map=[{"block_ids": [0], "xs": xs_super}])
+
+    sigma_a = xs_super.sigma_a[0]
+    nu_sigma_f = xs_super.nu_sigma_f[0]
+    inv_vel = xs_super.inv_velocity[0]
+    v = 1.0 / inv_vel
+
+    k_eff = nu_sigma_f / sigma_a
+    rho = (k_eff - 1.0) / k_eff
+    Lambda = 1.0 / (v * nu_sigma_f)
+
+    beta = read_block_values(
+        os.path.join(ASSET_XS_DIR, "xs1g_delayed_super_2p.cxs"),
+        "PRECURSOR_FRACTIONAL_YIELDS",
+    )
+    lambdas = read_block_values(
+        os.path.join(ASSET_XS_DIR, "xs1g_delayed_super_2p.cxs"),
+        "PRECURSOR_DECAY_CONSTANTS",
+    )
+
+    dt = 1.0e-2
+    solver.SetTimeStep(dt)
+    solver.SetTheta(1.0)
+
+    fp0 = phys.ComputeFissionProduction("new")
+
+    phi = 1.0
+    C = [beta[i] / (Lambda * lambdas[i]) for i in range(len(lambdas))]
+
+    t_end = 0.15
+    rel_tol = 2.0e-2
+    ok = True
+    while phys.GetTime() < t_end:
+        solver.Advance()
+        fp_new = phys.ComputeFissionProduction("new")
+
+        phi, C = prke_step(phi, C, dt, beta, lambdas, rho, Lambda)
+        ratio_num = fp_new / fp0
+        ratio_prke = phi
+        if abs(ratio_num - ratio_prke) > rel_tol * ratio_prke:
+            ok = False
+
+    if rank == 0:
+        print(f"PRKE_STK_PASS {1 if ok else 0}")
diff --git a/test/python/modules/linear_boltzmann_solvers/transport_transient/xs1g_delayed_crit_1p.h5 b/test/python/modules/linear_boltzmann_solvers/transport_transient/xs1g_delayed_crit_1p.h5
new file mode 100644
index 000000000..95ee272f8
Binary files /dev/null and b/test/python/modules/linear_boltzmann_solvers/transport_transient/xs1g_delayed_crit_1p.h5 differ
diff --git a/test/python/modules/linear_boltzmann_solvers/transport_transient/xs1g_delayed_crit_2p.h5 b/test/python/modules/linear_boltzmann_solvers/transport_transient/xs1g_delayed_crit_2p.h5
new file mode 100644
index 000000000..cfc9e0eeb
Binary files /dev/null and b/test/python/modules/linear_boltzmann_solvers/transport_transient/xs1g_delayed_crit_2p.h5 differ
diff --git a/test/python/modules/linear_boltzmann_solvers/transport_transient/xs1g_delayed_super_1p.h5 b/test/python/modules/linear_boltzmann_solvers/transport_transient/xs1g_delayed_super_1p.h5
new file mode 100644
index 000000000..9bd2d9b5d
Binary files /dev/null and b/test/python/modules/linear_boltzmann_solvers/transport_transient/xs1g_delayed_super_1p.h5 differ
diff --git a/test/python/modules/linear_boltzmann_solvers/transport_transient/xs1g_delayed_super_2p.h5 b/test/python/modules/linear_boltzmann_solvers/transport_transient/xs1g_delayed_super_2p.h5
new file mode 100644
index 000000000..c015e194b
Binary files /dev/null and b/test/python/modules/linear_boltzmann_solvers/transport_transient/xs1g_delayed_super_2p.h5 differ
diff --git a/test/unit/framework/materials/openmc_import_test.cc b/test/unit/framework/materials/openmc_import_test.cc
new file mode 100644
index 000000000..2e73200bb
--- /dev/null
+++ b/test/unit/framework/materials/openmc_import_test.cc
@@ -0,0 +1,191 @@
+#include "framework/materials/multi_group_xs/multi_group_xs.h"
+#include "framework/utils/hdf_utils.h"
+#include <chrono>
+#include <cstdint>
+#include <filesystem>
+#include <string>
+#include <vector>
+#include <gtest/gtest.h>
+
+using namespace opensn;
+
+namespace
+{
+
+bool
+WriteStringAttribute(hid_t id, const std::string& name, const std::string& value)
+{
+  const hid_t type = H5Tcopy(H5T_C_S1);
+  if (type < 0)
+    return false;
+  H5Tset_size(type, value.size());
+  H5Tset_strpad(type, H5T_STR_NULLTERM);
+
+  const hid_t dataspace = H5Screate(H5S_SCALAR);
+  if (dataspace < 0)
+  {
+    H5Tclose(type);
+    return false;
+  }
+
+  const hid_t attribute = H5Acreate2(id, name.c_str(), type, dataspace, H5P_DEFAULT, H5P_DEFAULT);
+  if (attribute < 0)
+  {
+    H5Sclose(dataspace);
+    H5Tclose(type);
+    return false;
+  }
+
+  const bool ok = H5Awrite(attribute, type, value.c_str()) >= 0;
+  H5Aclose(attribute);
+  H5Sclose(dataspace);
+  H5Tclose(type);
+  return ok;
+}
+
+bool
+WriteDataset2D(hid_t id,
+               const std::string& name,
+               const std::vector<double>& data,
+               const hsize_t dim0,
+               const hsize_t dim1)
+{
+  const hsize_t dims[2] = {dim0, dim1};
+  const hid_t dataspace = H5Screate_simple(2, dims, nullptr);
+  if (dataspace < 0)
+    return false;
+
+  const hid_t dataset = H5Dcreate2(
+    id, name.c_str(), H5T_NATIVE_DOUBLE, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+  if (dataset < 0)
+  {
+    H5Sclose(dataspace);
+    return false;
+  }
+
+  const bool ok =
+    H5Dwrite(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, data.data()) >= 0;
+  H5Dclose(dataset);
+  H5Sclose(dataspace);
+  return ok;
+}
+
+std::string
+MakeTempOpenMCDelayedXSFile()
+{
+  const auto now = std::chrono::steady_clock::now().time_since_epoch().count();
+  const auto path = std::filesystem::temp_directory_path() /
+                    ("opensn_openmc_delayed_test_" + std::to_string(now) + ".h5");
+
+  const H5FileHandle file(H5Fcreate(path.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT));
+  EXPECT_GE(file.Id(), 0);
+
+  EXPECT_TRUE(WriteStringAttribute(file.Id(), "filetype", "mgxs"));
+  EXPECT_TRUE(H5CreateAttribute<std::int64_t>(file.Id(), "energy_groups", 2));
+  EXPECT_TRUE(H5CreateAttribute<std::int64_t>(file.Id(), "delayed_groups", 2));
+  EXPECT_TRUE(H5WriteDataset1D<double>(file.Id(), "/group structure", {2.0e7, 1.0, 1.0e-5}));
+
+  EXPECT_TRUE(H5CreateGroup(file.Id(), "/set1"));
+  const hid_t set1 = H5Gopen2(file.Id(), "/set1", H5P_DEFAULT);
+  EXPECT_GE(set1, 0);
+  EXPECT_TRUE(H5CreateAttribute<unsigned int>(set1, "order", 0));
+  EXPECT_TRUE(H5CreateAttribute<bool>(set1, "fissionable", true));
+  EXPECT_TRUE(WriteStringAttribute(set1, "scatter_shape", "[G][G'][Order]"));
+  H5Gclose(set1);
+
+  EXPECT_TRUE(H5CreateGroup(file.Id(), "/set1/294K"));
+  const hid_t xs = H5Gopen2(file.Id(), "/set1/294K", H5P_DEFAULT);
+  EXPECT_GE(xs, 0);
+
+  EXPECT_TRUE(H5WriteDataset1D<double>(xs, "inverse-velocity", {1.0, 0.5}));
+  EXPECT_TRUE(H5WriteDataset1D<double>(xs, "total", {1.0, 2.0}));
+  EXPECT_TRUE(H5WriteDataset1D<double>(xs, "fission", {0.5, 0.25}));
+  EXPECT_TRUE(H5WriteDataset1D<double>(xs, "nu-fission", {1.4, 1.0}));
+  EXPECT_TRUE(H5WriteDataset1D<double>(xs, "prompt-nu-fission", {1.0, 0.7}));
+  EXPECT_TRUE(H5WriteDataset1D<double>(xs, "chi", {0.6, 0.4}));
+  EXPECT_TRUE(H5WriteDataset1D<double>(xs, "chi-prompt", {0.7, 0.3}));
+  EXPECT_TRUE(H5WriteDataset1D<double>(xs, "decay-rate", {0.08, 0.2}));
+  EXPECT_TRUE(WriteDataset2D(xs, "delayed-nu-fission", {0.2, 0.1, 0.2, 0.2}, 2, 2));
+  EXPECT_TRUE(WriteDataset2D(xs, "chi-delayed", {0.8, 0.2, 0.1, 0.9}, 2, 2));
+
+  H5Gclose(xs);
+
+  EXPECT_TRUE(H5CreateGroup(file.Id(), "/water"));
+  const hid_t water = H5Gopen2(file.Id(), "/water", H5P_DEFAULT);
+  EXPECT_GE(water, 0);
+  EXPECT_TRUE(H5CreateAttribute<unsigned int>(water, "order", 0));
+  EXPECT_TRUE(H5CreateAttribute<bool>(water, "fissionable", false));
+  EXPECT_TRUE(WriteStringAttribute(water, "scatter_shape", "[G][G'][Order]"));
+  H5Gclose(water);
+
+  EXPECT_TRUE(H5CreateGroup(file.Id(), "/water/294K"));
+  const hid_t water_xs = H5Gopen2(file.Id(), "/water/294K", H5P_DEFAULT);
+  EXPECT_GE(water_xs, 0);
+  EXPECT_TRUE(H5WriteDataset1D<double>(water_xs, "inverse-velocity", {1.0, 0.5}));
+  EXPECT_TRUE(H5WriteDataset1D<double>(water_xs, "total", {0.3, 0.4}));
+  H5Gclose(water_xs);
+
+  return path.string();
+}
+
+} // namespace
+
+TEST(OpenMCImportTest, ReadDelayedNeutronData)
+{
+  const std::string fname = MakeTempOpenMCDelayedXSFile();
+  const MultiGroupXS xs = MultiGroupXS::LoadFromOpenMC(fname, "set1", 294.0);
+
+  EXPECT_TRUE(xs.IsFissionable());
+  EXPECT_EQ(xs.GetNumGroups(), 2u);
+  EXPECT_EQ(xs.GetNumPrecursors(), 2u);
+
+  ASSERT_EQ(xs.GetNuSigmaF().size(), 2u);
+  EXPECT_NEAR(xs.GetNuSigmaF()[0], 1.4, 1.0e-12);
+  EXPECT_NEAR(xs.GetNuSigmaF()[1], 1.0, 1.0e-12);
+
+  ASSERT_EQ(xs.GetNuPromptSigmaF().size(), 2u);
+  EXPECT_NEAR(xs.GetNuPromptSigmaF()[0], 1.0, 1.0e-12);
+  EXPECT_NEAR(xs.GetNuPromptSigmaF()[1], 0.7, 1.0e-12);
+
+  ASSERT_EQ(xs.GetNuDelayedSigmaF().size(), 2u);
+  EXPECT_NEAR(xs.GetNuDelayedSigmaF()[0], 0.4, 1.0e-12);
+  EXPECT_NEAR(xs.GetNuDelayedSigmaF()[1], 0.3, 1.0e-12);
+
+  const auto& production = xs.GetProductionMatrix();
+  ASSERT_EQ(production.size(), 2u);
+  EXPECT_NEAR(production[0][0], 0.7, 1.0e-12);
+  EXPECT_NEAR(production[1][0], 0.3, 1.0e-12);
+  EXPECT_NEAR(production[0][1], 0.49, 1.0e-12);
+  EXPECT_NEAR(production[1][1], 0.21, 1.0e-12);
+
+  const auto& precursors = xs.GetPrecursors();
+  ASSERT_EQ(precursors.size(), 2u);
+  EXPECT_NEAR(precursors[0].decay_constant, 0.08, 1.0e-12);
+  EXPECT_NEAR(precursors[1].decay_constant, 0.2, 1.0e-12);
+  EXPECT_NEAR(precursors[0].fractional_yield, 0.3 / 0.7, 1.0e-12);
+  EXPECT_NEAR(precursors[1].fractional_yield, 0.4 / 0.7, 1.0e-12);
+  ASSERT_EQ(precursors[0].emission_spectrum.size(), 2u);
+  EXPECT_NEAR(precursors[0].emission_spectrum[0], 0.8, 1.0e-12);
+  EXPECT_NEAR(precursors[0].emission_spectrum[1], 0.2, 1.0e-12);
+  EXPECT_NEAR(precursors[1].emission_spectrum[0], 0.1, 1.0e-12);
+  EXPECT_NEAR(precursors[1].emission_spectrum[1], 0.9, 1.0e-12);
+
+  std::filesystem::remove(fname);
+}
+
+TEST(OpenMCImportTest, NonFissionableMaterialHasNoPrecursors)
+{
+  const std::string fname = MakeTempOpenMCDelayedXSFile();
+  const MultiGroupXS xs = MultiGroupXS::LoadFromOpenMC(fname, "water", 294.0);
+
+  EXPECT_FALSE(xs.IsFissionable());
+  EXPECT_EQ(xs.GetNumGroups(), 2u);
+  EXPECT_EQ(xs.GetNumPrecursors(), 0u);
+  EXPECT_TRUE(xs.GetSigmaFission().empty());
+  EXPECT_TRUE(xs.GetNuSigmaF().empty());
+  EXPECT_TRUE(xs.GetNuPromptSigmaF().empty());
+  EXPECT_TRUE(xs.GetNuDelayedSigmaF().empty());
+  EXPECT_TRUE(xs.GetPrecursors().empty());
+
+  std::filesystem::remove(fname);
+}