diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst
index 4411b2144b..7091757f23 100644
--- a/docs/source/io_formats/settings.rst
+++ b/docs/source/io_formats/settings.rst
@@ -1044,8 +1044,10 @@ attributes/sub-elements:
   :emission_density:
     A list of neutron emission densities :math:`S(r)` evaluated at each
     ``r_over_a`` grid point (arbitrary units, must be non-negative). Only the
-    shape matters, since the profile is normalized internally. Must have the
-    same length as ``r_over_a``.
+    shape matters, since the profile is normalized internally. Values are
+    interpolated linearly between grid points and the profile is refined on an
+    internal grid for radial sampling. Must have the same length as
+    ``r_over_a`` and contain at least one positive value.
 
   :phi_start:
     The starting toroidal angle in [rad].
@@ -1060,7 +1062,8 @@ attributes/sub-elements:
 
   :n_alpha:
     The number of poloidal-angle grid points used to build the sampling CDFs
-    (must be > 2). Larger values reduce discretization bias.
+    (must be > 2). Larger values reduce discretization bias; values below 51
+    produce a warning.
 
     *Default*: 101
 
diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst
index e9ef404c14..5e7007a1ab 100644
--- a/docs/source/usersguide/settings.rst
+++ b/docs/source/usersguide/settings.rst
@@ -306,7 +306,9 @@ The plasma shape is described by the major radius :math:`R_0`, minor radius
 Shafranov shift :math:`\Delta`. The neutron birth profile is given as an
 emission density :math:`S(r/a)` tabulated on a normalized minor-radius grid that
 runs from 0 (magnetic axis) to 1 (last closed flux surface); only the shape of
-the profile matters, since it is normalized internally. For example::
+the profile matters, since it is normalized internally. The emission density
+is linearly interpolated between the supplied points and refined internally for
+radial sampling. For example::
 
   import numpy as np
 
diff --git a/openmc/source.py b/openmc/source.py
index 228505058d..4e9f84dfc0 100644
--- a/openmc/source.py
+++ b/openmc/source.py
@@ -46,7 +46,7 @@ class SourceBase(ABC):
 
     Attributes
     ----------
-    type : {'independent', 'file', 'compiled', 'mesh'}
+    type : {'independent', 'file', 'compiled', 'mesh', 'tokamak'}
         Indicator of source type.
     strength : float
         Strength of the source
@@ -931,7 +931,9 @@ class TokamakSource(SourceBase):
         Normalized minor radius grid points, must start at 0 and end at 1
     emission_density : numpy.ndarray
         Emission density S(r) at each r/a point (arbitrary units, must be >= 0).
-        Must have the same length as ``r_over_a``.
+        Values are linearly interpolated between grid points and refined on an
+        internal grid for radial sampling. Must have the same length as
+        ``r_over_a`` and contain at least one positive value.
     energy : openmc.stats.Univariate or Sequence[openmc.stats.Univariate]
         Energy distribution(s). Either a single distribution used at all radii,
         or one distribution per ``r_over_a`` grid point. When one distribution
@@ -948,7 +950,9 @@ class TokamakSource(SourceBase):
     phi_extent : float
         Toroidal angle extent in [rad] (default: 2π)
     n_alpha : int
-        Number of poloidal angle grid points for CDF sampling (default: 101)
+        Number of poloidal angle grid points for CDF sampling (default: 101).
+        Values below 51 produce a warning because they may introduce noticeable
+        discretization bias.
     vertical_shift : float
         Vertical shift of the plasma center in [cm] (default: 0)
     strength : float
@@ -1027,8 +1031,10 @@ class TokamakSource(SourceBase):
         self.energy = energy
         self.time = time
 
-        # Cross-field consistency checks (each setter only validates its own
-        # field, so the relationships between fields are verified here)
+        self._validate()
+
+    def _validate(self):
+        """Validate relationships between tokamak source parameters."""
         if self.minor_radius >= self.major_radius:
             raise ValueError(
                 f"minor_radius ({self.minor_radius}) must be smaller than "
@@ -1041,6 +1047,8 @@ class TokamakSource(SourceBase):
             raise ValueError(
                 f"emission_density (length {len(self.emission_density)}) must "
                 f"have the same length as r_over_a (length {len(self.r_over_a)})")
+        if not np.any(self.emission_density > 0.0):
+            raise ValueError("emission_density must contain a positive value")
         if len(self.energy) not in (1, len(self.r_over_a)):
             raise ValueError(
                 f"Number of energy distributions ({len(self.energy)}) must be "
@@ -1182,6 +1190,10 @@ class TokamakSource(SourceBase):
     def n_alpha(self, value: int):
         cv.check_type('n_alpha', value, Integral)
         cv.check_greater_than('n_alpha', value, 2)
+        if value < 51:
+            warnings.warn(
+                "n_alpha values below 51 may introduce noticeable "
+                "discretization bias in tokamak source sampling", stacklevel=2)
         self._n_alpha = value
 
     @property
@@ -1202,6 +1214,8 @@ class TokamakSource(SourceBase):
             XML element containing source data
 
         """
+        self._validate()
+
         # Geometry parameters
         ET.SubElement(element, "major_radius").text = str(self.major_radius)
         ET.SubElement(element, "minor_radius").text = str(self.minor_radius)
diff --git a/src/source.cpp b/src/source.cpp
index b5f6665dd5..8106ac1780 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -789,6 +789,10 @@ TokamakSource::TokamakSource(pugi::xml_node node) : Source(node)
   if (n_alpha_ <= 2) {
     fatal_error("TokamakSource: n_alpha must be > 2.");
   }
+  if (n_alpha_ < 51) {
+    warning("TokamakSource: n_alpha values below 51 may introduce noticeable "
+            "discretization bias in source sampling.");
+  }
   if (energy_dists_.empty()) {
     fatal_error("TokamakSource: At least one energy distribution is required.");
   }
diff --git a/tests/cpp_unit_tests/test_math.cpp b/tests/cpp_unit_tests/test_math.cpp
index 983322e4b9..7467aa1fe7 100644
--- a/tests/cpp_unit_tests/test_math.cpp
+++ b/tests/cpp_unit_tests/test_math.cpp
@@ -46,6 +46,21 @@ TEST_CASE("Test t_percentile")
   }
 }
 
+TEST_CASE("Test cylindrical Bessel functions")
+{
+  constexpr double x = 2.0;
+  constexpr double expected[] {0.22389077914123567, 0.57672480775687339,
+    0.35283402861563773, 0.12894324947440205};
+
+  for (int n = 0; n < 4; ++n) {
+    REQUIRE_THAT(openmc::cyl_bessel_j(n, x),
+      Catch::Matchers::WithinRel(expected[n], 1.0e-14));
+    REQUIRE_THAT(openmc::cyl_bessel_j(n, -x),
+      Catch::Matchers::WithinRel(
+        (n % 2 == 0 ? 1.0 : -1.0) * expected[n], 1.0e-14));
+  }
+}
+
 TEST_CASE("Test calc_pn")
 {
   // The reference solutions come from scipy.special.eval_legendre
diff --git a/tests/unit_tests/test_source_tokamak.py b/tests/unit_tests/test_source_tokamak.py
index b950efac60..ec35cbbaea 100644
--- a/tests/unit_tests/test_source_tokamak.py
+++ b/tests/unit_tests/test_source_tokamak.py
@@ -85,6 +85,7 @@ def test_tokamak_source_multiple_energies():
     (dict(r_over_a=np.linspace(0.1, 1.0, 10)), "must start at 0"),
     (dict(r_over_a=np.linspace(0.0, 0.9, 10)), "must end at 1"),
     (dict(emission_density=-np.linspace(0.0, 1.0, 10)), "cannot be negative"),
+    (dict(emission_density=np.zeros(10)), "must contain a positive value"),
 ])
 def test_tokamak_source_invalid(kwargs, match):
     with pytest.raises(ValueError, match=match):
@@ -102,6 +103,25 @@ def test_tokamak_source_invalid_n_alpha():
         make_source(n_alpha=2)
 
 
+def test_tokamak_source_low_n_alpha_warning():
+    with pytest.warns(UserWarning, match="below 51"):
+        make_source(n_alpha=50)
+
+
+@pytest.mark.parametrize(("attribute", "value", "match"), [
+    ("minor_radius", 700.0, "smaller than major_radius"),
+    ("shafranov_shift", 150.0, "half the minor_radius"),
+    ("emission_density", np.ones(5), "same length as r_over_a"),
+    ("energy", [openmc.stats.Discrete([1.0], [1.0])] * 2,
+     "Number of energy distributions"),
+])
+def test_tokamak_source_mutation_validation(attribute, value, match):
+    src = make_source()
+    setattr(src, attribute, value)
+    with pytest.raises(ValueError, match=match):
+        src.to_xml_element()
+
+
 @pytest.mark.parametrize(("emission_density", "expected_mean"), [
     ([1.0, 1.0], 2.0 / 3.0),
     ([1.0, 0.0], 0.5),