diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst
index 3b1031769..2309643dc 100644
--- a/docs/source/io_formats/statepoint.rst
+++ b/docs/source/io_formats/statepoint.rst
@@ -149,6 +149,8 @@ The current version of the statepoint file format is 18.1.
                tallies will have a value of 0 unless otherwise instructed.
              - **multiply_density** (*int*) -- Flag indicating whether reaction
                rates should be multiplied by atom density (1) or not (0).
+             - **higher_moments** (*int*) -- Flag indicating whether
+               higher-order tally moments are enabled (1) or not (0).
 
 :Datasets: - **n_realizations** (*int*) -- Number of realizations.
            - **n_filters** (*int*) -- Number of filters used.
diff --git a/docs/source/methods/tallies.rst b/docs/source/methods/tallies.rst
index 79a63fbdd..27a3f873a 100644
--- a/docs/source/methods/tallies.rst
+++ b/docs/source/methods/tallies.rst
@@ -387,6 +387,101 @@ of this is that the longer you run a simulation, the better you know your
 results. Therefore, by running a simulation long enough, it is possible to
 reduce the stochastic uncertainty to arbitrarily low levels.
 
+Skewness
+++++++++
+
+The `skewness`_ of a population quantifies the asymmetry of the probability
+distribution around its mean. Positive and negative skewness indicate a
+longer/heavier right and left tail respectively. Let :math:`x_1,\ldots,x_n` be
+the per-realization values for a bin, with sample mean :math:`\bar{x}` and
+sample central moments:
+
+.. math::
+
+   m_k \;=\; \frac{1}{n}\sum_{i=1}^{n}\bigl(x_i-\bar{x}\bigr)^k.
+
+OpenMC reports the *adjusted Fisher-Pearson skewness* (defined for :math:`n \ge
+3`), which is commonly used in many statistical packages:
+
+.. math::
+
+   G_1 \;=\; \frac{\sqrt{n \cdot (n-1)}}{\,n-2\,}\cdot\frac{m_3}{m_2^{3/2}}.
+
+where :math:`m_2` and :math:`m_3` correspond to the biased sample second and
+third central moment respectively.
+
+Kurtosis
+++++++++
+
+The `kurtosis`_ of a population quantifies tail weight (also called tailedness)
+of the probability distribution relative to a normal distribution. Positive
+excess kurtosis indicates *heavier tails* whereas negative excess kurtosis
+indicates *lighter tails*. Kurtosis is especially useful for identifying bins
+where occasional extreme scores dominate uncertainty. OpenMC reports the
+*adjusted excess kurtosis* (defined for :math:`n \ge 4`):
+
+.. math::
+
+   G_2 \;=\; \frac{(n-1)}{(n-2)(n-3)}
+             \left[(n+1)\,\frac{m_4}{m_2^{2}} \;-\; 3(n-1)\right].
+
+where :math:`m_2` and :math:`m_4` correspond to the biased sample second and
+fourth central moment respectively. For a perfectly normal distribution, the
+excess kurtosis is :math:`0`.
+
+Variance of Variance
+++++++++++++++++++++
+
+The variance of the variance (also known as the coefficient of variation
+squared) measures *stability of the sample variance* :math:`s^2` and, by
+extension, the reliability of reported relative errors. High VOV means that
+error bars themselves are noisy—often due to heavy tails, skewness, or too few
+realizations.
+
+.. math::
+
+    VOV = \frac{s^2(s_{\bar{X}}^2)}{s_{\bar{X}}^4 } = \frac{m_4}{m_2^2} - \frac{1}{n}
+
+where :math:`s_{\bar{X}}^2` is the estimated variance of the mean and
+:math:`s^2(s_{\bar{X}}^2)` is the estimated variance in :math:`s_{\bar{X}}^2`.
+The MCNP manual suggests a hard threshold such that :math:`VOV < 0.1` to improve
+the probability of forming a reliable confidence interval. However, OpenMC does
+not enforce an universal cut-off because the suitability of any single threshold
+depends strongly on problem specifics (estimator choice, variance-reduction
+settings, tally binning, or even effective sample size).
+
+
+Normality Tests (D'Agostino-Pearson)
+++++++++++++++++++++++++++++++++++++
+
+These normality test verify the hypothesis that fluctuations are *approximately
+normal*, a working assumption behind many Monte Carlo diagnostics and
+`confidence-interval heuristics`_. Tests are provided for: (i) skewness-only,
+(ii) kurtosis-only, and (iii) the *omnibus* combination. OpenMC uses the
+finite-sample-adjusted skewness :math:`G_1` and excess kurtosis :math:`G_2`
+above to construct standardized normal scores :math:`Z_1` (from :math:`G_1`) and
+:math:`Z_2` (from :math:`G_2`) via the D'Agostino-Pearson transformations. The
+omnibus statistic is
+
+.. math::
+
+   K^2 \;=\; Z_1^{\,2} \;+\; Z_2^{\,2}
+   \;\sim\; \chi^2_{(2)} \quad \text{under } H_0:\ \text{normality}.
+
+OpenMC reports :math:`Z_1`, :math:`Z_2`, :math:`K^2`, and their p-values when
+prerequisites are met (skewness for :math:`n\ge 3`, kurtosis and omnibus for
+:math:`n\ge 4`). Given a user-chosen significance level :math:`\alpha` (default
+is :math:`0.05`), reject :math:`H_0` if :math:`\text{p-value}<\alpha`; otherwise
+fail to reject. OpenMC leaves the interpretation to the user, who should
+consider VOV together with skewness, kurtosis, and normality tests results when
+judging whether reported confidence intervals are credible for their application
+[#norm-tests]_.
+
+.. [#norm-tests]
+   Higher-moments accumulation must be enabled with ``higher_moments = True``
+   for running these diagnostics including the skewness, kurtosis, and normality
+   tests.
+
 Figure of Merit
 +++++++++++++++
 
@@ -405,14 +500,16 @@ defined as
 .. math::
     :label: relative_error
 
-    r = \frac{s_\bar{X}}{\bar{x}}.
+    r = \frac{s_{\bar{X}}}{\bar{x}}.
 
 Based on this definition, one can see that a higher FOM is desirable. The FOM is
 useful as a comparative tool. For example, if a variance reduction technique is
 being applied to a simulation, the FOM with variance reduction can be compared
 to the FOM without variance reduction to ascertain whether the reduction in
 variance outweighs the potential increase in execution time (e.g., due to
-particle splitting).
+particle splitting). It is important to note that MCNP reports the FOM using CPU
+time (wall-clock time multiplied by the number of threads/cores), whereas OpenMC
+reports the FOM using only the wall-clock time :math:`t`.
 
 Confidence Intervals
 ++++++++++++++++++++
@@ -521,6 +618,8 @@ improve the estimate of the percentile.
 
    .. rubric:: References
 
+.. _confidence-interval heuristics: https://doi.org/10.1080/00031305.1990.10475751
+
 .. _following approximation: https://doi.org/10.1080/03610918708812641
 
 .. _Bessel's correction: https://en.wikipedia.org/wiki/Bessel's_correction
@@ -541,6 +640,10 @@ improve the estimate of the percentile.
 
 .. _converges in distribution: https://en.wikipedia.org/wiki/Convergence_of_random_variables#Convergence_in_distribution
 
+.. _skewness: https://en.wikipedia.org/wiki/Skewness
+
+.. _kurtosis: https://en.wikipedia.org/wiki/Kurtosis
+
 .. _confidence intervals: https://en.wikipedia.org/wiki/Confidence_interval
 
 .. _Student's t-distribution: https://en.wikipedia.org/wiki/Student%27s_t-distribution
diff --git a/include/openmc/constants.h b/include/openmc/constants.h
index a0d164613..b66193481 100644
--- a/include/openmc/constants.h
+++ b/include/openmc/constants.h
@@ -291,7 +291,7 @@ enum class MgxsType {
 // ============================================================================
 // TALLY-RELATED CONSTANTS
 
-enum class TallyResult { VALUE, SUM, SUM_SQ, SIZE };
+enum class TallyResult { VALUE, SUM, SUM_SQ, SUM_THIRD, SUM_FOURTH };
 
 enum class TallyType { VOLUME, MESH_SURFACE, SURFACE, PULSE_HEIGHT };
 
diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h
index 0092c08f8..28b0d2b11 100644
--- a/include/openmc/hdf5_interface.h
+++ b/include/openmc/hdf5_interface.h
@@ -100,8 +100,8 @@ void read_llong(hid_t obj_id, const char* name, long long* buffer, bool indep);
 void read_string(
   hid_t obj_id, const char* name, size_t slen, char* buffer, bool indep);
 
-void read_tally_results(
-  hid_t group_id, hsize_t n_filter, hsize_t n_score, double* results);
+void read_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score,
+  hsize_t n_results, double* results);
 void write_attr_double(hid_t obj_id, int ndim, const hsize_t* dims,
   const char* name, const double* buffer);
 void write_attr_int(hid_t obj_id, int ndim, const hsize_t* dims,
@@ -114,9 +114,9 @@ void write_int(hid_t group_id, int ndim, const hsize_t* dims, const char* name,
 void write_llong(hid_t group_id, int ndim, const hsize_t* dims,
   const char* name, const long long* buffer, bool indep);
 void write_string(hid_t group_id, int ndim, const hsize_t* dims, size_t slen,
-  const char* name, char const* buffer, bool indep);
-void write_tally_results(
-  hid_t group_id, hsize_t n_filter, hsize_t n_score, const double* results);
+  const char* name, const char* buffer, bool indep);
+void write_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score,
+  hsize_t n_results, const double* results);
 } // extern "C"
 
 //==============================================================================
diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h
index 3beeb9d5a..374daff92 100644
--- a/include/openmc/tallies/tally.h
+++ b/include/openmc/tallies/tally.h
@@ -106,6 +106,8 @@ public:
 
   bool writable() const { return writable_; }
 
+  bool higher_moments() const { return higher_moments_; }
+
   //----------------------------------------------------------------------------
   // Other methods.
 
@@ -190,6 +192,9 @@ private:
   //! Whether to multiply by atom density for reaction rates
   bool multiply_density_ {true};
 
+  //! Whether to accumulate higher moments (third and fourth)
+  bool higher_moments_ {false};
+
   int64_t index_;
 };
 
diff --git a/openmc/statepoint.py b/openmc/statepoint.py
index a763db397..11986841f 100644
--- a/openmc/statepoint.py
+++ b/openmc/statepoint.py
@@ -434,6 +434,10 @@ class StatePoint:
                     if "multiply_density" in group.attrs:
                         tally.multiply_density = group.attrs["multiply_density"].item() > 0
 
+                    # Check if tally has higher_moments attribute
+                    if 'higher_moments' in group.attrs:
+                        tally.higher_moments = bool(group.attrs['higher_moments'][()])
+
                     # Read the number of realizations
                     n_realizations = group['n_realizations'][()]
 
diff --git a/openmc/tallies.py b/openmc/tallies.py
index 075b1e991..25ec29a58 100644
--- a/openmc/tallies.py
+++ b/openmc/tallies.py
@@ -3,6 +3,7 @@ from collections.abc import Iterable, MutableSequence
 import copy
 from functools import partial, reduce, wraps
 from itertools import product
+from math import sqrt, log
 from numbers import Integral, Real
 import operator
 from pathlib import Path
@@ -12,6 +13,7 @@ import h5py
 import numpy as np
 import pandas as pd
 import scipy.sparse as sps
+from scipy.stats import chi2, norm
 
 import openmc
 import openmc.checkvalue as cv
@@ -91,10 +93,20 @@ class Tally(IDManagerMixin):
     sum_sq : numpy.ndarray
         An array containing the sum of each independent realization squared for
         each bin
+    sum_third : numpy.ndarray
+        An array containing the sum of each independent realization to the third power for
+        each bin
+    sum_fourth : numpy.ndarray
+        An array containing the sum of each independent realization to the fourth power for
+        each bin
     mean : numpy.ndarray
         An array containing the sample mean for each bin
     std_dev : numpy.ndarray
         An array containing the sample standard deviation for each bin
+    vov : numpy.ndarray
+        An array containing the variance of the variance for each tally bin
+    higher_moments : bool
+        Whether or not the tally accumulates the sums third and fourth to compute higher-order moments
     figure_of_merit : numpy.ndarray
         An array containing the figure of merit for each bin
 
@@ -129,8 +141,12 @@ class Tally(IDManagerMixin):
 
         self._sum = None
         self._sum_sq = None
+        self._sum_third = None
+        self._sum_fourth = None
         self._mean = None
         self._std_dev = None
+        self._vov = None
+        self._higher_moments = False
         self._simulation_time = None
         self._with_batch_statistics = False
         self._derived = False
@@ -221,6 +237,15 @@ class Tally(IDManagerMixin):
         cv.check_type('multiply density', value, bool)
         self._multiply_density = value
 
+    @property
+    def higher_moments(self) -> bool:
+        return self._higher_moments
+
+    @higher_moments.setter
+    def higher_moments(self, value):
+        cv.check_type("higher_moments", value, bool)
+        self._higher_moments = value
+
     @property
     def filters(self):
         return self._filters
@@ -371,6 +396,11 @@ class Tally(IDManagerMixin):
             # Update nuclides
             nuclide_names = group['nuclides'][()]
             self._nuclides = [name.decode().strip() for name in nuclide_names]
+            # Check for higher_moments attribute
+            if "higher_moments" in group.attrs:
+                self._higher_moments = bool(group.attrs["higher_moments"][()])
+            else:
+                self._higher_moments = False
 
             # Extract Tally data from the file
             data = group['results']
@@ -385,10 +415,25 @@ class Tally(IDManagerMixin):
             self._sum = sum_
             self._sum_sq = sum_sq
 
+            if self._higher_moments:
+                # Extract additional Tally data when higher moments enabled
+                sum_third = data[:, :, 2]
+                sum_fourth = data[:, :, 3]
+
+                # Reshape the results arrays
+                sum_third = np.reshape(sum_third, self.shape)
+                sum_fourth = np.reshape(sum_fourth, self.shape)
+
+                # Set the additional data for this Tally
+                self._sum_third = sum_third
+                self._sum_fourth = sum_fourth
+
             # Convert NumPy arrays to SciPy sparse LIL matrices
             if self.sparse:
                 self._sum = sps.lil_matrix(self._sum.flatten(), self._sum.shape)
                 self._sum_sq = sps.lil_matrix(self._sum_sq.flatten(), self._sum_sq.shape)
+                self._sum_third = sps.lil_matrix(self._sum_third.flatten(), self._sum_third.shape)
+                self._sum_fourth = sps.lil_matrix(self.sum_fourth.flatten(), self._sum_fourth.shape)
 
             # Read simulation time (needed for figure of merit)
             self._simulation_time = f["runtime"]["simulation"][()]
@@ -428,6 +473,52 @@ class Tally(IDManagerMixin):
         cv.check_type('sum_sq', sum_sq, Iterable)
         self._sum_sq = sum_sq
 
+    @property
+    @ensure_results
+    def sum_third(self):
+        if not self._higher_moments:
+            raise ValueError(
+                "Higher moments have not been enabled for this tally. To make "
+                "higher moments available, set the higher_moments attribute to "
+                "True before running a simulation."
+            )
+
+        if not self._sp_filename or self.derived:
+            return None
+
+        if self.sparse:
+            return np.reshape(self._sum_third.toarray(), self.shape)
+        else:
+            return self._sum_third
+
+    @sum_third.setter
+    def sum_third(self, sum_third):
+        cv.check_type("sum_third", sum_third, Iterable)
+        self._sum_third = sum_third
+
+    @property
+    @ensure_results
+    def sum_fourth(self):
+        if not self._higher_moments:
+            raise ValueError(
+                "Higher moments have not been enabled for this tally. To make "
+                "higher moments available, set the higher_moments attribute to "
+                "True before running a simulation."
+            )
+
+        if not self._sp_filename or self.derived:
+            return None
+
+        if self.sparse:
+            return np.reshape(self._sum_fourth.toarray(), self.shape)
+        else:
+            return self._sum_fourth
+
+    @sum_fourth.setter
+    def sum_fourth(self, sum_fourth):
+        cv.check_type("sum_fourth", sum_fourth, Iterable)
+        self._sum_fourth = sum_fourth
+
     @property
     def mean(self):
         if self._mean is None:
@@ -470,14 +561,370 @@ class Tally(IDManagerMixin):
         else:
             return self._std_dev
 
+    @property
+    def vov(self):
+        if self._vov is None:
+            n = self.num_realizations
+            sum1 = self.sum
+            sum2 = self.sum_sq
+            sum3 = self.sum_third
+            sum4 = self.sum_fourth
+            self._vov = np.zeros_like(sum1, dtype=float)
+
+            # Calculate the variance of the variance (Eq. 2.232 in
+            # https://doi.org/10.2172/2372634)
+            numerator = (sum4 - (4.0*sum3*sum1)/n
+                        + (6.0*sum2*(sum1**2))/(n**2)
+                        - (3.0*(sum1)**4)/(n**3))
+            denominator = (sum2 - (1.0/n)*(sum1**2))**2
+
+            mask = denominator > 0.0
+
+            self._vov[mask] = numerator[mask]/denominator[mask] - 1.0/n
+
+            if self.sparse:
+                self._vov = sps.lil_matrix(self._vov.flatten(), self._vov.shape)
+
+        if self.sparse:
+            return np.reshape(self._vov.toarray(), self.shape)
+        else:
+            return self._vov
+
+    @property
+    def m2(self):
+        n = self.num_realizations
+        return self.sum_sq/n - self.mean**2
+
+    @property
+    def m3(self):
+        n = self.num_realizations
+        mean = self.mean
+        sum2 = self.sum_sq/n
+        sum3 = self.sum_third/n
+
+        return sum3 - 3.0*mean*sum2 + 2.0*mean**3
+
+    @property
+    def m4(self):
+        n = self.num_realizations
+        mean = self.mean
+        sum2 = self.sum_sq/n
+        sum3 = self.sum_third/n
+        sum4 = self.sum_fourth/n
+
+        return sum4 - 4.0*mean*sum3 + 6.0*(mean**2)*sum2 - 3.0*mean**4
+
+    def skew(self, bias=False) -> np.ndarray:
+        """Return the sample skewness of each tally bin.
+
+        This method computes and returns the unadjusted or adjusted
+        Fisher-Pearson coefficient of skewness.
+
+        Parameters
+        ----------
+        bias : bool
+            If False, calculations are corrected for bias and the adjusted
+            Fisher-Pearson skewness (:math:`G_1`) is returned. If True,
+            calculations are not corrected for bias and the unadjusted skewness
+            (:math:`g_1`) is returned.
+
+        Returns
+        -------
+        float
+            The skewness of each tally bin
+        """
+        n = self.num_realizations
+        m2 = self.m2
+        m3 = self.m3
+
+        with np.errstate(divide="ignore", invalid="ignore"):
+            g1 = np.where(m2 > 0.0, m3/(m2**1.5), 0.0)
+
+        if bias:
+            return g1
+        else:
+            if n <= 2:
+                raise ValueError("Insufficient number of independent realizations"
+                                 f"for bias-corrected skewness: need n >= 3, got {n=}.")
+            else:
+                return sqrt(n*(n - 1))/(n - 2)*g1
+
+    def kurtosis(self, fisher=True, bias=False) -> np.ndarray:
+        r"""Return the sample kurtosis of each tally bin.
+
+        This method computes and returns the sample kurtosis using either
+        Pearson's or Fisher's definition, with or without finite-sample bias
+        correction. The value returned depends on the `bias` and `fisher`
+        arguments as follows:
+
+        - **bias=True, fisher=False**: Returns :math:`b_2` (Pearson's kurtosis)
+          This is the raw fourth standardized moment: :math:`m_4/m_2^2`. For a
+          normal distribution, :math:`b_2\approx 3`.
+
+        - **bias=True, fisher=True**: Returns :math:`g_2` (excess kurtosis) This
+          is :math:`b_2 - 3`, centered at 0 for normal distributions. Positive
+          values indicate heavier tails, negative values lighter tails.
+
+        - **bias=False, fisher=True** (default): Returns :math:`G_2` (adjusted
+          excess kurtosis). This applies finite-sample bias correction to
+          :math:`g_2`. This is the recommended estimator for statistical
+          inference.
+
+        - **bias=False, fisher=False**: Returns bias-corrected Pearson's
+          kurtosis. This is :math:`G_2 + 3`.
+
+        Parameters
+        ----------
+        fisher : bool, optional
+            If True (default), Fisher's definition is used (excess kurtosis). If
+            False, Pearson's definition is used.
+        bias : bool, optional
+            If False (default), calculations are corrected for statistical bias
+            using finite-sample adjustments. If True, calculations use the
+            biased estimator (population formulas).
+
+        Returns
+        -------
+        numpy.ndarray
+            The kurtosis of each tally bin
+
+        """
+        n = self.num_realizations
+        m2 = self.m2
+        m4 = self.m4
+
+        with np.errstate(divide="ignore", invalid="ignore"):
+            b2 = np.where(m2 > 0.0, m4/(m2**2), 0.0)
+        g2 = b2 - 3.0
+
+        if bias:
+            # Biased estimator (g2 or b2)
+            return g2 if fisher else b2
+        else:
+            # Unbiased estimator with finite-sample correction
+            if n <= 3:
+                raise ValueError("Insufficient number of independent realizations"
+                                 f"for bias-corrected kurtosis: need n >= 4, got {n=}.")
+            else:
+                G2 = ((n - 1)/((n - 2)*(n - 3)))*((n + 1)*g2 + 6.0)
+                return G2 if fisher else G2 + 3.0
+
+    def skewtest(self, alternative: str = "two-sided"):
+        """Perform D'Agostino and Pearson's test for skewness.
+
+        This method tests the null hypothesis that the skewness of the
+        population that the sample was drawn from is the same as that of a
+        corresponding normal distribution.
+
+        Parameters
+        ----------
+        alternative : {'two-sided', 'less', 'greater'}, optional
+            Defines the alternative hypothesis. The following options are
+            available:
+
+            * 'two-sided': the skewness of the distribution is different from
+              that of the normal distribution (i.e., non-zero)
+            * 'less': the skewness of the distribution is less than that of the
+              normal distribution
+            * 'greater': the skewness of the distribution is greater than that
+              of the normal distribution
+
+        Returns
+        -------
+        statistic : np.ndarray
+            The computed z-score for the skewness test for each tally bin
+        pvalue : np.ndarray
+            The p-value for the hypothesis test for each tally bin
+
+        Notes
+        -----
+        This test is based on D'Agostino and Pearson's test [1]_. The test
+        requires at least 8 realizations to produce valid results.
+
+        References
+        ----------
+        .. [1] D'Agostino, R. B. (1971), "An omnibus test of normality for
+               moderate and large sample size", Biometrika, 58, 341-348
+
+        """
+        n = self.num_realizations
+        if n < 8:
+            raise ValueError("Skewness test is not well-defined for n < 8.")
+
+        g1 = self.skew(bias=True)
+
+        # --- Z1 (skewness) ---
+        y = g1 * sqrt(((n + 1.0)*(n + 3.0))/(6.0*(n - 2.0)))
+        beta2 = (3.0*(n**2 + 27.0*n - 70.0)*(n + 1.0)*(n + 3.0)
+                )/((n - 2.0)*(n + 5.0)*(n + 7.0)*(n + 9.0))
+        W2 = -1.0 + sqrt(2.0*(beta2 - 1.0))
+        delta = 1.0 / sqrt(log(sqrt(W2)))
+        alpha = sqrt(2.0 / (W2 - 1.0))
+        Zb1 = np.where(
+            y >= 0.0,
+            delta*np.log((y/alpha) + np.sqrt((y/alpha)**2 + 1.0)),
+            -delta*np.log((-y/alpha) + np.sqrt((y/alpha)**2 + 1.0))
+        )
+
+        # p-value
+        if alternative == "two-sided":
+            p = 2.0 * (1.0 - norm.cdf(np.abs(Zb1)))
+        elif alternative == "greater":
+            p = 1.0 - norm.cdf(Zb1)
+        elif alternative == "less":
+            p = norm.cdf(Zb1)
+        else:
+            raise ValueError("alternative must be 'two-sided', 'greater', or 'less'")
+
+        return Zb1, p
+
+    def kurtosistest(self, alternative: str = "two-sided"):
+        """Perform D'Agostino and Pearson's test for kurtosis.
+
+        This method tests the null hypothesis that the kurtosis of the
+        population that the sample was drawn from is the same as that of a
+        corresponding normal distribution.
+
+        Parameters
+        ----------
+        alternative : {'two-sided', 'less', 'greater'}, optional
+            Defines the alternative hypothesis. Default is 'two-sided'.
+            The following options are available:
+
+            * 'two-sided': the kurtosis of the distribution is different from
+              that of the normal distribution
+            * 'less': the kurtosis of the distribution is less than that of the
+              normal distribution
+            * 'greater': the kurtosis of the distribution is greater than that
+              of the normal distribution
+
+        Returns
+        -------
+        statistic : np.ndarray
+            The computed z-score for the kurtosis test for each tally bin
+        pvalue : np.ndarray
+            The p-value for the hypothesis test for each tally bin
+
+        Raises
+        ------
+        ValueError
+            If the number of realizations is less than 20, or if an invalid
+            alternative hypothesis is specified.
+
+        Notes
+        -----
+        This test is based on D'Agostino and Pearson's test [1]_. The test
+        is typically recommended for at least 20 realizations to produce
+        valid results.
+
+        References
+        ----------
+        .. [1] D'Agostino, R. B. (1971), "An omnibus test of normality for
+               moderate and large sample size", Biometrika, 58, 341-348
+
+        """
+        n = self.num_realizations
+        if n < 20:
+            raise ValueError("Kurtosis test is typically recommended for n >= 20.")
+
+        b2 = self.kurtosis(bias=True, fisher=False)
+
+        # --- Z2 (kurtosis) ---
+        mean_b2 = 3.0 * (n - 1.0) / (n + 1.0)
+        var_b2 = (24.0*n*(n - 2.0)*(n - 3.0)/(
+                    (n + 1.0)**2*(n + 3.0)*(n + 5.0)))
+        x = (b2 - mean_b2)/np.sqrt(var_b2)
+        moment = ((6.0*(n**2 - 5.0*n + 2.0))/((n + 7.0)*(n + 9.0))
+                    )*sqrt((6.0*(n + 3.0)*(n + 5.0))/(n*(n - 2.0)*(n - 3.0)))
+        A = 6.0 + (8.0/moment)*((2.0/moment) + sqrt(1.0 + 4.0/(moment**2)))
+        Zb2 = (1.0- 2.0/(9.0*A) - ((1.0 - 2.0/A) / (1.0 + (x
+                )*sqrt(2.0/(A - 4.0))))**(1.0/3.0)) / sqrt(2.0/(9.0*A))
+
+        # p-value
+        if alternative == "two-sided":
+            p = 2.0 * (1.0 - norm.cdf(np.abs(Zb2)))
+        elif alternative == "greater":
+            p = 1.0 - norm.cdf(Zb2)
+        elif alternative == "less":
+            p = norm.cdf(Zb2)
+        else:
+            raise ValueError("alternative must be 'two-sided', 'greater', or 'less'")
+
+        return Zb2, p
+
+    def normaltest(self, alternative: str = "two-sided"):
+        """Perform D'Agostino and Pearson's omnibus test for normality.
+
+        This method tests the null hypothesis that a sample comes from a
+        normal distribution. It combines skewness and kurtosis to produce an
+        omnibus test of normality.
+
+        Parameters
+        ----------
+        alternative : {'two-sided', 'less', 'greater'}, optional
+            Defines the alternative hypothesis used for the component skewness
+            and kurtosis tests. Default is 'two-sided'. The following options
+            are available:
+
+            * 'two-sided': the distribution is different from normal
+            * 'less': used for the component tests
+            * 'greater': used for the component tests
+
+        Returns
+        -------
+        statistic : np.ndarray
+            The computed z-score for the normality test for each tally bin
+        pvalue : np.ndarray
+            The p-value for the hypothesis test for each tally bin
+
+        Raises
+        ------
+        ValueError
+            If the number of realizations is less than 20, or if an invalid
+            alternative hypothesis is specified.
+
+        Notes
+        -----
+        This test combines a test for skewness and a test for kurtosis to
+        produce an omnibus test [1]_. The test statistic is:
+
+        .. math::
+
+            K^2 = Z_1^2 + Z_2^2
+
+        where :math:`Z_1` is the z-score from the skewness test and
+        :math:`Z_2` is the z-score from the kurtosis test. This statistic
+        follows a chi-square distribution with 2 degrees of freedom.
+
+        The test requires at least 20 realizations to produce valid results.
+
+        References
+        ----------
+        .. [1] D'Agostino, R. B. and Pearson, E. S. (1973), "Tests for
+               departure from normality", Biometrika, 60, 613-622
+
+        """
+        n = self.num_realizations
+        if n < 20:
+            raise ValueError("normaltest requires n >= 20 (per D'Agostino-Pearson).")
+
+        # Use the component tests
+        Z1, _ = self.skewtest(alternative)
+        Z2, _ = self.kurtosistest(alternative)
+
+        # Combine as chi-square with df=2 since we have skewness and kurtosis
+        K2 = Z1*Z1 + Z2*Z2
+        p = chi2.sf(K2, df=2)
+        return K2, p
+
     @property
     def figure_of_merit(self):
         mean = self.mean
         std_dev = self.std_dev
         fom = np.zeros_like(mean)
         nonzero = np.abs(mean) > 0
-        fom[nonzero] = 1.0 / (
-            (std_dev[nonzero] / mean[nonzero])**2 * self._simulation_time)
+        rel_err = std_dev[nonzero] / mean[nonzero]
+        fom[nonzero] = 1.0 / (rel_err**2 * self._simulation_time)
         return fom
 
     @property
@@ -528,6 +975,12 @@ class Tally(IDManagerMixin):
             if self._sum_sq is not None:
                 self._sum_sq = sps.lil_matrix(self._sum_sq.flatten(),
                                               self._sum_sq.shape)
+            if self._sum_third is not None:
+                self._sum_third = sps.lil_matrix(self._sum_third.flatten(),
+                                                 self._sum_third.shape)
+            if self._sum_fourth is not None:
+                self._sum_fourth = sps.lil_matrix(self._sum_fourth.flatten(),
+                                                  self._sum_fourth.shape)
             if self._mean is not None:
                 self._mean = sps.lil_matrix(self._mean.flatten(),
                                             self._mean.shape)
@@ -543,6 +996,10 @@ class Tally(IDManagerMixin):
                 self._sum = np.reshape(self._sum.toarray(), self.shape)
             if self._sum_sq is not None:
                 self._sum_sq = np.reshape(self._sum_sq.toarray(), self.shape)
+            if self._sum_third is not None:
+                self._sum_third = np.reshape(self._sum_third.toarray(), self.shape)
+            if self._sum_fourth is not None:
+                self._sum_fourth = np.reshape(self._sum_fourth.toarray(), self.shape)
             if self._mean is not None:
                 self._mean = np.reshape(self._mean.toarray(), self.shape)
             if self._std_dev is not None:
@@ -869,6 +1326,34 @@ class Tally(IDManagerMixin):
 
             merged_tally._sum_sq = np.reshape(merged_sum_sq, merged_tally.shape)
 
+        # Concatenate sum_third arrays if present in both tallies
+        if self._sum_third is not None and other._sum_third is not None:
+            self_sum_third = self.get_reshaped_data(value="sum_third")
+            other_sum_third = other_copy.get_reshaped_data(value="sum_third")
+
+            if join_right:
+                merged_sum_third = np.concatenate((self_sum_third, other_sum_third),
+                                                  axis=merge_axis)
+            else:
+                merged_sum_third = np.concatenate((other_sum_third, self_sum_third),
+                                                  axis=merge_axis)
+
+            merged_tally._sum_third = np.reshape(merged_sum_third, merged_tally.shape)
+
+        # Concatenate sum_fourth arrays if present in both tallies
+        if self._sum_fourth is not None and other._sum_fourth is not None:
+            self_sum_fourth = self.get_reshaped_data(value="sum_fourth")
+            other_sum_fourth = other_copy.get_reshaped_data(value="sum_fourth")
+
+            if join_right:
+                merged_sum_fourth = np.concatenate((self_sum_fourth, other_sum_fourth),
+                                                   axis=merge_axis)
+            else:
+                merged_sum_fourth = np.concatenate((other_sum_fourth, self_sum_fourth),
+                                                   axis=merge_axis)
+
+            merged_tally._sum_fourth = np.reshape(merged_sum_fourth, merged_tally.shape)
+
         # Concatenate mean arrays if present in both tallies
         if self.mean is not None and other.mean is not None:
             self_mean = self.get_reshaped_data(value='mean')
@@ -958,6 +1443,11 @@ class Tally(IDManagerMixin):
             subelement = ET.SubElement(element, "derivative")
             subelement.text = str(self.derivative.id)
 
+        # Optional higher moments accumulation
+        if self.higher_moments:
+            subelement = ET.SubElement(element, "higher_moments")
+            subelement.text = str(self.higher_moments).lower()
+
         return element
 
     def add_results(self, statepoint: cv.PathLike | openmc.StatePoint):
@@ -984,8 +1474,12 @@ class Tally(IDManagerMixin):
         # point are based on the current statepoint file
         self._sum = None
         self._sum_sq = None
+        self._sum_third = None
+        self._sum_fourth = None
         self._mean = None
         self._std_dev = None
+        self._vov = None
+        self._higher_moments = False
         self._num_realizations = 0
         self._results_read = False
 
@@ -1355,7 +1849,9 @@ class Tally(IDManagerMixin):
            (value == 'std_dev' and self.std_dev is None) or \
            (value == 'rel_err' and self.mean is None) or \
            (value == 'sum' and self.sum is None) or \
-           (value == 'sum_sq' and self.sum_sq is None):
+           (value == 'sum_sq' and self.sum_sq is None) or \
+           (value == "sum_third" and self.sum_third is None) or \
+           (value == "sum_fourth" and self.sum_fourth is None):
             msg = f'The Tally ID="{self.id}" has no data to return'
             raise ValueError(msg)
 
@@ -1378,10 +1874,14 @@ class Tally(IDManagerMixin):
             data = self.sum[indices]
         elif value == 'sum_sq':
             data = self.sum_sq[indices]
+        elif value == "sum_third":
+            data = self.sum_third[indices]
+        elif value == "sum_fourth":
+            data = self.sum_fourth[indices]
         else:
             msg = f'Unable to return results from Tally ID="{value}" since ' \
                   f'the requested value "{self.id}" is not \'mean\', ' \
-                  '\'std_dev\', \'rel_err\', \'sum\', or \'sum_sq\''
+                  '\'std_dev\', \'rel_err\', \'sum\', \'sum_sq\', \'sum_third\' or \'sum_fourth\''
             raise LookupError(msg)
 
         return data
@@ -2711,6 +3211,16 @@ class Tally(IDManagerMixin):
             new_sum_sq = self.get_values(scores, filters, filter_bins,
                                          nuclides, 'sum_sq')
             new_tally.sum_sq = new_sum_sq
+        if not self.derived and self._sum_third is not None:
+            new_sum_third = self.get_values(
+                scores, filters, filter_bins, nuclides, "sum_third"
+            )
+            new_tally._sum_third = new_sum_third
+        if not self.derived and self._sum_fourth is not None:
+            new_sum_fourth = self.get_values(
+                scores, filters, filter_bins, nuclides, "sum_fourth"
+            )
+            new_tally._sum_fourth = new_sum_fourth
         if self.mean is not None:
             new_mean = self.get_values(scores, filters, filter_bins,
                                        nuclides, 'mean')
@@ -3151,6 +3661,12 @@ class Tally(IDManagerMixin):
         if not self.derived and self.sum_sq is not None:
             new_tally._sum_sq = np.zeros(new_tally.shape, dtype=np.float64)
             new_tally._sum_sq[diag_indices, :, :] = self.sum_sq
+        if not self.derived and self._sum_third is not None:
+            new_tally._sum_third = np.zeros(new_tally.shape, dtype=np.float64)
+            new_tally._sum_third[diag_indices, :, :] = self.sum_third
+        if not self.derived and self._sum_fourth is not None:
+            new_tally._sum_fourth = np.zeros(new_tally.shape, dtype=np.float64)
+            new_tally._sum_fourth[diag_indices, :, :] = self.sum_fourth
         if self.mean is not None:
             new_tally._mean = np.zeros(new_tally.shape, dtype=np.float64)
             new_tally._mean[diag_indices, :, :] = self.mean
diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp
index bf1f79549..c56d485e2 100644
--- a/src/hdf5_interface.cpp
+++ b/src/hdf5_interface.cpp
@@ -536,14 +536,14 @@ void read_complex(
   H5Tclose(complex_id);
 }
 
-void read_tally_results(
-  hid_t group_id, hsize_t n_filter, hsize_t n_score, double* results)
+void read_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score,
+  hsize_t n_results, double* results)
 {
   // Create dataspace for hyperslab in memory
   constexpr int ndim = 3;
-  hsize_t dims[ndim] {n_filter, n_score, 3};
+  hsize_t dims[ndim] {n_filter, n_score, n_results};
   hsize_t start[ndim] {0, 0, 1};
-  hsize_t count[ndim] {n_filter, n_score, 2};
+  hsize_t count[ndim] {n_filter, n_score, n_results - 1};
   hid_t memspace = H5Screate_simple(ndim, dims, nullptr);
   H5Sselect_hyperslab(memspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
 
@@ -686,15 +686,15 @@ void write_string(
     group_id, 0, nullptr, buffer.length(), name, buffer.c_str(), indep);
 }
 
-void write_tally_results(
-  hid_t group_id, hsize_t n_filter, hsize_t n_score, const double* results)
+void write_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score,
+  hsize_t n_results, const double* results)
 {
   // Set dimensions of sum/sum_sq hyperslab to store
   constexpr int ndim = 3;
-  hsize_t count[ndim] {n_filter, n_score, 2};
+  hsize_t count[ndim] {n_filter, n_score, n_results - 1};
 
   // Set dimensions of results array
-  hsize_t dims[ndim] {n_filter, n_score, 3};
+  hsize_t dims[ndim] {n_filter, n_score, n_results};
   hsize_t start[ndim] {0, 0, 1};
   hid_t memspace = H5Screate_simple(ndim, dims, nullptr);
   H5Sselect_hyperslab(memspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
diff --git a/src/state_point.cpp b/src/state_point.cpp
index 8195c4865..0b0fed132 100644
--- a/src/state_point.cpp
+++ b/src/state_point.cpp
@@ -201,6 +201,12 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source)
           write_attribute(tally_group, "multiply_density", 0);
         }
 
+        if (tally->higher_moments()) {
+          write_attribute(tally_group, "higher_moments", 1);
+        } else {
+          write_attribute(tally_group, "higher_moments", 0);
+        }
+
         if (tally->estimator_ == TallyEstimator::ANALOG) {
           write_dataset(tally_group, "estimator", "analog");
         } else if (tally->estimator_ == TallyEstimator::TRACKLENGTH) {
@@ -264,12 +270,13 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source)
         for (const auto& tally : model::tallies) {
           if (!tally->writable_)
             continue;
-          // Write sum and sum_sq for each bin
+
+          // Write results for each bin
           std::string name = "tally " + std::to_string(tally->id_);
           hid_t tally_group = open_group(tallies_group, name.c_str());
           auto& results = tally->results_;
           write_tally_results(tally_group, results.shape()[0],
-            results.shape()[1], results.data());
+            results.shape()[1], results.shape()[2], results.data());
           close_group(tally_group);
         }
       } else {
@@ -509,7 +516,8 @@ extern "C" int openmc_statepoint_load(const char* filename)
         } else {
           auto& results = tally->results_;
           read_tally_results(tally_group, results.shape()[0],
-            results.shape()[1], results.data());
+            results.shape()[1], results.shape()[2], results.data());
+
           read_dataset(tally_group, "n_realizations", tally->n_realizations_);
           close_group(tally_group);
         }
@@ -1001,7 +1009,8 @@ void write_tally_results_nr(hid_t file_id)
 
       // Write reduced tally results to file
       auto shape = results_copy.shape();
-      write_tally_results(tally_group, shape[0], shape[1], results_copy.data());
+      write_tally_results(
+        tally_group, shape[0], shape[1], shape[2], results_copy.data());
 
       close_group(tally_group);
     } else {
diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp
index 12ee3427e..9daeb1d69 100644
--- a/src/tallies/tally.cpp
+++ b/src/tallies/tally.cpp
@@ -107,6 +107,9 @@ Tally::Tally(pugi::xml_node node)
     multiply_density_ = get_node_value_bool(node, "multiply_density");
   }
 
+  if (check_for_node(node, "higher_moments")) {
+    higher_moments_ = get_node_value_bool(node, "higher_moments");
+  }
   // =======================================================================
   // READ DATA FOR FILTERS
 
@@ -800,7 +803,11 @@ void Tally::init_triggers(pugi::xml_node node)
 void Tally::init_results()
 {
   int n_scores = scores_.size() * nuclides_.size();
-  results_ = xt::empty<double>({n_filter_bins_, n_scores, 3});
+  if (higher_moments_) {
+    results_ = xt::empty<double>({n_filter_bins_, n_scores, 5});
+  } else {
+    results_ = xt::empty<double>({n_filter_bins_, n_scores, 3});
+  }
 }
 
 void Tally::reset()
@@ -838,14 +845,33 @@ void Tally::accumulate()
       norm = 1.0;
     }
 
-// Accumulate each result
+    // Accumulate each result
+    if (higher_moments_) {
 #pragma omp parallel for
-    for (int i = 0; i < results_.shape()[0]; ++i) {
-      for (int j = 0; j < results_.shape()[1]; ++j) {
-        double val = results_(i, j, TallyResult::VALUE) * norm;
-        results_(i, j, TallyResult::VALUE) = 0.0;
-        results_(i, j, TallyResult::SUM) += val;
-        results_(i, j, TallyResult::SUM_SQ) += val * val;
+      // filter bins (specific cell, energy bins)
+      for (int i = 0; i < results_.shape()[0]; ++i) {
+        // score bins (flux, total reaction rate, fission reaction rate, etc.)
+        for (int j = 0; j < results_.shape()[1]; ++j) {
+          double val = results_(i, j, TallyResult::VALUE) * norm;
+          double val2 = val * val;
+          results_(i, j, TallyResult::VALUE) = 0.0;
+          results_(i, j, TallyResult::SUM) += val;
+          results_(i, j, TallyResult::SUM_SQ) += val2;
+          results_(i, j, TallyResult::SUM_THIRD) += val2 * val;
+          results_(i, j, TallyResult::SUM_FOURTH) += val2 * val2;
+        }
+      }
+    } else {
+#pragma omp parallel for
+      // filter bins (specific cell, energy bins)
+      for (int i = 0; i < results_.shape()[0]; ++i) {
+        // score bins (flux, total reaction rate, fission reaction rate, etc.)
+        for (int j = 0; j < results_.shape()[1]; ++j) {
+          double val = results_(i, j, TallyResult::VALUE) * norm;
+          results_(i, j, TallyResult::VALUE) = 0.0;
+          results_(i, j, TallyResult::SUM) += val;
+          results_(i, j, TallyResult::SUM_SQ) += val * val;
+        }
       }
     }
   }
diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp
index 674fae49c..9333800bc 100644
--- a/src/weight_windows.cpp
+++ b/src/weight_windows.cpp
@@ -547,8 +547,10 @@ void WeightWindows::update_weights(const Tally* tally, const std::string& value,
 
   // build a shape for a view of the tally results, this will always be
   // dimension 5 (3 filter dimensions, 1 score dimension, 1 results dimension)
-  std::array<int, 5> shape = {
-    1, 1, 1, tally->n_scores(), static_cast<int>(TallyResult::SIZE)};
+  // Look for the size of the last dimension of the results array
+  const auto& results_arr = tally->results();
+  const int results_dim = static_cast<int>(results_arr.shape()[2]);
+  std::array<int, 5> shape = {1, 1, 1, tally->n_scores(), results_dim};
 
   // set the shape for the filters applied on the tally
   for (int i = 0; i < tally->filters().size(); i++) {
@@ -586,7 +588,7 @@ void WeightWindows::update_weights(const Tally* tally, const std::string& value,
 
   // get a fully reshaped view of the tally according to tally ordering of
   // filters
-  auto tally_values = xt::reshape_view(tally->results(), shape);
+  auto tally_values = xt::reshape_view(results_arr, shape);
 
   // get a that is (particle, energy, mesh, scores, values)
   auto transposed_view = xt::transpose(tally_values, transpose);
diff --git a/tests/unit_tests/test_tallies.py b/tests/unit_tests/test_tallies.py
index c38f067d5..7b1bf0a2f 100644
--- a/tests/unit_tests/test_tallies.py
+++ b/tests/unit_tests/test_tallies.py
@@ -1,6 +1,8 @@
+from math import sqrt
 import numpy as np
 import pytest
 import openmc
+import scipy.stats as sps
 
 
 def test_xml_roundtrip(run_in_tmpdir):
@@ -163,3 +165,214 @@ def test_tally_application(sphere_model, run_in_tmpdir):
     assert (sp_tally.std_dev == tally.std_dev).all()
     assert (sp_tally.mean == tally.mean).all()
     assert sp_tally.nuclides == tally.nuclides
+
+def _tally_from_data(x, *, higher_moments=True, normality=True):
+    t = openmc.Tally()
+    t.scores = ["flux"]  # 1 score
+    t.nuclides = [openmc.Nuclide("H1")]  # 1 nuclide
+    t._sp_filename = "dummy.h5"  # mark "results available"
+    t._results_read = True  # don't try to read from disk
+    t._num_realizations = int(len(x))  # n
+    t.higher_moments = bool(higher_moments)
+
+    x = np.asarray(x, dtype=float)
+    # (num_filter_bins=1, num_nuclides=1, num_scores=1) -> (1,1,1) arrays
+    t._sum = np.array([[[np.sum(x)]]], dtype=float)
+    t._sum_sq = np.array([[[np.sum(x**2)]]], dtype=float)
+    if higher_moments:
+        t._sum_third = np.array([[[np.sum(x**3)]]], dtype=float)
+        t._sum_fourth = np.array([[[np.sum(x**4)]]], dtype=float)
+    return t
+
+@pytest.mark.parametrize(
+    "x, skew_true, kurt_true",
+    [   # Rademacher distribution
+        (np.array([1.0, -1.0] * 200), 0.0, 1.0),
+        # Two-point {0,3} with p(0)=3/4, p(3)=1/4
+        (np.concatenate([np.zeros(600), np.full(200, 3.0)]), 2.0 / sqrt(3.0), 7.0 / 3.0),
+        # Bernoulli distribution
+        (np.concatenate([np.ones(300), np.zeros(700)]), (1 - 2 * 0.3) / sqrt(0.3 * 0.7), (1 - 3 * 0.3 + 3 * 0.3**2) / (0.3 * 0.7)),
+    ],
+)
+def test_b1_b2_analytical_against_tally(x, skew_true, kurt_true):
+    t = _tally_from_data(x, higher_moments=True, normality=False)
+
+    g1 = t.skew(bias=True)[0, 0, 0]
+    b2 = t.kurtosis(bias=True, fisher=False)[0, 0, 0]
+
+    assert np.isclose(g1, skew_true, rtol=0, atol=1e-12)
+    assert np.isclose(b2, kurt_true, rtol=0, atol=1e-12)
+
+@pytest.mark.parametrize(
+    "draw, skew_true, kurt_true",
+    [(lambda rng, n: rng.normal(0, 1, n), 0.0, 3.0),  # Normal
+     (lambda rng, n: rng.random(n), 0.0, 1.8),  # Uniform(0,1)
+     (lambda rng, n: rng.exponential(1.0, n), 2.0, 9.0),  # Exp(1)
+     (lambda rng, n: (rng.random(n) < 0.3).astype(float),
+            (1 - 2 * 0.3) / sqrt(0.3 * 0.7),
+            (1 - 3 * 0.3 + 3 * 0.3**2) / (0.3 * 0.7),),],)
+
+def test_b1_b2_scipy_and_theory(draw, skew_true, kurt_true):
+    rng = np.random.default_rng(12345)
+    N = 200_000
+    x = draw(rng, N)
+
+    # Tally outputs
+    t = _tally_from_data(x, higher_moments=True, normality=False)
+    g1_t = t.skew(bias=True)[0, 0, 0]
+    b2_t = t.kurtosis(bias=True, fisher=False)[0, 0, 0]
+
+    # SciPy (population, bias=True to match population-moment style)
+    skew_sp = sps.skew(x, bias=True)
+    kurt_sp = sps.kurtosis(x, fisher=False, bias=True)
+
+    # Compare to SciPy numerically
+    assert np.isclose(g1_t, skew_sp, rtol=0, atol=5e-3)
+    assert np.isclose(b2_t, kurt_sp, rtol=0, atol=5e-3)
+
+    # Compare to analytical targets with size-dependent tolerances
+    tol_skew = 0.02 if abs(skew_true) < 0.5 else 0.05
+    tol_kurt = 0.03 if kurt_true < 4 else 0.1
+    assert abs(g1_t - skew_true) < tol_skew
+    assert abs(b2_t - kurt_true) < tol_kurt
+
+
+def test_kurtosis_bias_fisher_combinations():
+    """Test that all combinations of bias and fisher match scipy.stats.kurtosis"""
+    rng = np.random.default_rng(42)
+    x = rng.normal(0, 1, 10000)
+
+    t = _tally_from_data(x, higher_moments=True, normality=False)    # Test all four combinations
+    # 1. bias=True, fisher=False (Pearson's kurtosis, b2)
+    b2_tally = t.kurtosis(bias=True, fisher=False)[0, 0, 0]
+    b2_scipy = sps.kurtosis(x, fisher=False, bias=True)
+    assert np.isclose(b2_tally, b2_scipy, rtol=0, atol=1e-10)
+    assert np.isclose(b2_tally, 3.0, rtol=0.05, atol=0.1)  # Should be ~3 for normal
+
+    # 2. bias=True, fisher=True (excess kurtosis, g2)
+    g2_tally = t.kurtosis(bias=True, fisher=True)[0, 0, 0]
+    g2_scipy = sps.kurtosis(x, fisher=True, bias=True)
+    assert np.isclose(g2_tally, g2_scipy, rtol=0, atol=1e-10)
+    assert np.isclose(g2_tally, 0.0, rtol=0, atol=0.1)  # Should be ~0 for normal
+    assert np.isclose(g2_tally, b2_tally - 3.0, rtol=0, atol=1e-10)  # g2 = b2 - 3
+
+    # 3. bias=False, fisher=True (adjusted excess kurtosis, G2)
+    G2_tally = t.kurtosis(bias=False, fisher=True)[0, 0, 0]
+    G2_tally_default = t.kurtosis()[0, 0, 0]  # Should be same as default
+    G2_scipy = sps.kurtosis(x, fisher=True, bias=False)
+    assert np.isclose(G2_tally, G2_tally_default, rtol=0, atol=1e-10)
+    assert np.isclose(G2_tally, G2_scipy, rtol=0, atol=1e-10)
+    assert np.isclose(G2_tally, 0.0, rtol=0, atol=0.1)  # Should be ~0 for normal
+
+    # 4. bias=False, fisher=False (adjusted Pearson's kurtosis)
+    adj_b2_tally = t.kurtosis(bias=False, fisher=False)[0, 0, 0]
+    adj_b2_scipy = sps.kurtosis(x, fisher=False, bias=False)
+    assert np.isclose(adj_b2_tally, adj_b2_scipy, rtol=0, atol=1e-10)
+    assert np.isclose(adj_b2_tally, 3.0, rtol=0.05, atol=0.1)  # Should be ~3 for normal
+    assert np.isclose(adj_b2_tally, G2_tally + 3.0, rtol=0, atol=1e-10)  # adj_b2 = G2 + 3
+
+
+def test_ztests_scipy_comparison():
+    rng = np.random.default_rng(987)
+    x_norm = rng.normal(size=50_000)
+    x_exp = rng.exponential(size=50_000)
+
+    # -------- Normal dataset (should not reject) --------
+    t0 = _tally_from_data(x_norm, higher_moments=True, normality=True)
+    Zb1_0, p_skew_0 = t0.skewtest(alternative="two-sided")
+    Zb2_0, p_kurt_0 = t0.kurtosistest(alternative="two-sided")
+    K2_0, p_omni_0 = t0.normaltest(alternative="two-sided")
+
+    Zb1_0 = Zb1_0.ravel()[0]
+    p_skew_0 = p_skew_0.ravel()[0]
+    Zb2_0 = Zb2_0.ravel()[0]
+    p_kurt_0 = p_kurt_0.ravel()[0]
+    K2_0 = K2_0.ravel()[0]
+    p_omni_0 = p_omni_0.ravel()[0]
+
+    z_skew_sp0, p_skew_sp0 = sps.skewtest(x_norm)
+    z_kurt_sp0, p_kurt_sp0 = sps.kurtosistest(x_norm)
+    k2_sp0, p_omni_sp0 = sps.normaltest(x_norm)
+
+    assert np.isclose(Zb1_0, z_skew_sp0, atol=0.15)
+    assert np.isclose(Zb2_0, z_kurt_sp0, atol=0.15)
+    assert np.isclose(K2_0, k2_sp0, atol=0.30)
+    assert np.isclose(p_skew_0, p_skew_sp0, atol=5e-3)
+    assert np.isclose(p_kurt_0, p_kurt_sp0, atol=5e-3)
+    assert np.isclose(p_omni_0, p_omni_sp0, atol=5e-3)
+
+    # -------- Exponential dataset (should strongly reject) --------
+    t1 = _tally_from_data(x_exp, higher_moments=True, normality=True)
+
+    Zb1_1, p_skew_1 = t1.skewtest(alternative="two-sided")
+    Zb2_1, p_kurt_1 = t1.kurtosistest(alternative="two-sided")
+    K2_1, p_omni_1 = t1.normaltest(alternative="two-sided")
+
+    Zb1_1 = Zb1_1.ravel()[0]
+    p_skew_1 = p_skew_1.ravel()[0]
+    Zb2_1 = Zb2_1.ravel()[0]
+    p_kurt_1 = p_kurt_1.ravel()[0]
+    K2_1 = K2_1.ravel()[0]
+    p_omni_1 = p_omni_1.ravel()[0]
+
+    z_skew_sp1, p_skew_sp1 = sps.skewtest(x_exp)
+    z_kurt_sp1, p_kurt_sp1 = sps.kurtosistest(x_exp)
+    k2_sp1, p_omni_sp1 = sps.normaltest(x_exp)
+
+    # Both pipelines should reject very strongly
+    assert p_skew_1 < 1e-6 and p_skew_sp1 < 1e-6
+    assert p_kurt_1 < 1e-6 and p_kurt_sp1 < 1e-6
+    assert p_omni_1 < 1e-6 and p_omni_sp1 < 1e-6
+
+    # Right-skewed and heavy-tailed → large positive Z-statistics
+    assert Zb1_1 > 30 and z_skew_sp1 > 30
+    assert Zb2_1 > 30 and z_kurt_sp1 > 30
+    assert K2_1 > 2000 and k2_sp1 > 2000
+
+def test_vov_stochastic(sphere_model, run_in_tmpdir):
+    tally = openmc.Tally(name="test tally")
+    ef = openmc.EnergyFilter([0.0, 0.1, 1.0, 10.0e6])
+    mesh = openmc.RegularMesh.from_domain(sphere_model.geometry, (2, 2, 2))
+    mf = openmc.MeshFilter(mesh)
+    tally.filters = [ef, mf]
+    tally.scores = ["flux", "absorption", "fission", "scatter"]
+    tally.higher_moments = True
+    sphere_model.tallies = [tally]
+
+    sp_file = sphere_model.run(apply_tally_results=True)
+
+    assert tally._mean is None
+    assert tally._std_dev is None
+    assert tally._sum is None
+    assert tally._sum_sq is None
+    assert tally._sum_third is None
+    assert tally._sum_fourth is None
+    assert tally._num_realizations == 0
+    assert tally._sp_filename == sp_file
+
+    with openmc.StatePoint(sp_file) as sp:
+        assert tally in sp.tallies.values()
+        sp_tally = sp.tallies[tally.id]
+
+    assert np.all(sp_tally.std_dev == tally.std_dev)
+    assert np.all(sp_tally.mean == tally.mean)
+    assert np.all(sp_tally.vov == tally.vov)
+    assert sp_tally.nuclides == tally.nuclides
+
+    n = sp_tally.num_realizations
+    mean = sp_tally.mean
+    sum_ = sp_tally._sum
+    sum_sq = sp_tally._sum_sq
+    sum_third = sp_tally._sum_third
+    sum_fourth = sp_tally._sum_fourth
+
+    expected_vov = np.zeros_like(mean)
+    nonzero = np.abs(mean) > 0
+
+    num = (sum_fourth - (4.0*sum_third*sum_)/n + (6.0*sum_sq*sum_**2)/(n**2)
+           - (3.0*sum_**4)/(n**3))
+    den = (sum_sq - (1.0/n)*sum_**2)**2
+
+    expected_vov[nonzero] = num[nonzero]/den[nonzero] - 1.0/n
+
+    assert np.allclose(expected_vov, sp_tally.vov, rtol=1e-7, atol=0.0)