© Y. Jo and D. Lee, Published by EDP Sciences, 2026

1. Introduction

The Monte Carlo method [1] is widely used for nuclear reactor analysis with the increase of computing power. It simulates particle histories by using pseudo-random numbers to achieve statistically converged results. It can use high-fidelity physics phenomena through stochastic particle simulations. In the Monte Carlo method, the exact model of problem geometry can be used with continuous energy cross section data. The reactor analysis requires computing the effect of small changes in input parameter. The number of perturbed reactor states can exceed a thousand, depending on the range of interests such as fuel temperature, moderator temperature, xenon concentration and boron concentration. Unfortunately, the computational burden of Monte Carlo method is too huge to run every perturbed reactor state. Additionally, the Monte Carlo method struggles to compute the effect of small perturbations due to its statistical uncertainty.

The Monte Carlo perturbation theory is used to solve this problem. By using perturbation methods, which are the differential operator sampling (DOS) and correlated sampling (CS) methods [2], it is possible to estimate the k-eigenvalue change in a perturbed system during a forward Monte Carlo calculation of an unperturbed system. However, it was observed that the neglection of perturbed fission source effect may cause a large error [2]. To improve its accuracy, Nagaya and Mori presented perturbation methods calculating the perturbed source effect explicitly [3, 4]. Additionally, the first-order adjoint-weighted perturbation (AWP) method was developed [5–9]. The AWP method uses the adjoint flux of unperturbed system as a weight function to consider a fission source perturbation effect implicitly. However, those conventional first-order perturbation methods still have a limitation of first-order accuracy. Therefore, the error of k-eigenvalue change estimated by the conventional first-order perturbation method becomes large for the perturbed system where product of perturbations are significant. It is because the product of perturbations is assumed to be negligible during a derivation of the conventional first-order perturbation methods. Examples showing large errors due to large perturbations are presented in Section 3. To overcome this problem, Truchet et al. [10, 11], and Tuya and Nagaya [12] presented perturbation methods based on the exact perturbation theory (EPT). However, their approach requires additional particle simulations in the perturbed system for L + 1 generations, where L is the number of latent generations, to compute the adjoint flux of the perturbed system. The computational time increases by a factor of L + 2 due to the computational burden of additional particle simulations.

In this work, a calculation technique named as a perturbation-included iterated fission probability (PIFP) method is developed to overcome the limitation of first-order accuracy in conventional perturbation methods and the computational burden of conventional exact perturbation (EP) methods. The PIFP method computes the adjoint flux of perturbed system for the EPT in Monte Carlo transport. The EP method is implemented in the Monte Carlo code MCS. In Section 2, the conventional perturbation theory and the PIFP method are presented. In Section 3, the results of perturbation methods for k-eigenvalue change are presented for the Godiva benchmark and PLUS7 fuel assembly.

2. Material and methods

2.1. Perturbation theory

In the perturbation theory, two reactor systems are considered. One is the unperturbed system (system 1), and the other is the perturbed system (system 2) where there is a change in input parameters. The steady-state Boltzmann transport equations for the unperturbed and perturbed systems are expressed as equations (1) and (2), respectively, with a subscript representing the system index.

$$\begin{array}{c}\hfill \left(\mathbf{I}-{\lambda }_1{\mathbf{H}}_1\right)S_1=0,\end{array}$$(1)

$$\begin{array}{c}\hfill \left(\mathbf{I}-{\lambda }_2{\mathbf{H}}_2\right)S_2=0,\end{array}$$(2)

where S is the fission source density, H is the fission operator, and λ is the eigenvalue. Equation (2) can be re-written as follows by using Δλ = λ₂ − λ₁, ΔH = H₂ − H₁, and ΔS = S₂ − S₁.

$$\begin{array}{c}\hfill \left(\mathbf{I}-\left({\lambda }_1+\mathrm{\Delta }\lambda \right)\left({\mathbf{H}}_1+\mathrm{\Delta }\mathbf{H}\right)\right)\left(S_1+\mathrm{\Delta }S\right)=0.\end{array}$$(3)

One can obtain the first-order perturbation equation by taking inner product of Equation (3) with the weight function τ which is a function of ones and neglecting terms having product of perturbations through first-order approximation. Additionally, by assuming the fission source density change ΔS is negligible, the first-order perturbation equation is expressed as follows:

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }\mathrm{k}}{k_2}=-\frac{\mathrm{\Delta }\lambda }{{\lambda }_1}\cong \frac{\langle \tau ,{\lambda }_1\mathrm{\Delta }\mathbf{H}{S}_1\rangle }{\langle \tau ,{\lambda }_1{\mathbf{H}}_1{S}_1\rangle }.\end{array}$$(4)

By taking inner product of Equation (3) with the adjoint flux of unperturbed system ${\varphi }_1^{*}$ and neglecting terms having product of perturbations through first-order approximation, one can obtain the first-order AWP equation for k-eigenvalue as follows:

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }\mathrm{k}}{k_2}\cong \frac{\langle {\varphi }_1^{\ast },{\lambda }_1\mathrm{\Delta }\mathbf{H}{S}_1\rangle }{\langle {\varphi }_1^{\ast },{\lambda }_1\left({\mathbf{H}}_1+\mathrm{\Delta }\mathbf{H}\right)S_1\rangle },\end{array}$$(5)

where the bracket denotes an integration over the phase space, which are position, direction, and energy. The first-order AWP equation also can be expressed as follows [5]:

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }\mathrm{k}}{k_2}\cong -\frac{\langle {\varphi }_1^{\ast },\left(\mathrm{\Delta }\mathbf{L}-{\lambda }_1\mathrm{\Delta }\mathbf{F}\right){\varphi }_1\rangle }{\langle {\varphi }_1^{\ast },{\lambda }_1\left({\mathbf{F}}_1+\mathrm{\Delta }\mathbf{F}\right){\varphi }_1\rangle },\end{array}$$(6)

where L and F are the net loss and net fission production operators, respectively. ϕ is the neutron flux. The first-order AWP method has a limitation in its accuracy due to the first-order approximation used to neglect terms having product of perturbations. On the other hand, the EPT allows to estimate the perturbed k-eigenvalue exactly without applying any approximations. The EPT equation for k-eigenvalue is expressed as follows by taking inner product of Equation (3) with the adjoint flux of perturbed system ${\varphi }_2^{*}$.

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }\mathrm{k}}{k_2}=\frac{\langle {\varphi }_2^{\ast },{\lambda }_1\mathrm{\Delta }\mathbf{H}{S}_1\rangle }{\langle {\varphi }_2^{\ast },{\lambda }_1\left({\mathbf{H}}_1+\mathrm{\Delta }\mathbf{H}\right)S_1\rangle }.\end{array}$$(7)

The unknown term on the right-hand side of Equation (7) is the adjoint flux of perturbed system ${\varphi }_2^{*}$ which can be calculated by the PIFP method. The PIFP method is introduced in the next section.

2.2. Perturbation-included iterated fission probability method

The PIFP method is based on the iterated fission probability (IFP) method [13]. By using CS and IFP methods together, the adjoint flux of perturbed system ${\varphi }_2^{*}$ is calculated during a forward Monte Carlo calculation of unperturbed system.

In the PIFP method, no additional particle simulation is required to compute the adjoint flux of perturbed system. In this section, the algorithm of PIFP method is presented. The HS term in Equation (1) can be expressed as follows in terms of kernels:

$$\begin{array}{cc}& \int d{\mathbf{P}}_0\mathbf{H}\left({\mathbf{P}}_0\to \mathbf{P}\right)S\left({\mathbf{P}}_0\right)\hfill \\ \hfill & =\sum _{n=0}^{\infty }\int d{{\mathbf{P}}^{\prime }}_n\cdots \int d{\mathbf{P}}_0{C}_f\left({{\mathbf{P}}^{\prime }}_n\to \mathbf{P}\right)T_{n+1}\left({\mathbf{P}}_n\to {{\mathbf{P}}^{\prime }}_n\right)\hfill \\ \hfill & ·C_{s,n}\left({{\mathbf{P}}^{\prime }}_{n-1}\to {\mathbf{P}}_n\right)T_n\left({\mathbf{P}}_{n-1}\to {{\mathbf{P}}^{\prime }}_{n-1}\right)\hfill \\ \hfill & \cdots C_{s,1}\left({{\mathbf{P}}^{\prime }}_0\to {\mathbf{P}}_1\right)T_1\left({\mathbf{P}}_0\to {{\mathbf{P}}^{\prime }}_0\right)S\left({\mathbf{P}}_0\right),\hfill \end{array}$$(8)

where C_f is the fission collision kernel, T_i is the i-th transport kernel, C_s, i is the i-th scattering collision kernel, and n is the number of scattering collisions before fission collision. In the Monte Carlo calculation, each kernel of Equation (8) is defined as follows:

$$\begin{array}{cc}\hfill T_i\left({\mathbf{P}}_{i-1}\to {{\mathbf{P}}^{\prime }}_{i-1}\right)& =T_i\left(E_{i-1},{\mathrm{\Omega }}_{i-1};{\mathbf{r}}_{i-1}\to {\mathbf{r}}_i\right)\hfill \\ \hfill & =exp\left(-{\mathrm{\Sigma }}_t\left({\mathbf{r}}_{i-1},E_{i-1}\right){\ell }_i\right),\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill C_{s,i}\left({{\mathbf{P}}^{\prime }}_{i-1}\to {\mathbf{P}}_i\right)& =C_{s,i}\left({\mathbf{r}}_i;E_{i-1},{\mathrm{\Omega }}_{i-1}\to E_i,{\mathrm{\Omega }}_i\right)\hfill \\ \hfill & ={\mathrm{\Sigma }}_s\left({\mathbf{r}}_i,E_{i-1}\right)f_s\left(E_{i-1},{\mathrm{\Omega }}_{i-1}\to E_i,{\mathrm{\Omega }}_i\right),\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill C_f\left({{\mathbf{P}}^{\prime }}_n\to \mathbf{P}\right)& =C_f\left({\mathbf{r}}_{n+1};E_n,{\mathrm{\Omega }}_n\to E,\mathrm{\Omega }\right)\hfill \\ \hfill & =\frac{\chi \left(E_n\to E\right)}{4\pi }\nu \left({\mathbf{r}}_{n+1},E_n\right){\mathrm{\Sigma }}_f\left({\mathbf{r}}_{n+1},E_n\right),\hfill \end{array}$$(11)

where Σ_t, Σ_s and Σ_f are the total, scattering and fission cross sections, respectively. l_i is the i-th neutron flight distance, f_s(E_i − 1, Ω_i − 1 → E_i, Ω_i) is the scattering probability that a neutron with incident energy E_i − 1 and incident direction Ω_i − 1 is scattered to the outgoing energy E_i and outgoing direction Ω_i. χ is the fission neutron energy spectrum, and ν is the average number of fission neutrons produced from a fission reaction. In the current implementation, the perturbation effect of scattering probability f_s is assumed to be neglected.

The ΔHS term can be calculated by using one of DOS and CS methods. The DOS method computes the target response change by using Taylor series. By applying the first-order Taylor series of ΔHS for Δx change in the input parameter x and simplifying phase space notations, the ΔHS term is expressed as follows:

$$\begin{array}{cc}\hfill \mathrm{\Delta }\mathbf{H}S& =\sum _{n=0}^{\infty }\int d{{\mathbf{P}}^{\prime }}_n\cdots \int d{\mathbf{P}}_0\mathrm{\Delta }x[\frac{1}{C_f}\frac{\partial C_f}{\partial x}\hfill \\ \hfill & +\frac{1}{T_{n+1}}\frac{\partial T_{n+1}}{\partial x}+\sum _{i=1}^n\left(\frac{1}{C_{s,i}}\frac{\partial C_{s,i}}{\partial x}+\frac{1}{T_i}\frac{\partial T_i}{\partial x}\right)]\hfill \\ \hfill & ·C_f{T}_{n+1}\prod _{i=1}^n\left(C_{s,i}{T}_i\right)S.\hfill \end{array}$$(12)

Equation (12) represents that a derivative of fission operator H can be estimated by scoring derivatives of each kernel, which are expressions in square brackets, during a forward Monte Carlo calculation of the unperturbed system. In the CS method, it is assumed that neutron histories in the perturbed system follows the neutron tracks in the unperturbed system. ΔHS is computed by scoring the ratio of perturbed and unperturbed kernels as follows:

$$\begin{array}{cc}\hfill \mathrm{\Delta }\mathbf{H}S& =\left(\frac{{\mathbf{H}}_2}{{\mathbf{H}}_1}-1\right){\mathbf{H}}_{1}S\hfill \\ \hfill & =\sum _{n=0}^{\infty }\int d{{\mathbf{P}}^{\prime }}_n\cdots \int d{\mathbf{P}}_0\left[\frac{C_{f,2}{T}_{n+1,2}\prod _{i=1}^n\left(C_{s,i,2}{T}_{i,2}\right)}{C_{f,1}{T}_{n+1,1}\prod _{i=1}^n\left(C_{s,i,1}{T}_{i,1}\right)}-1\right]\hfill \\ \hfill & ·C_{f,1}{T}_{n+1,1}\prod _{i=1}^n\left(C_{s,i,1}{T}_{i,1}\right)S.\hfill \end{array}$$(13)

The numerator and denominator in square brackets consist of perturbed kernels and unperturbed kernels, respectively. The expression in square brackets are used as a weight adjusting factor to represent particle weights of the perturbed system during a forward Monte Carlo calculation of unperturbed system. The CS method has the advantage that no approximation is applied unlike the DOS method using n-th order Taylor series. However, as it was mentioned in reference [2], the CS method should be carefully used since it may cause a large or unbounded variance of perturbed k-eigenvalue depending on the amount of perturbation. The result showing the large variance is presented in Section 3.

The adjoint flux of unperturbed system can be estimated as follows by using the power method and IFP interpretation [5, 7].

$$\begin{array}{c}\hfill {\varphi }_{1,n}^{\ast }\left({\mathbf{P}}_0\right)=C\int d{\mathbf{P}}_n\frac{1}{k_1^n}{H}_1^n\left({\mathbf{P}}_0\to {\mathbf{P}}_n\right){\varphi }_0^{\ast }\left({\mathbf{P}}_0\right),\end{array}$$(14)

where n is treated as the number of latent generations in the IFP method, and C is an arbitrary constant. ${\varphi }_{1,n}^{*}$ is the n-th iterative solution of adjoint flux in the unperturbed system from the initial guess ${\varphi }_0^{*}$. As the number of latent generations increases, the iterative solution converges to the true solution of adjoint flux ${\varphi }_1^{*}$. By assuming ${\varphi }_0^{*}({\mathbf{P}}_{\mathbf{0}})=S({\mathbf{P}}_{\mathbf{0}})$ in Equation (14) as a reference [7], it can be expressed as

$$\begin{array}{cc}\hfill {\varphi }_{1,n}^{\ast }\left({\mathbf{P}}_0\right)& =C\int d{\mathbf{P}}_n\cdots \int d{\mathbf{P}}_1\prod _{i=0}^{n-1}\hfill \\ \hfill & \times \left[\frac{1}{k_1^{\left(i\right)}}{H}_1^{\left(i\right)}\left({\mathbf{P}}_i\to {\mathbf{P}}_{i+1}\right)\right]S\left({\mathbf{P}}_0\right),\hfill \end{array}$$(15)

where $k_1^{(i)}$ and $H_1^{\left(i\right)}$ are the k-eigenvalue and fission operator of the unperturbed system in i-th generation from the fission source density at phase space P₀, respectively. Equation (15) represents the adjoint flux at phase space P₀ is proportional to the number of progeny fission neutrons generated in the n-th generation from a unit source at P₀. The adjoint flux at phase space P₀ in the perturbed system ${\varphi }_{2,n}^{*}({\mathbf{P}}_{\mathbf{0}})$ is expressed as

$$\begin{array}{cc}\hfill {\varphi }_{2,n}^{\ast }\left({\mathbf{P}}_0\right)& =C^{\prime }\int d{\mathbf{P}}_n\cdots \int d{\mathbf{P}}_1\prod _{i=0}^{n-1}\hfill \\ \hfill & \times \left[\frac{1}{k_2^{\left(i\right)}}{H}_2^{\left(i\right)}\left({\mathbf{P}}_i\to {\mathbf{P}}_{i+1}\right)\right]S\left({\mathbf{P}}_0\right).\hfill \end{array}$$(16)

Equation (16) also can be expressed as

$$\begin{array}{cc}& {\varphi }_{2,n}^{\ast }\left({\mathbf{P}}_0\right)=C^{″}\int d{\mathbf{P}}_n\cdots \int d{\mathbf{P}}_1\prod _{i=0}^{n-1}\hfill \\ \hfill & \times \left[\frac{1}{k_1^{\left(i\right)}}\left(H_1^{\left(i\right)}\left({\mathbf{P}}_i\to {\mathbf{P}}_{i+1}\right)+\mathrm{\Delta }{H}^{\left(i\right)}\left({\mathbf{P}}_i\to {\mathbf{P}}_{i+1}\right)\right)\right]S\left({\mathbf{P}}_0\right),\hfill \end{array}$$(17)

where $C″=C\prime {\prod }_{i=0}^{n-1}(k_1^{\left(i\right)}/k_2^{\left(i\right)})$. Equation (17) represents the adjoint flux of the perturbed system can be calculated during a forward Monte Carlo calculation of the unperturbed system by performing a cumulative product of perturbed fission operator (H₁ + ΔH) during latent generations. The arbitrary constant C″ is unknown because it contains the k-eigenvalue of the perturbed system. However, it can be neglected since it will be cancelled out from the adjoint flux multiplied both in the numerator and denominator of Equation (7). By substituting equations (13) and (17) into Equation (7), one can obtain the estimation of k-eigenvalue change for the EP method as follows:

$$\begin{array}{cc}& {\left(\frac{\mathrm{\Delta }\mathrm{k}}{k_2}\right)}^{\left(i+n\right)}=\hfill \\ \hfill & \frac{\sum _j\left\{\prod _{l=0}^{n-1}\left[\frac{1}{k_1^{\left(i+1+l\right)}}\left({\mathbf{H}}_1^{\left(i+1+l\right)}+\mathrm{\Delta }{\mathbf{H}}^{\left(i+1+l\right)}\right)\right]·\frac{1}{k_1^{\left(i\right)}}\mathrm{\Delta }{\mathbf{H}}^{\left(i\right)}{S}_j^{\left(i\right)}\right\}}{\sum _j\left\{\prod _{l=0}^{n-1}\left[\frac{1}{k_1^{\left(i+1+l\right)}}\left({\mathbf{H}}_1^{\left(i+1+l\right)}+\mathrm{\Delta }{\mathbf{H}}^{\left(i+1+l\right)}\right)\right]·\frac{1}{k_1^{\left(i\right)}}\left({\mathbf{H}}_1^{\left(i\right)}+\mathrm{\Delta }{\mathbf{H}}^{\left(i\right)}\right)S_j^{\left(i\right)}\right\}},\hfill \end{array}$$(18)

where the superscript represents the cycle index, j is a particle index, and l is the number of latent generations in IFP method. ${\left(\mathrm{\Delta }k/k_2\right)}^{(i+n)}$ is the estimation of k-eigenvalue change calculated at cycle (i + n). $k_1^{(i+1+l)}$ is the k-eigenvalue of unperturbed system calculated at cycle (i + 1 + l). ${\mathbf{H}}_{\mathbf{1}}^{(\mathbf{i})}$ and ΔH ⁽ⁱ⁾ are the fission operator of the unperturbed system and fission operator change at cycle (i). $S_j^{\left(i\right)}$ is the j-th particle in cycle (i).

Even though the large or unbounded variance problem is solved by using the DOS method instead of CS method for PIFP, it is determined to use the CS method for the PIFP method to compute the perturbation term ΔHS. It is because no n-th order approximation is applied in the CS method unlike the DOS method using n-th order Taylor series.

2.3. Monte Carlo code MCS

A Monte Carlo code MCS has been developed at Ulsan National Institute of Science and Technology (UNIST) since 2013 [14]. Both criticality and fixed-source modes are available in MCS for reactor analysis and shielding calculation. The criticality and multi-cycle depletion calculation capability of MCS has been verified and validated by solving various benchmark problems including the BEAVRS benchmark [15], and about 300 cases from the International Criticality Safety Benchmark Experimental Problems (ICSBEP) [16, 17].

3. Perturbation results for k-eigenvalue

This section presents perturbation results to estimate the k-eigenvalue change for various types of perturbation problems. Section 3.1 presents the perturbation results of 2-group infinite homogeneous problems. For Sections 3.2 and 3.3, results of perturbation methods implemented in the Monte Carlo code MCS are compared for Godiva benchmark and PLUS7 fuel assembly model. The continuous energy ENDF/B-VII.1 nuclear data library is used for Monte Carlo calculations. The reference result, which is the direct perturbation (DP) result, is calculated by running two independent inputs for unperturbed and perturbed systems. DOS, AWP, and exact perturbation using PIFP method (EPP), are compared against the reference DP results.

3.1. 2-group infinite homogeneous problem

In this section, the k-eigenvalue change due to the cross-section change is analytically calculated by solving 2-group infinite homogeneous problem. Table 1 shows the 2-group cross section data of unperturbed system. Σ_a is the absorption cross section. Σ_{s, g → g′} and χ_g→g′ are the scattering cross section and the fission neutron spectrum from group g to group g′, respectively. The analytic solution of k-eigenvalue for the 2-group infinite homogeneous problem is expressed as Equation (19). The analytic solution of unperturbed system having the cross section data in Table 1 is 1.09628.

$$\begin{array}{c}\hfill k=\frac{-b+\sqrt{b^2-4ac}}{2a},\end{array}$$(19)

where

$$\begin{array}{c}\hfill \begin{array}{c}a=\left({\mathrm{\Sigma }}_{\mathrm{a}1}+{\mathrm{\Sigma }}_{\mathrm{f}1}+{\mathrm{\Sigma }}_{\mathrm{s}12}\right)\left({\mathrm{\Sigma }}_{\mathrm{a}2}+{\mathrm{\Sigma }}_{\mathrm{f}2}+{\mathrm{\Sigma }}_{\mathrm{s}21}\right)-{\mathrm{\Sigma }}_{s21}{\mathrm{\Sigma }}_{s12},\hfill \\ b=-[\left({\mathrm{\Sigma }}_{\mathrm{a}1}+{\mathrm{\Sigma }}_{\mathrm{f}1}+{\mathrm{\Sigma }}_{\mathrm{s}12}\right){\chi }_{22}{\nu }_2{\mathrm{\Sigma }}_{f2}\hfill \\ +\left({\mathrm{\Sigma }}_{\mathrm{a}2}+{\mathrm{\Sigma }}_{\mathrm{f}2}+{\mathrm{\Sigma }}_{\mathrm{s}21}\right){\chi }_{11}{\nu }_1{\mathrm{\Sigma }}_{f1}+{\mathrm{\Sigma }}_{s21}{\chi }_{12}{\nu }_1{\mathrm{\Sigma }}_{f1}\hfill \\ +{\mathrm{\Sigma }}_{s12}{\chi }_{21}{\nu }_2{\mathrm{\Sigma }}_{f2}],\hfill \\ c={\chi }_{11}{\nu }_1{\mathrm{\Sigma }}_{f1}{\chi }_{22}{\nu }_2{\mathrm{\Sigma }}_{f2}-{\chi }_{21}{\nu }_2{\mathrm{\Sigma }}_{f2}{\chi }_{12}{\nu }_1{\mathrm{\Sigma }}_{f1}.\hfill \end{array}\end{array}$$

Figures 1 and 2 show the analytic results of perturbation methods for Σ_a1 change and Σ_{s, 1 → 2} change, respectively. In Figure 1, maximum errors of DOS and AWP methods are −5876 pcm and −580 pcm, respectively. In Figure 2, maximum errors of DOS and AWP methods are 14122 pcm and 6281 pcm, respectively. The AWP method shows more accurate results compared with the DOS method. It is because the adjoint flux of unperturbed system is used as the weight function in AWP method. The EP method estimates the perturbed k-eigenvalue exactly since it uses the EPT.

Fig. 1. Analytic results of perturbation methods for Σa1 change in the 2-group infinite homogeneous problem.

Fig. 2. Analytic results of perturbation methods for Σs, 1 → 2 change in the 2-group infinite homogeneous problem.

3.2. Godiva benchmark

Godiva is a highly enriched uranium sphere of 8.7407 cm radius. The density of uranium is 18.74 g/cm³, and it consists of 1.02 wt.% U-234, 93.71 wt.% U-235, and 5.27 wt.% U-238. For the Godiva benchmark, three density perturbation cases are considered as follows: 1) density perturbation in the outermost shell region, 2) density perturbation in the center sphere region, and 3) density perturbation in the global region. The Godiva sphere is divided into 10 equi-volume regions. In the outermost shell region case, the density perturbation is applied only to the outermost shell region of 0.3016 cm thickness. In the center sphere region case, the density perturbation is applied to the center sphere of 4.0571 cm radius. In the global region case, the density perturbation is applied to the whole region of Godiva. The reference calculations for DP are performed by using 40 active cycles with 10 million histories per a cycle. The calculations for perturbation methods are performed by using 40 active cycles with 40 million histories per a cycle. The k-eigenvalue of base case is 0.99980 ± 0.00002.

Fig. 3. Results of perturbation methods for local density change in the shell region of Godiva: k-eigenvalue change (up) and error (down).

The local density change ranging 0.2–1.8 of nominal density is considered in the outermost shell region. Figure 3 shows the results of perturbation methods for local density change in the outermost shell region of Godiva. The error bars represent 1-σ statistical errors. The k-eigenvalue changes −2400–2800 pcm depending on the local density change. The maximum error of DOS method is larger than 2000 pcm due to the neglection of the perturbed source effect. The AWP(DOS) and AWP(CS) methods shows maximum errors less than 200 pcm. On the other hand, the EP method shows the most accurate results with errors less than 10 pcm. This is because EPP method is based on the EPT.

The local density change ranging 0.2–1.8 of nominal density is considered in the center sphere region. Figure 4 shows the results of perturbation methods for local density change in the center sphere region of Godiva. As it was observed in the outermost shell region case, large errors are observed in the DOS method due to the neglection of the perturbed source effect. The maximum errors of AWP(DOS) and AWP(CS) methods are −4083 pcm and −6219 pcm, respectively. Those large errors of AWP cases are due to the neglection of terms having product of perturbations. On the other hand, the maximum error of EP method is 397 pcm. For the density change ranging 0.2–1.5 of nominal density, the maximum error of EPP method is 13 pcm. For cases having a density perturbation larger than 0.5 of nominal density, the statistical uncertainty of EPP method becomes large or unbounded. As it was mentioned in reference [2], this is due to the limitation of CS method using the same neutron track histories for the estimation of the perturbed system. In this work, no reset condition is applied for the CS method to prevent the large variance.

Fig. 4. Results of perturbation methods for local density change in center sphere region of Godiva: k-eigenvalue change (up) and error (down).

Fig. 5. Results of perturbation methods for global density change in Godiva: k-eigenvalue change (up) and error (down).

Fig. 6. Results of perturbation methods for boron concentration change in PLUS7 fuel assembly: k-eigenvalue change (up) and error (down).

The global density change ranging 0.2–1.8 of nominal density is considered in the global region. Figure 5 shows the results of perturbation methods for global density change in Godiva. The maximum errors of DOS, AWP(DOS), and AWP(CS) methods are 14981, 7485 and −12002 pcm, respectively. The maximum error of EPP method is −1203 pcm for the case having 1.8 of nominal density. For the density change ranging 0.5–1.4 of nominal density, the maximum error of EPP method is 11 pcm. As it was observed in Figure 4, the statistical uncertainty of EPP method becomes large or unbounded for the cases having a density perturbation larger than 0.5 of nominal density. The error of EPP method for those cases also increases.

3.3. PLUS7 fuel assembly

A 2-D 16 × 16 PLUS7 fuel assembly model is used to check the effectiveness of EP method for a practical problem. A 16 × 16 PLUS7 fuel assembly contains 236 fuel rods and 5 large guide tubes. The fuel enrichments of normal and zoned fuel rods are 4.50 wt.% and 4.00 wt.%, respectively. The geometry specification is referred from the APR1400 benchmark [18]. The calculations are performed by using 25 active cycles with 8 million histories per a cycle. The k-eigenvalue of base case is 1.40896 ± 0.00005.

Fig. 7. PLUS7 fuel assembly models having different number of gadolinia rods.

Fig. 8. Results of perturbation methods for the number of gadolinia rods change in PLUS7 fuel assembly.

The soluble boron concentration change ranging 0–5000 ppm is considered in PLUS7 fuel assembly. The boron concentration is 0 ppm in the unperturbed model. Figure 6 shows the results of perturbation methods for soluble boron concentration change in PLUS7 fuel assembly. Maximum errors of DOS, AWP(DOS) and AWP(CS) methods are −22015, −2246, and 9481 pcm, respectively. The EPP method results show a significantly improved accuracy by means of EPT. Even when the boron concentration is changed from 0 ppm to 5000 ppm, the error of EPP method is 6 pcm. The maximum error of EPP method is 14 pcm. The computational time of perturbation methods are also compared for the soluble boron concentration change in PLUS7 fuel assembly. As shown in Figure 6, the 28 perturbation cases changing soluble boron concentration from 0 ppm are considered. 155 processors of Intel Xeon Gold 6242R are used for the MPI parallel calculations. The DP case requires 29 Monte Carlo runs (one base case and 28 perturbation cases), and it takes 7.54 hours. When a perturbation method is used, one can obtain k-eigenvalue changes of 28 perturbation cases just by running a single unperturbed input. The computational time of AWP(CS) and EPP method is 0.50 hour which is 15.1 times faster than the DP case.

The k-eigenvalue change from the number of burnable poison rods change is calculated by perturbation methods. Figure 7 shows PLUS7 fuel assembly models having different numbers of gadolinia rods. The fuel assembly containing no gadolinia rod, shown in Figure 7a, is set as a base fuel assembly model. The k-eigenvalue of the base model is 1.40889 ± 0.00005. Figure 8 shows the results of perturbation methods for the number of gadolinia rods change in the PLUS7 fuel assembly. The case name has a format of “[Gd₂O₃ enrichment (wt.%)]wt.#[number of gadolinia rods]”. For example, 6 wt.#08 case represents the fuel assembly having 8 gadolinia rods with 6 wt.% of Gd₂O₃ enrichment. The error bars represent 1-σ statistical errors in Monte Carlo simulations. DOS results show a negative k-eigenvalue due to the first-order approximation and neglected fission source perturbation during the derivation of equation (4). Maximum errors of AWP(DOS) and AWP(CS) methods are −94176 pcm and 4846 pcm, respectively, due to their first-order approximation. In linear perturbation methods, higher-order perturbation terms are assumed to be negligible. Therefore, the error of AWP methods become large for the perturbation cases where the higher-order perturbation terms are significant. On the other hand, the EPP method shows significantly improved accuracy with a maximum error of 23 pcm. This is because the adjoint flux of perturbed system is used in the EPP method.

The k-eigenvalue changes due to the change in multiple parameters are calculated for PLUS7 fuel assembly. Table 2 shows perturbation parameters for PLUS7 fuel assembly. By selecting one of the perturbed values for each parameter in Table 2, 48 fuel assembly cases, which are one base case and 47 perturbed cases, are generated. The 48 fuel assembly cases are shown in Table 3. Additionally, four reactivity coefficients, which are fuel temperature coefficients (FTC), xenon worth, moderator density coefficients (MDC), and boron worth, are also calculated for each of the 48 fuel assembly cases. Therefore, when the perturbation methods are used, 240 perturbation cases (=[1 base case +47 perturbed cases] × [k_∞ + 4 reactivity coefficients]) are calculated through a single Monte Carlo calculation.

Figures 9–13 show the k-eigenvalue and reactivity coefficient results of perturbation methods. Table 4 shows the summary of perturbation results for 240 perturbed states of PLUS7 fuel assembly. The calculations are performed by using 25 active cycles with 8 million histories per cycle. The statistical uncertainty of k-eigenvalue is about 10 pcm. In every case, the root mean square (RMS) and maximum errors of EPP method are smaller than the error of conventional perturbation methods, which are DOS, AWP(DOS), and AWP(CS). Especially, the EPP method result shows a significantly improved accuracy for the k-eigenvalue change in Figure 9. For the k-eigenvalue change of 48 branch cases, the RMS and maximum errors of EPP method are 22 pcm and 59 pcm, respectively. However, in several cases, it is also observed that the error of EPP method is larger than the one of AWP method. It is due to the statistical uncertainty of Monte Carlo simulations which also affect to the DP results. By comparing the results of EPP method and AWP(CS) method, one can observe the effect of using adjoint flux of perturbed system instead of the adjoint flux of unperturbed system.

Fig. 9. k-eigenvalue change results of perturbation methods for multiple parameter perturbation.

Fig. 10. Fuel temperature coefficients results of perturbation methods for multiple parameter perturbation.

Fig. 11. Xenon worth results of perturbation methods for the multiple parameter perturbation.

Fig. 12. Moderator density coefficients results of perturbation methods for the multiple parameter perturbation.

Fig. 13. Boron worth results of perturbation methods for the multiple parameter perturbation.

Table 5 shows the computational costs of perturbation methods for the multiple parameter perturbation in PLUS7 fuel assembly. When perturbation methods are used, the memory usage increased from 291 MB to 433 MB. It is because additional cross section data are used in perturbation cases to consider fuel temperature changes in perturbed systems. When perturbation methods are used, the computational time is effectively reduced. The reduction of computational time is larger in cases using DOS compared to cases using CS. It is because, in terms of computational burden, multiplication and division operations used in CS method are more expensive than addition and subtraction operations used in DOS method.

4. Conclusion

In this work, the PIFP method is developed for the EPT in Monte Carlo transport. It allows to use the EPT by computing the adjoint flux of perturbed system during a forward Monte Carlo calculation. The EPP method is implemented in the Monte Carlo code MCS. For the Godiva and PLUS7 fuel assembly, k-eigenvalue change results of perturbation methods, which are DOS, AWP(DOS), AWP(CS), and EPP, are compared against the reference DP result. The EPP result shows a significantly improved accuracy compared with the conventional perturbation methods. It is also observed that the EPP method takes 15.1 times faster than the DP method to compute 29 boron concentration change cases in PLUS7 fuel assembly. For the 240 perturbation cases considering multiple parameter perturbation, the computational time is reduced by 31.7 times when the EPP method is used instead of DP method. It is because DP method requires 240 Monte Carlo calculations, while the EPP method can obtain k-eigenvalues of 240 perturbation cases by running a single Monte Carlo calculation. Additionally, EPP method shows the smallest error compared with the conventional perturbation methods especially for the estimation of k-eigenvalue change, xenon worth and boron worth. For the future work, the EPP method can be applied to the stochastic sampling method to reduce the number of Monte Carlo simulations effectively during the sensitivity analysis and uncertainty quantification.

Funding

This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. RS-2025-02353161).

Conflicts of interest

The authors have nothing to disclose.

Data availability statement

This article has no associated data.

Author contribution statement

Conceptualization, Yunki Jo and Deokjung Lee; Methodology, Yunki Jo; Software, Yunki Jo; Validation, Yunki Jo; Formal Analysis, Yunki Jo; Investigation, Yunki Jo; Writing – Original Draft Preparation, Yunki Jo; Writing – Review & Editing, Yunki Jo and Deokjung Lee; Visualization, Yunki Jo; Supervision, Deokjung Lee;

References

All Tables

All Figures

	Fig. 1. Analytic results of perturbation methods for Σa1 change in the 2-group infinite homogeneous problem.
In the text

	Fig. 2. Analytic results of perturbation methods for Σs, 1 → 2 change in the 2-group infinite homogeneous problem.
In the text

	Fig. 3. Results of perturbation methods for local density change in the shell region of Godiva: k-eigenvalue change (up) and error (down).
In the text

	Fig. 4. Results of perturbation methods for local density change in center sphere region of Godiva: k-eigenvalue change (up) and error (down).
In the text

	Fig. 5. Results of perturbation methods for global density change in Godiva: k-eigenvalue change (up) and error (down).
In the text

	Fig. 6. Results of perturbation methods for boron concentration change in PLUS7 fuel assembly: k-eigenvalue change (up) and error (down).
In the text

	Fig. 7. PLUS7 fuel assembly models having different number of gadolinia rods.
In the text

	Fig. 8. Results of perturbation methods for the number of gadolinia rods change in PLUS7 fuel assembly.
In the text

	Fig. 9. k-eigenvalue change results of perturbation methods for multiple parameter perturbation.
In the text

	Fig. 10. Fuel temperature coefficients results of perturbation methods for multiple parameter perturbation.
In the text

	Fig. 11. Xenon worth results of perturbation methods for the multiple parameter perturbation.
In the text

	Fig. 12. Moderator density coefficients results of perturbation methods for the multiple parameter perturbation.
In the text

	Fig. 13. Boron worth results of perturbation methods for the multiple parameter perturbation.
In the text