diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst
index c7c65e14a1d..55289864f23 100644
--- a/docs/source/pythonapi/data.rst
+++ b/docs/source/pythonapi/data.rst
@@ -63,6 +63,7 @@ Core Functions
     atomic_weight
     combine_distributions
     decay_constant
+    decay_photon_energy
     dose_coefficients
     gnd_name
     half_life
diff --git a/openmc/config.py b/openmc/config.py
index 628df7edc17..087a8c27e98 100644
--- a/openmc/config.py
+++ b/openmc/config.py
@@ -4,6 +4,7 @@
 import warnings
 
 from openmc.data import DataLibrary
+from openmc.data.decay import _DECAY_PHOTON_ENERGY
 
 __all__ = ["config"]
 
@@ -22,6 +23,10 @@ def __delitem__(self, key):
             del os.environ['OPENMC_CROSS_SECTIONS']
         elif key == 'mg_cross_sections':
             del os.environ['OPENMC_MG_CROSS_SECTIONS']
+        elif key == 'chain_file':
+            del os.environ['OPENMC_CHAIN_FILE']
+            # Reset photon source data since it relies on chain file
+            _DECAY_PHOTON_ENERGY.clear()
 
     def __setitem__(self, key, value):
         if key == 'cross_sections':
@@ -34,6 +39,8 @@ def __setitem__(self, key, value):
         elif key == 'chain_file':
             self._set_path(key, value)
             os.environ['OPENMC_CHAIN_FILE'] = str(value)
+            # Reset photon source data since it relies on chain file
+            _DECAY_PHOTON_ENERGY.clear()
         else:
             raise KeyError(f'Unrecognized config key: {key}')
 
diff --git a/openmc/data/decay.py b/openmc/data/decay.py
index 2c2505bf56f..d7ededf33dd 100644
--- a/openmc/data/decay.py
+++ b/openmc/data/decay.py
@@ -1,15 +1,19 @@
 from collections.abc import Iterable
 from io import StringIO
 from math import log
+from pathlib import Path
 import re
 from warnings import warn
+from xml.etree import ElementTree as ET
 
 import numpy as np
 from uncertainties import ufloat, UFloat
 
+import openmc
 import openmc.checkvalue as cv
+from openmc.exceptions import DataError
 from openmc.mixin import EqualityMixin
-from openmc.stats import Discrete, Tabular, combine_distributions
+from openmc.stats import Discrete, Tabular, Univariate, combine_distributions
 from .data import ATOMIC_SYMBOL, ATOMIC_NUMBER
 from .function import INTERPOLATION_SCHEME
 from .endf import Evaluation, get_head_record, get_list_record, get_tab1_record
@@ -570,3 +574,48 @@ def sources(self):
 
         self._sources = merged_sources
         return self._sources
+
+
+_DECAY_PHOTON_ENERGY = {}
+
+
+def decay_photon_energy(nuclide: str) -> Univariate:
+    """Get photon energy distribution resulting from the decay of a nuclide
+
+    This function relies on data stored in a depletion chain. Before calling it
+    for the first time, you need to ensure that a depletion chain has been
+    specified in openmc.config['chain_file'].
+
+    .. versionadded:: 0.14.0
+
+    Parameters
+    ----------
+    nuclide : str
+        Name of nuclide, e.g., 'Co58'
+
+    Returns
+    -------
+    openmc.stats.Univariate
+        Distribution of energies in [eV] of photons emitted from decay. Note
+        that the probabilities represent intensities, given as [decay/sec].
+    """
+    if not _DECAY_PHOTON_ENERGY:
+        chain_file = openmc.config.get('chain_file')
+        if chain_file is None:
+            raise DataError(
+                "A depletion chain file must be specified with "
+                "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:
+                _DECAY_PHOTON_ENERGY[nuc.name] = nuc.sources['photon']
+
+        # If the chain file contained no sources at all, warn the user
+        if not _DECAY_PHOTON_ENERGY:
+            warn(f"Chain file '{chain_file}' does not have any decay photon "
+                 "sources listed.")
+
+    return _DECAY_PHOTON_ENERGY.get(nuclide)
diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py
index 217c9cac05e..ce376a18159 100644
--- a/openmc/deplete/chain.py
+++ b/openmc/deplete/chain.py
@@ -420,6 +420,8 @@ def from_endf(cls, decay_files, fpy_files, neutron_files,
                     # Append decay mode
                     nuclide.add_decay_mode(type_, target, br)
 
+                nuclide.sources = data.sources
+
             fissionable = False
             if parent in reactions:
                 reactions_available = set(reactions[parent].keys())
diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py
index 04e0c26c589..b84b91a7ff7 100644
--- a/openmc/deplete/nuclide.py
+++ b/openmc/deplete/nuclide.py
@@ -16,6 +16,7 @@
 import numpy as np
 
 from openmc.checkvalue import check_type
+from openmc.stats import Univariate
 
 __all__ = [
     "DecayTuple", "ReactionTuple", "Nuclide", "FissionYield",
@@ -101,6 +102,9 @@ class Nuclide:
     reactions : list of openmc.deplete.ReactionTuple
         Reaction information. Each element of the list is a named tuple with
         attribute 'type', 'target', 'Q', and 'branching_ratio'.
+    sources : dict
+        Dictionary mapping particle type as string to energy distribution of
+        decay source represented as :class:`openmc.stats.Univariate`
     yield_data : FissionYieldDistribution or None
         Fission product yields at tabulated energies for this nuclide. Can be
         treated as a nested dictionary ``{energy: {product: yield}}``
@@ -120,6 +124,9 @@ def __init__(self, name=None):
         # Reaction paths
         self.reactions = []
 
+        # Decay sources
+        self.sources = {}
+
         # Neutron fission yields, if present
         self._yield_data = None
 
@@ -237,6 +244,12 @@ def from_xml(cls, element, root=None, fission_q=None):
             branching_ratio = float(decay_elem.get('branching_ratio'))
             nuc.decay_modes.append(DecayTuple(d_type, target, branching_ratio))
 
+        # Check for sources
+        for src_elem in element.iter('source'):
+            particle = src_elem.get('particle')
+            distribution = Univariate.from_xml_element(src_elem)
+            nuc.sources[particle] = distribution
+
         # Check for reaction paths
         for reaction_elem in element.iter('reaction'):
             r_type = reaction_elem.get('type')
@@ -302,6 +315,19 @@ def to_xml_element(self):
                     mode_elem.set('target', daughter)
                 mode_elem.set('branching_ratio', str(br))
 
+        # Write decay sources
+        if self.sources:
+            for particle, source in self.sources.items():
+                # TODO: Ugly hack to deal with the fact that
+                # 'source.to_xml_element' will return an xml.etree object
+                # whereas here lxml is being used preferentially. We should just
+                # switch to use lxml everywhere
+                import xml.etree.ElementTree as etree
+                src_elem_xmletree = source.to_xml_element('source')
+                src_elem = ET.fromstring(etree.tostring(src_elem_xmletree))
+                src_elem.set('particle', particle)
+                elem.append(src_elem)
+
         elem.set('reactions', str(len(self.reactions)))
         for rx, daughter, Q, br in self.reactions:
             rx_elem = ET.SubElement(elem, 'reaction')
diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py
index ca1516803c2..fd9e665248d 100644
--- a/openmc/deplete/openmc_operator.py
+++ b/openmc/deplete/openmc_operator.py
@@ -334,7 +334,7 @@ def _set_number_from_results(self, mat, prev_res):
         mat_id = str(mat.id)
 
         # Get nuclide lists from geometry and depletion results
-        depl_nuc = prev_res[-1].nuc_to_ind
+        depl_nuc = prev_res[-1].index_nuc
         geom_nuc_densities = mat.get_nuclide_atom_densities()
 
         # Merge lists of nuclides, with the same order for every calculation
diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py
index b2c8e77332e..d9aa52b79e4 100644
--- a/openmc/deplete/results.py
+++ b/openmc/deplete/results.py
@@ -435,7 +435,7 @@ def export_to_materials(
         # results, and save them to the new depleted xml file.
         for mat in mat_file:
             mat_id = str(mat.id)
-            if mat_id in result.mat_to_ind:
+            if mat_id in result.index_mat:
                 mat.volume = result.volume[mat_id]
 
                 # Change density of all nuclides in material to atom/b-cm
@@ -447,7 +447,7 @@ def export_to_materials(
 
                 # For nuclides in chain that have cross sections, replace
                 # density in original material with new density from results
-                for nuc in result.nuc_to_ind:
+                for nuc in result.index_nuc:
                     if nuc not in available_cross_sections:
                         continue
                     atoms = result[0, mat_id, nuc]
diff --git a/openmc/deplete/stepresult.py b/openmc/deplete/stepresult.py
index db5260c8080..e462f175420 100644
--- a/openmc/deplete/stepresult.py
+++ b/openmc/deplete/stepresult.py
@@ -6,6 +6,7 @@
 
 from collections import OrderedDict
 import copy
+import warnings
 
 import h5py
 import numpy as np
@@ -42,9 +43,9 @@ class StepResult:
         The reaction rates for each substep.
     volume : OrderedDict of str to float
         Dictionary mapping mat id to volume.
-    mat_to_ind : OrderedDict of str to int
+    index_mat : OrderedDict of str to int
         A dictionary mapping mat ID as string to index.
-    nuc_to_ind : OrderedDict of str to int
+    index_nuc : OrderedDict of str to int
         A dictionary mapping nuclide name as string to index.
     mat_to_hdf5_ind : OrderedDict of str to int
         A dictionary mapping mat ID as string to global index.
@@ -67,8 +68,8 @@ def __init__(self):
         self.volume = None
         self.proc_time = None
 
-        self.mat_to_ind = None
-        self.nuc_to_ind = None
+        self.index_mat = None
+        self.index_nuc = None
         self.mat_to_hdf5_ind = None
 
         self.data = None
@@ -98,9 +99,9 @@ def __getitem__(self, pos):
         if isinstance(mat, openmc.Material):
             mat = str(mat.id)
         if isinstance(mat, str):
-            mat = self.mat_to_ind[mat]
+            mat = self.index_mat[mat]
         if isinstance(nuc, str):
-            nuc = self.nuc_to_ind[nuc]
+            nuc = self.index_nuc[nuc]
 
         return self.data[stage, mat, nuc]
 
@@ -120,19 +121,19 @@ def __setitem__(self, pos, val):
         """
         stage, mat, nuc = pos
         if isinstance(mat, str):
-            mat = self.mat_to_ind[mat]
+            mat = self.index_mat[mat]
         if isinstance(nuc, str):
-            nuc = self.nuc_to_ind[nuc]
+            nuc = self.index_nuc[nuc]
 
         self.data[stage, mat, nuc] = val
 
     @property
     def n_mat(self):
-        return len(self.mat_to_ind)
+        return len(self.index_mat)
 
     @property
     def n_nuc(self):
-        return len(self.nuc_to_ind)
+        return len(self.index_nuc)
 
     @property
     def n_hdf5_mats(self):
@@ -160,8 +161,8 @@ def allocate(self, volume, nuc_list, burn_list, full_burn_list, stages):
 
         """
         self.volume = copy.deepcopy(volume)
-        self.nuc_to_ind = {nuc: i for i, nuc in enumerate(nuc_list)}
-        self.mat_to_ind = {mat: i for i, mat in enumerate(burn_list)}
+        self.index_nuc = {nuc: i for i, nuc in enumerate(nuc_list)}
+        self.index_mat = {mat: i for i, mat in enumerate(burn_list)}
         self.mat_to_hdf5_ind = {mat: i for i, mat in enumerate(full_burn_list)}
 
         # Create storage array
@@ -186,10 +187,10 @@ def distribute(self, local_materials, ranges):
         """
         new = StepResult()
         new.volume = {lm: self.volume[lm] for lm in local_materials}
-        new.mat_to_ind = {mat: idx for (idx, mat) in enumerate(local_materials)}
+        new.index_mat = {mat: idx for (idx, mat) in enumerate(local_materials)}
 
         # Direct transfer
-        direct_attrs = ("time", "k", "source_rate", "nuc_to_ind",
+        direct_attrs = ("time", "k", "source_rate", "index_nuc",
                         "mat_to_hdf5_ind", "proc_time")
         for attr in direct_attrs:
             setattr(new, attr, getattr(self, attr))
@@ -198,6 +199,34 @@ def distribute(self, local_materials, ranges):
         new.rates = [r[ranges] for r in self.rates]
         return new
 
+    def get_material(self, mat_id):
+        """Return material object for given depleted composition
+
+        .. versionadded:: 0.14.0
+
+        Parameters
+        ----------
+        mat_id : str
+            Material ID as a string
+
+        Returns
+        -------
+        openmc.Material
+            Equivalent material
+        """
+        with warnings.catch_warnings():
+            warnings.simplefilter('ignore', openmc.IDWarning)
+            material = openmc.Material(material_id=int(mat_id))
+        vol = self.volume[mat_id]
+        for nuc, _ in sorted(self.index_nuc.items(), key=lambda x: x[1]):
+            atoms = self[0, mat_id, nuc]
+            if atoms < 0.0:
+                continue
+            atom_per_bcm = atoms / vol * 1e-24
+            material.add_nuclide(nuc, atom_per_bcm)
+        material.volume = vol
+        return material
+
     def export_to_hdf5(self, filename, step):
         """Export results to an HDF5 file
 
@@ -239,11 +268,11 @@ def _write_hdf5_metadata(self, handle):
         """
         # Create and save the 5 dictionaries:
         # quantities
-        #   self.mat_to_ind -> self.volume (TODO: support for changing volumes)
-        #   self.nuc_to_ind
+        #   self.index_mat -> self.volume (TODO: support for changing volumes)
+        #   self.index_nuc
         # reactions
-        #   self.rates[0].nuc_to_ind (can be different from above, above is superset)
-        #   self.rates[0].react_to_ind
+        #   self.rates[0].index_nuc (can be different from above, above is superset)
+        #   self.rates[0].index_rx
         # these are shared by every step of the simulation, and should be deduplicated.
 
         # Store concentration mat and nuclide dictionaries (along with volumes)
@@ -252,7 +281,7 @@ def _write_hdf5_metadata(self, handle):
         handle.attrs['filetype'] = np.string_('depletion results')
 
         mat_list = sorted(self.mat_to_hdf5_ind, key=int)
-        nuc_list = sorted(self.nuc_to_ind)
+        nuc_list = sorted(self.index_nuc)
         rxn_list = sorted(self.rates[0].index_rx)
 
         n_mats = self.n_hdf5_mats
@@ -272,7 +301,7 @@ def _write_hdf5_metadata(self, handle):
 
         for nuc in nuc_list:
             nuc_single_group = nuc_group.create_group(nuc)
-            nuc_single_group.attrs["atom number index"] = self.nuc_to_ind[nuc]
+            nuc_single_group.attrs["atom number index"] = self.index_nuc[nuc]
             if nuc in self.rates[0].index_nuc:
                 nuc_single_group.attrs["reaction rate index"] = self.rates[0].index_nuc[nuc]
 
@@ -367,13 +396,13 @@ def _to_hdf5(self, handle, index, parallel=False):
             proc_time_dset.resize(proc_shape)
 
         # If nothing to write, just return
-        if len(self.mat_to_ind) == 0:
+        if len(self.index_mat) == 0:
             return
 
         # Add data
         # Note, for the last step, self.n_stages = 1, even if n_stages != 1.
         n_stages = self.n_stages
-        inds = [self.mat_to_hdf5_ind[mat] for mat in self.mat_to_ind]
+        inds = [self.mat_to_hdf5_ind[mat] for mat in self.index_mat]
         low = min(inds)
         high = max(inds)
         for i in range(n_stages):
@@ -427,8 +456,8 @@ def from_hdf5(cls, handle, step):
 
         # Reconstruct dictionaries
         results.volume = OrderedDict()
-        results.mat_to_ind = OrderedDict()
-        results.nuc_to_ind = OrderedDict()
+        results.index_mat = OrderedDict()
+        results.index_nuc = OrderedDict()
         rxn_nuc_to_ind = OrderedDict()
         rxn_to_ind = OrderedDict()
 
@@ -437,11 +466,11 @@ def from_hdf5(cls, handle, step):
             ind = mat_handle.attrs["index"]
 
             results.volume[mat] = vol
-            results.mat_to_ind[mat] = ind
+            results.index_mat[mat] = ind
 
         for nuc, nuc_handle in handle["/nuclides"].items():
             ind_atom = nuc_handle.attrs["atom number index"]
-            results.nuc_to_ind[nuc] = ind_atom
+            results.index_nuc[nuc] = ind_atom
 
             if "reaction rate index" in nuc_handle.attrs:
                 rxn_nuc_to_ind[nuc] = nuc_handle.attrs["reaction rate index"]
@@ -452,7 +481,7 @@ def from_hdf5(cls, handle, step):
         results.rates = []
         # Reconstruct reactions
         for i in range(results.n_stages):
-            rate = ReactionRates(results.mat_to_ind, rxn_nuc_to_ind, rxn_to_ind, True)
+            rate = ReactionRates(results.index_mat, rxn_nuc_to_ind, rxn_to_ind, True)
 
             rate[:] = handle["/reaction rates"][step, i, :, :, :]
             results.rates.append(rate)
diff --git a/openmc/material.py b/openmc/material.py
index 1b2b46a73a1..16b5a3048aa 100644
--- a/openmc/material.py
+++ b/openmc/material.py
@@ -18,6 +18,7 @@
 from ._xml import clean_indentation, reorder_attributes
 from .mixin import IDManagerMixin
 from openmc.checkvalue import PathLike
+from openmc.stats import Univariate
 
 
 # Units for density supported by OpenMC
@@ -34,10 +35,10 @@ class Material(IDManagerMixin):
     To create a material, one should create an instance of this class, add
     nuclides or elements with :meth:`Material.add_nuclide` or
     :meth:`Material.add_element`, respectively, and set the total material
-    density with :meth:`Material.set_density()`. Alternatively, you can
-    use :meth:`Material.add_components()` to pass a dictionary
-    containing all the component information. The material can then be
-    assigned to a cell using the :attr:`Cell.fill` attribute.
+    density with :meth:`Material.set_density()`. Alternatively, you can use
+    :meth:`Material.add_components()` to pass a dictionary containing all the
+    component information. The material can then be assigned to a cell using the
+    :attr:`Cell.fill` attribute.
 
     Parameters
     ----------
@@ -91,6 +92,13 @@ class Material(IDManagerMixin):
     fissionable_mass : float
         Mass of fissionable nuclides in the material in [g]. Requires that the
         :attr:`volume` attribute is set.
+    decay_photon_energy : openmc.stats.Univariate
+        Energy distribution of photons emitted from decay of unstable nuclides
+        within the material. The integral of this distribution is the total
+        intensity of the photon source in [decay/sec].
+
+        .. versionadded:: 0.14.0
+
     """
 
     next_id = 1
@@ -255,6 +263,18 @@ def fissionable_mass(self):
                            / openmc.data.AVOGADRO
         return density*self.volume
 
+    @property
+    def decay_photon_energy(self) -> Univariate:
+        atoms = self.get_nuclide_atoms()
+        dists = []
+        probs = []
+        for nuc, num_atoms in atoms.items():
+            source_per_atom = openmc.data.decay_photon_energy(nuc)
+            if source_per_atom is not None:
+                dists.append(source_per_atom)
+                probs.append(num_atoms)
+        return openmc.data.combine_distributions(dists, probs)
+
     @classmethod
     def from_hdf5(cls, group: h5py.Group):
         """Create material from HDF5 group
@@ -820,7 +840,7 @@ def get_nuclides(self, element: Optional[str] = None):
             for nuclide in self._nuclides:
                 if nuclide.name not in matching_nuclides:
                     matching_nuclides.append(nuclide.name)
-                        
+
         return matching_nuclides
 
     def get_nuclide_densities(self):
diff --git a/tests/chain_simple.xml b/tests/chain_simple.xml
index c2e50a370f3..28c5739c9fc 100644
--- a/tests/chain_simple.xml
+++ b/tests/chain_simple.xml
@@ -3,10 +3,16 @@
   <nuclide name="I135" decay_modes="1" reactions="1" half_life="2.36520E+04">
     <decay type="beta" target="Xe135" branching_ratio="1.0" />
     <reaction type="(n,gamma)" Q="0.0" target="Xe136" /> <!-- Not precisely true, but whatever -->
+    <source type="discrete" particle="photon">
+      <parameters>3696.125 4095.822 4477.27 5097.122 29452.1 29781.3 33566.5 33629.4 33865.1 33878.5 34395.3 34408.0 34486.2 34488.2 112780.0 113150.0 162650.0 165740.0 184490.0 197190.0 220502.0 229720.0 247500.0 254740.0 264260.0 288451.0 290270.0 304910.0 305830.0 326000.0 333600.0 342520.0 361850.0 403030.0 414830.0 417633.0 429930.0 433741.0 451630.0 530800.0 546557.0 575970.0 588280.0 616900.0 649850.0 656090.0 679220.0 684600.0 690130.0 707920.0 785480.0 795500.0 797710.0 807200.0 836804.0 960290.0 961430.0 971960.0 972620.0 995090.0 1038760.0 1096860.0 1101580.0 1124000.0 1131511.0 1151510.0 1159900.0 1169040.0 1225600.0 1240470.0 1254800.0 1260409.0 1315770.0 1334800.0 1343660.0 1367890.0 1441800.0 1448350.0 1457560.0 1502790.0 1521990.0 1543700.0 1566410.0 1678027.0 1706459.0 1791196.0 1830690.0 1845300.0 1927300.0 1948490.0 2045880.0 2112400.0 2151500.0 2189400.0 2255457.0 2408650.0 2466070.0 2477100.0 9.714352819815078e-10 7.460941651551526e-09 6.047882056201745e-09 8.510389107205747e-10 3.7979729633684727e-08 7.033747061198817e-08 6.602815946851458e-09 1.2800909198245071e-08 6.513060244589454e-11 8.894667887850525e-11 1.3843783302275766e-09 2.700005498135763e-09 1.2141115256573815e-11 1.654893875618862e-11 3.7007705885806645e-09 2.0186021392258173e-09 2.859686363903241e-09 9.16781804898392e-09 6.896890642354875e-09 9.588360161322631e-09 5.130613770532286e-07 7.06510748729036e-08 8.410842246774237e-09 6.7286737974193905e-09 5.3829390379355124e-08 9.083709626516178e-07 8.915492781580693e-08 9.251926471451662e-09 2.783988783682273e-08 6.728673797419391e-10 1.093409492080651e-08 2.5232526740322717e-10 5.467047460403255e-08 6.812782219887132e-08 8.83138435911295e-08 1.0345335963532313e-06 8.915492781580693e-08 1.623292553627428e-07 9.251926471451663e-08 9.251926471451662e-09 2.0942997194467853e-06 3.784879011048407e-08 1.513951604419363e-08 1.093409492080651e-08 1.337323917237104e-07 2.186818984161302e-08 1.5980600268871052e-08 6.7286737974193905e-09 3.784879011048407e-08 1.9344937167580748e-07 4.457746390790346e-08 6.7286737974193905e-09 5.0465053480645434e-08 1.3457347594838781e-08 1.9597262434983976e-06 1.0093010696129087e-08 4.289529545854862e-08 2.607361096500014e-07 3.53255374364518e-07 4.541854813258089e-08 2.329803302356464e-06 2.607361096500014e-08 4.7100716581935733e-07 1.059766123093554e-06 6.6193328482113254e-06 4.205421123387119e-10 3.027903208838726e-08 2.565306885266143e-07 1.2616263370161358e-08 2.649415307733885e-07 3.3643368987096953e-09 8.410842246774237e-06 1.934493716758075e-08 9.251926471451662e-09 2.2709274066290446e-08 1.7830985563161385e-07 5.046505348064544e-09 9.251926471451663e-08 2.5400743585258203e-06 3.154065842540339e-07 1.093409492080651e-08 7.569758022096815e-09 3.784879011048407e-07 2.8008104681758214e-06 1.2027504412887161e-06 2.26251656438227e-06 1.6989901338483962e-07 1.6821684493548476e-09 8.663167514177466e-08 1.8503852942903323e-08 2.5568960430193683e-07 2.0186021392258175e-08 6.560456952483906e-09 3.784879011048407e-09 1.7999202408096872e-07 2.800810468175822e-07 2.1027105616935597e-08 4.205421123387119e-10</parameters>
+    </source>
   </nuclide>
   <nuclide name="Xe135" decay_modes="1" reactions="1" half_life="3.29040E+04">
     <decay type=" beta" target="Cs135" branching_ratio="1.0" />
     <reaction type="(n,gamma)" Q="0.0" target="Xe136" />
+    <source type="tabular" interpolation="histogram" particle="photon">
+      <parameters>0.0 10000.0 20000.0 50000.0 100000.0 200000.0 300000.0 400000.0 600000.0 800000.0 1000000.0 1220000.0 1.1612176249914943e-11 0.0 3.1976524142990123e-11 0.0 6.418466676129901e-13 1.894551923065874e-10 5.099949011192198e-13 3.170644340540729e-13 2.756798470867507e-12 6.67627508515641e-14 3.711746246444946e-15 0.0</parameters>
+    </source>
   </nuclide>
   <nuclide name="Xe136" decay_modes="0" reactions="0" />
   <nuclide name="Cs135" decay_modes="0" reactions="0" />
@@ -28,6 +34,9 @@
   </nuclide>
   <nuclide name="U235" decay_modes="0" reactions="1">
     <reaction type="fission" Q="193410000."/>
+    <source type="discrete" particle="photon">
+      <parameters>3065.349 12960.11 13197.49 16125.43 19185.05 19590.0 31600.0 34700.0 41400.0 41960.0 51220.0 54100.0 54250.0 64350.0 72700.0 75020.0 76198.0 90330.0 93795.0 96090.0 105278.0 106074.0 106608.0 106771.0 108948.0 109154.0 109160.0 109395.0 109433.0 115450.0 120350.0 136550.0 140760.0 142400.0 143760.0 150930.0 163330.0 173300.0 182610.0 185715.0 194940.0 198900.0 202110.0 205311.0 215280.0 221380.0 228780.0 233500.0 240870.0 246840.0 266450.0 275129.0 275430.0 281420.0 282920.0 289560.0 291650.0 301700.0 317100.0 343500.0 345900.0 356030.0 387820.0 410290.0 428710.0 448400.0 1.4211820389290126e-18 3.695705148646686e-18 4.752472005970374e-19 4.0825252294602915e-18 8.998127221991115e-19 1.9035285506819052e-21 5.305446177665699e-21 1.1547147563154756e-20 9.362552078233586e-21 1.86835843588509e-20 1.0610892355331398e-20 2.736290732002608e-22 4.776811520523089e-21 3.977552295559136e-21 3.432935762018982e-20 1.872510415646717e-20 2.4966805541956234e-21 1.0265454754703584e-18 1.7575878976520235e-18 2.839974130397521e-20 2.1057468800839102e-19 4.1342252420362975e-19 7.098565272539688e-21 8.20050332279016e-21 5.276915360632628e-20 1.0602918581811435e-19 4.806110066826575e-19 2.0521431485853304e-21 2.375669880669523e-21 9.362552078233586e-21 8.267152210184412e-21 3.745020831293435e-21 6.865871524037964e-20 1.5604253463722645e-21 3.420452359248004e-18 2.4966805541956232e-20 1.5853921519142206e-18 1.8725104156467174e-21 1.0610892355331397e-19 1.7851265962498703e-17 1.966135936429053e-19 1.3107572909527022e-20 3.370518748164091e-19 1.5635461970650089e-18 9.050467008959134e-21 3.745020831293434e-20 2.18459548492117e-21 9.050467008959134e-21 2.3406380195583967e-20 1.6540508671546e-20 1.8725104156467174e-21 1.622842360227155e-20 2.18459548492117e-21 1.8725104156467174e-21 1.8725104156467174e-21 2.18459548492117e-21 1.2483402770978116e-20 1.5604253463722645e-21 3.120850692744529e-22 9.362552078233587e-22 1.2483402770978116e-20 1.5604253463722645e-21 1.2483402770978116e-20 9.362552078233587e-22 3.120850692744529e-22 3.120850692744529e-22</parameters>
+    </source>
     <neutron_fission_yields>
       <energies>2.53000e-02</energies>
       <fission_yields energy="2.53000e-02">
@@ -38,6 +47,9 @@
   </nuclide>
   <nuclide name="U238" decay_modes="0" reactions="1">
     <reaction type="fission" Q="197790000."/>
+    <source type="discrete" particle="photon">
+      <parameters>3061.32 12959.8 13440.07 16150.05 19148.83 49550.0 90330.0 93795.0 105278.0 106074.0 106608.0 106771.0 108948.0 109154.0 109395.0 109433.0 113500.0 5.3130501476983077e-20 1.4936931167066328e-19 1.943509007980312e-20 1.8224649880365157e-19 4.0452400597373603e-20 3.146222282132249e-21 3.300026341588747e-23 5.444173877278485e-23 6.3486292118621375e-24 1.2400176384733938e-23 2.124537722121886e-25 2.4542156071495764e-25 1.5792541400719878e-24 3.1731991718126067e-24 6.141568457919309e-26 7.109808533300877e-26 5.014291762148272e-22</parameters>
+    </source>
     <neutron_fission_yields>
       <energies>2.53000e-02</energies>
       <fission_yields energy="2.53000e-02">
diff --git a/tests/regression_tests/deplete_no_transport/test.py b/tests/regression_tests/deplete_no_transport/test.py
index f5ee4bf1dbb..384653ff7e6 100644
--- a/tests/regression_tests/deplete_no_transport/test.py
+++ b/tests/regression_tests/deplete_no_transport/test.py
@@ -161,23 +161,23 @@ def _create_operator(from_nuclides,
     return op
 
 def _assert_same_mats(res_ref, res_test):
-    for mat in res_ref[0].mat_to_ind:
-            assert mat in res_test[0].mat_to_ind, \
+    for mat in res_ref[0].index_mat:
+            assert mat in res_test[0].index_mat, \
                 "Material {} not in new results.".format(mat)
-    for nuc in res_ref[0].nuc_to_ind:
-        assert nuc in res_test[0].nuc_to_ind, \
+    for nuc in res_ref[0].index_nuc:
+        assert nuc in res_test[0].index_nuc, \
             "Nuclide {} not in new results.".format(nuc)
 
-    for mat in res_test[0].mat_to_ind:
-        assert mat in res_ref[0].mat_to_ind, \
+    for mat in res_test[0].index_mat:
+        assert mat in res_ref[0].index_mat, \
             "Material {} not in old results.".format(mat)
-    for nuc in res_test[0].nuc_to_ind:
-        assert nuc in res_ref[0].nuc_to_ind, \
+    for nuc in res_test[0].index_nuc:
+        assert nuc in res_ref[0].index_nuc, \
             "Nuclide {} not in old results.".format(nuc)
 
 def _assert_atoms_equal(res_ref, res_test, tol):
-    for mat in res_test[0].mat_to_ind:
-        for nuc in res_test[0].nuc_to_ind:
+    for mat in res_test[0].index_mat:
+        for nuc in res_test[0].index_nuc:
             _, y_test = res_test.get_atoms(mat, nuc)
             _, y_old = res_ref.get_atoms(mat, nuc)
 
diff --git a/tests/regression_tests/deplete_with_transport/test.py b/tests/regression_tests/deplete_with_transport/test.py
index 58065035f3b..9b37e43510e 100644
--- a/tests/regression_tests/deplete_with_transport/test.py
+++ b/tests/regression_tests/deplete_with_transport/test.py
@@ -81,23 +81,23 @@ def test_full(run_in_tmpdir, problem, multiproc):
     res_ref = openmc.deplete.Results(path_reference)
 
     # Assert same mats
-    for mat in res_ref[0].mat_to_ind:
-        assert mat in res_test[0].mat_to_ind, \
+    for mat in res_ref[0].index_mat:
+        assert mat in res_test[0].index_mat, \
             "Material {} not in new results.".format(mat)
-    for nuc in res_ref[0].nuc_to_ind:
-        assert nuc in res_test[0].nuc_to_ind, \
+    for nuc in res_ref[0].index_nuc:
+        assert nuc in res_test[0].index_nuc, \
             "Nuclide {} not in new results.".format(nuc)
 
-    for mat in res_test[0].mat_to_ind:
-        assert mat in res_ref[0].mat_to_ind, \
+    for mat in res_test[0].index_mat:
+        assert mat in res_ref[0].index_mat, \
             "Material {} not in old results.".format(mat)
-    for nuc in res_test[0].nuc_to_ind:
-        assert nuc in res_ref[0].nuc_to_ind, \
+    for nuc in res_test[0].index_nuc:
+        assert nuc in res_ref[0].index_nuc, \
             "Nuclide {} not in old results.".format(nuc)
 
     tol = 1.0e-6
-    for mat in res_test[0].mat_to_ind:
-        for nuc in res_test[0].nuc_to_ind:
+    for mat in res_test[0].index_mat:
+        for nuc in res_test[0].index_nuc:
             _, y_test = res_test.get_atoms(mat, nuc)
             _, y_old = res_ref.get_atoms(mat, nuc)
 
diff --git a/tests/unit_tests/test_data_decay.py b/tests/unit_tests/test_data_decay.py
index f0bda1bd978..8900e7eace5 100644
--- a/tests/unit_tests/test_data_decay.py
+++ b/tests/unit_tests/test_data_decay.py
@@ -1,13 +1,14 @@
 #!/usr/bin/env python
 
-from collections.abc import Mapping
 import os
 from math import log
+from pathlib import Path
 
 import numpy as np
 import pytest
 from uncertainties import ufloat
 import openmc.data
+from openmc.exceptions import DataError
 
 
 def ufloat_close(a, b):
@@ -126,3 +127,25 @@ def test_sources(ba137m, nb90):
     # Nb90 decays by β+ and should emit positrons, electrons, and photons
     sources = nb90.sources
     assert len(set(sources.keys()) ^ {'positron', 'electron', 'photon'}) == 0
+
+
+def test_decay_photon_energy():
+    # If chain file is not set, we should get a data error
+    if 'chain_file' in openmc.config:
+        del openmc.config['chain_file']
+    with pytest.raises(DataError):
+        openmc.data.decay_photon_energy('I135')
+
+    # Set chain file to simple chain
+    openmc.config['chain_file'] = Path(__file__).parents[1] / "chain_simple.xml"
+
+    # Check strength of I135 source and presence of specific spectral line
+    src = openmc.data.decay_photon_energy('I135')
+    assert isinstance(src, openmc.stats.Discrete)
+    assert src.integral() == pytest.approx(3.920996223799345e-05)
+    assert 1260409. in src.x
+
+    # Check Xe135 source, which should be tabular
+    src = openmc.data.decay_photon_energy('Xe135')
+    assert isinstance(src, openmc.stats.Tabular)
+    assert src.integral() == pytest.approx(2.076506258964966e-05)
diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py
index 4bfdaa63e58..32c3e28098c 100644
--- a/tests/unit_tests/test_deplete_chain.py
+++ b/tests/unit_tests/test_deplete_chain.py
@@ -9,6 +9,7 @@
 import numpy as np
 from openmc.mpi import comm
 from openmc.deplete import Chain, reaction_rates, nuclide, cram, pool
+from openmc.stats import Discrete
 import pytest
 
 from tests import cdtemp
@@ -476,6 +477,15 @@ def gnd_simple_chain():
     return Chain.from_xml(chainfile)
 
 
+def test_chain_sources(gnd_simple_chain):
+    i135 = gnd_simple_chain['I135']
+    assert isinstance(i135.sources, dict)
+    assert list(i135.sources.keys()) == ['photon']
+    photon_src = i135.sources['photon']
+    assert isinstance(photon_src, Discrete)
+    assert photon_src.integral() == pytest.approx(3.920996223799345e-05)
+
+
 def test_reduce(gnd_simple_chain, endf_chain):
     ref_U5 = gnd_simple_chain["U235"]
     ref_iodine = gnd_simple_chain["I135"]
diff --git a/tests/unit_tests/test_deplete_resultslist.py b/tests/unit_tests/test_deplete_resultslist.py
index 63073cf57a3..e13944e0358 100644
--- a/tests/unit_tests/test_deplete_resultslist.py
+++ b/tests/unit_tests/test_deplete_resultslist.py
@@ -133,3 +133,17 @@ def test_get_steps(unit):
     actual = results.get_step_where(
         times[-1] * 100, time_units=unit, atol=inf, rtol=inf)
     assert actual == times.size - 1
+
+
+def test_stepresult_get_material(res):
+    # Get material at first timestep
+    step_result = res[0]
+    mat1 = step_result.get_material("1")
+    assert mat1.id == 1
+    assert mat1.volume == step_result.volume["1"]
+
+    # Spot check number densities
+    densities = mat1.get_nuclide_atom_densities()
+    assert densities['Xe135'] == pytest.approx(1e-14)
+    assert densities['I135'] == pytest.approx(1e-21)
+    assert densities['U234'] == pytest.approx(1.00506e-05)
diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py
index 7b8a60629e6..58df246ddd7 100644
--- a/tests/unit_tests/test_material.py
+++ b/tests/unit_tests/test_material.py
@@ -1,8 +1,10 @@
 from collections import defaultdict
+from pathlib import Path
 
 import pytest
 
 import openmc
+from openmc.data import decay_photon_energy
 import openmc.examples
 import openmc.model
 import openmc.stats
@@ -541,3 +543,33 @@ def test_get_activity():
     # volume is required to calculate total activity
     m4.volume = 10.
     assert pytest.approx(m4.get_activity(units='Bq')) == 355978108155965.94*3/2*10 # [Bq]
+
+
+def test_decay_photon_energy():
+    # Set chain file for testing
+    openmc.config['chain_file'] = Path(__file__).parents[1] / 'chain_simple.xml'
+
+    # Material representing single atom of I135 and Cs135
+    m = openmc.Material()
+    m.add_nuclide('I135', 1.0e-24)
+    m.add_nuclide('Cs135', 1.0e-24)
+    m.volume = 1.0
+
+    # Get decay photon source and make sure it's the right type
+    src = m.decay_photon_energy
+    assert isinstance(src, openmc.stats.Discrete)
+
+    # If we add Xe135 (which has a tabular distribution), the photon source
+    # should be a mixture distribution
+    m.add_nuclide('Xe135', 1.0e-24)
+    src = m.decay_photon_energy
+    assert isinstance(src, openmc.stats.Mixture)
+
+    # With a single atom of each, the intensity of the photon source should be
+    # equal to the sum of the intensities for each nuclide
+    def intensity(src):
+        return src.integral() if src is not None else 0.0
+
+    assert src.integral() == pytest.approx(sum(
+        intensity(decay_photon_energy(nuc)) for nuc in m.get_nuclides()
+    ))