diff --git a/openmc/data/decay.py b/openmc/data/decay.py
index f358dadf74..96c3f562af 100644
--- a/openmc/data/decay.py
+++ b/openmc/data/decay.py
@@ -316,6 +316,12 @@ class Decay(EqualityMixin):
         'excited_state', 'mass', 'stable', 'spin', and 'parity'.
     spectra : dict
         Resulting radiation spectra for each radiation type.
+    sources : dict
+        Radioactive decay source distributions represented as a dictionary
+        mapping particle types (e.g., 'photon') to instances of
+        :class:`openmc.stats.Univariate`.
+
+        .. versionadded:: 0.13.1
 
     """
     def __init__(self, ev_or_filename):
@@ -498,16 +504,14 @@ class Decay(EqualityMixin):
         """
         return cls(ev_or_filename)
 
-    def get_sources(self):
-        """Get radioactive decay source distributions
+    @property
+    def sources(self):
+        """Radioactive decay source distributions"""
+        # If property has been computed already, return it
+        # TODO: Replace with functools.cached_property when support is Python 3.9+
+        if self._sources is not None:
+            return self._sources
 
-        Returns
-        -------
-        sources : dict
-            Dictionary mapping particle types (e.g., 'photon') to instances of
-            :class:`openmc.stats.Univariate`
-
-        """
         sources = {}
         name = self.nuclide['name']
         decay_constant = self.decay_constant.n
@@ -563,6 +567,5 @@ class Decay(EqualityMixin):
             merged_sources[particle_type] = combine_distributions(
                 dist_list, [1.0]*len(dist_list))
 
-        return merged_sources
-
-
+        self._sources = merged_sources
+        return self._sources
diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py
index e11c9b50ee..66c6e5f064 100644
--- a/openmc/stats/univariate.py
+++ b/openmc/stats/univariate.py
@@ -224,6 +224,8 @@ class Discrete(Univariate):
     def integral(self):
         """Return integral of distribution
 
+        .. versionadded:: 0.13.1
+
         Returns
         -------
         float
@@ -1041,6 +1043,8 @@ class Tabular(Univariate):
     def integral(self):
         """Return integral of distribution
 
+        .. versionadded: 0.13.1
+
         Returns
         -------
         float
@@ -1230,6 +1234,8 @@ class Mixture(Univariate):
     def integral(self):
         """Return integral of the distribution
 
+        .. versionadded:: 0.13.1
+
         Returns
         -------
         float
@@ -1250,6 +1256,8 @@ def combine_distributions(dists, probs):
     distribution and the remainder of the distributions are put into a
     :class:`~openmc.stats.Mixture` distribution.
 
+    .. versionadded:: 0.13.1
+
     Parameters
     ----------
     dists : iterable of openmc.stats.Univariate
diff --git a/tests/unit_tests/test_data_decay.py b/tests/unit_tests/test_data_decay.py
index 06b8e6bedd..f0bda1bd97 100644
--- a/tests/unit_tests/test_data_decay.py
+++ b/tests/unit_tests/test_data_decay.py
@@ -23,6 +23,14 @@ def nb90():
     return openmc.data.Decay.from_endf(filename)
 
 
+@pytest.fixture(scope='module')
+def ba137m():
+    """Ba137_m1 decay data."""
+    endf_data = os.environ['OPENMC_ENDF_DATA']
+    filename = os.path.join(endf_data, 'decay', 'dec-056_Ba_137m1.endf')
+    return openmc.data.Decay.from_endf(filename)
+
+
 @pytest.fixture(scope='module')
 def u235_yields():
     """U235 fission product yield data."""
@@ -48,6 +56,7 @@ def test_nb90_halflife(nb90):
     ufloat_close(nb90.decay_constant, log(2.)/nb90.half_life)
     ufloat_close(nb90.decay_energy, ufloat(2265527.5, 25159.400474401213))
 
+
 def test_nb90_nuclide(nb90):
     assert nb90.nuclide['atomic_number'] == 41
     assert nb90.nuclide['mass_number'] == 90
@@ -91,3 +100,29 @@ def test_fpy(u235_yields):
     assert len(u235_yields.independent) == 3
     thermal = u235_yields.independent[0]
     ufloat_close(thermal['I135'], ufloat(0.0292737, 0.000819663))
+
+
+def test_sources(ba137m, nb90):
+    # Running .sources twice should give same objects
+    sources = ba137m.sources
+    sources2 = ba137m.sources
+    for key in sources:
+        assert sources[key] is sources2[key]
+
+    # Each source should be a univariate distribution
+    for dist in sources.values():
+        assert isinstance(dist, openmc.stats.Univariate)
+
+    # Check for presence of 662 keV gamma ray in decay of Ba137m
+    gamma_source = ba137m.sources['photon']
+    assert isinstance(gamma_source, openmc.stats.Discrete)
+    b = np.isclose(gamma_source.x, 661657.)
+    assert np.count_nonzero(b) == 1
+
+    # Check value of decay/s/atom
+    idx = np.flatnonzero(b)[0]
+    assert gamma_source.p[idx] == pytest.approx(0.004069614)
+
+    # Nb90 decays by β+ and should emit positrons, electrons, and photons
+    sources = nb90.sources
+    assert len(set(sources.keys()) ^ {'positron', 'electron', 'photon'}) == 0