Mesh-free Semi-quantitative Variance Underestimation Elimination Method in Monte Caro algorithm

Mesh-free Semi-quantitative Variance Underestimation Elimination Method in Monte Caro algorithm 1 2 first-author https://orcid.org/0000-0002-5504-4436 2 1 corresp * lizeguang@mail.tsinghua.edu.cn * lizeguang@mail.tsinghua.edu.cn lizeguang@mail.tsinghua.edu.cn https://orcid.org/0000-0002-3457-8471 2 2 2 1 Department of Engineering Physics, Tsinghua University, Beijing 100084, China Department of Engineering Physics, Tsinghua University, Beijing 100084, China China Nuclear Power Engineering Corporation, Beijing 100840, China China Nuclear Power Engineering Corporation, Beijing 100840, China The inter-cycle correlation of fission source distributions (FSDs) in the Monte Carlo Power Iteration process results in variance underestimation of tallied physical quantities, especially in large local tallies. This study provides a mesh-free Semi-quantitative Variance Underestimation Elimination (SeVUE) method to obtain a credible confidence interval for the tallied results. This method comprises two procedures: Estimation and Elimination. The FSD inter-cycle correlation length is estimated in the Estimation procedure using the Sliced Wasserstein distance algorithm. The batch method was then used in the Elimination procedure. The FSD inter-cycle correlation length was proved to be the optimum batch length to eliminate the variance underestimation problem. We exemplified this method using the OECD sphere array model and 3D PWR BEAVRS model. The results showed that the average variance underestimation ratios of local tallies declined from 37% and 87% to within ±5% in these models. Monte Carlo algorithm Power Iteration process Inter-cycle correlation Variance Underestimation Sliced Wasserstein distance Introduction 1 The Monte Carlo (MC) algorithm is commonly used in high-resolution modeling and high-fidelity analysis of reactor core models and other nuclear facilities. The Power Iteration (PI) process was adopted to solve the k-eigenvalue equations [ 1 1 ] and transient fixed-source equations [ 2 2 ]. Physical quantities such as k-eigenvalue and power densities are tallied cycle by cycle after source convergence. Commonly, these tally results are processed by assuming that they are independent and uncorrelated between cycles. However, in the PI process, the Fission Source Bank is transmitted between cycles, resulting in an inter-cycle correlation of the fission source distributions (FSDs). This inter-cycle correlation of the FSDs further led to the inter-cycle correlation of the tally results [ 3 3 ]. Thus, the calculated variance of the tally results is lower than the real variance owing to this inter-cycle correlation. This phenomenon is also called the MC variance underestimation phenomenon. In the PWR full-core power distribution calculations, the inter-cycle correlations can result in a variance underestimation ratio of 90%. Local tallies with relatively large sizes suffer the most from this problem [ 4 4 ]. The variance underestimation phenomenon belongs to the MC source convergence problem. Researchers have attempted to solve these problems [ 5 5 - 9 9 5-9 ]. To solve the variance underestimation problem, there are various proposed methods in literature, which can be divided into three categories. The first category involves modifying the PI process to reduce the inter-cycle correlations of the tally results. Brissenden and Garrick (1986) proposed the super history method [ 10 10 ]. This method simulates neutrons for i generations per cycle, and the physical quantities are tallied at the cycle end. This method is intuitively effective, but the number of neutron generations per cycle is set based on experience. Herman et al. (2014) proposed using CMFD feedback to adjust the neutron weight in the PI process to reduce FSD inter-cycle correlations [ 11 11 - 12 12 11-12 ]. However, it changes the fundamental mode of the FSD in the MC simulations. Gelbard and Prael (1990) proposed a batch method to de-correlate tally results by combining a certain number of active cycles into batches [ 13 13 - 14 14 13-14 ]. The physical quantities were tallied at the end of each batch. As observed with the super history method, setting an appropriate batch length is difficult. The second category quantifies the underestimation ratio by post-processing the auto-correlated tally results. Several mathematical models are used to analyze the auto-correlated tally results, such as the ARMA(2,1) model, orthonormal weighted standardized time series model (OWSTS), and fractional Brownian motion process [ 15 15 - 17 17 15-17 ]. These mathematical models are efficient in variance underestimation correction but lack a theoretical basis in the MC PI theory. The third category theoretically predicts the inter-cycle correlation. The prediction methods approximate MC simulations by using a surrogate model, calculating the correlation coefficients, and correcting the calculated variances. These methods are typically based on specific mesh grids or homogenized models [ 18 18 - 19 19 18-19 ]. Thus, the approaches to solving the variance underestimation problem are always demanding for users by setting appropriate mesh grids or tricky coefficients. This study proposes a mesh-free Semi-quantitative Variance Underestimation Elimination (SeVUE) method to obtain a credible confidence interval for the tallied results. This method comprised two procedures: Estimation and Elimination. The FSD inter-cycle correlation length was estimated in the Estimation procedure using the Sliced Wasserstein (SW) distance algorithm, the batch method was used in the Elimination procedure, and the FSD inter-cycle correlation length was set as the batch length. Local tallies with relatively large sizes were of interest in this study. The research was conducted on the RMC platform [ 20 20 ]. The paper proceeds as follows: Sect. 2 briefly introduces the batch method and SW-based FSD inter-cycle correlation observation method; Sect. 3 illustrates the SeVUE algorithm; Sect. 4 presents the validation calculations for the OECD sphere array model and the 3D PWR BEAVRS model; and Sect. 5 presents the conclusions. Theory 1 The batch method is summarized in Sect. 2.1, and the SW distance algorithm and mesh-free FSD inter-cycle correlation observation method are introduced in Sects. 2.2 and 2.3. The optimum batch length is related to the FSD inter-cycle correlation length, discussed in Sect. 2.4. The batch method 2 The batch method de-correlates the tally results by combining a certain number of active cycles into batches. In the MC simulation, <math id="M1"><mi>Q</mi></math> is assumed to be the physical quantity of interest. Next, its conventional MC estimate value <math id="M2"><mover accent="true"><mi>Q</mi><mo>¯</mo></mover></math> and estimate variance <math id="M3"><mrow><msubsup><mi>σ</mi><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow></mrow></math> are assumed to be [ 3 3 ] <math display="block" id="M4"><mrow><mtable><mtr><mtd><mrow><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>=</mo><mfrac><mn>1</mn><mi>N</mi></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><msup><mi>Q</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup></mrow></mstyle></mrow></mtd></mtr></mtable></mrow></math> <math display="block" id="M5"><mrow><mtable><mtr><mtd><mrow><msubsup><mi>σ</mi><mtext>s</mtext><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><mi>N</mi><mrow><mo>(</mo><mrow><mi>N</mi><mo>−</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><msup><mrow><mrow><mo>(</mo><mrow><msup><mi>Q</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup><mo>−</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></mstyle></mrow></mtd></mtr></mtable></mrow></math> where <math id="M6"><mrow><msup><mi>Q</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup></mrow></math> is the estimate in cycle i , and <math id="M7"><mi>N</mi></math> is the number of active cycles. The apparent variance of <math id="M8"><mover accent="true"><mi>Q</mi><mo>¯</mo></mover></math> , <math id="M9"><mrow><msubsup><mi>σ</mi><mtext>a</mtext><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow></mrow></math> is defined as the mathematical expectation of the estimated variance. <math display="block" id="M10"><mrow><mtable><mtr><mtd><mrow><msubsup><mi>σ</mi><mtext>a</mtext><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow><mo>=</mo><mi>E</mi><mrow><mo>(</mo><mrow><msubsup><mi>σ</mi><mtext>s</mtext><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mi>N</mi></mfrac><mi>E</mi><mrow><mo>(</mo><mrow><msup><mrow><mrow><mo>(</mo><mrow><msup><mi>Q</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow><mo>)</mo></mrow><mo>−</mo><mfrac><mn>1</mn><mi>N</mi></mfrac><msup><mi>E</mi><mn>2</mn></msup><mrow><mo>(</mo><mrow><msup><mi>Q</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup></mrow><mo>)</mo></mrow><mo>−</mo><mfrac><mn>1</mn><mrow><msup><mi>N</mi><mn>2</mn></msup></mrow></mfrac><mfrac><mn>1</mn><mrow><mi>N</mi><mo>−</mo><mn>1</mn></mrow></mfrac><mstyle displaystyle="true"><munder><mo>∑</mo><mi>i</mi></munder><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>j</mi><mo>≠</mo><mi>i</mi></mrow></munder><mrow><mtext>cov</mtext></mrow></mstyle></mrow></mstyle><mrow><mo>(</mo><mrow><msup><mi>Q</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup><mo>,</mo><msup><mi>Q</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup></mrow><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mrow></math> If we define <math id="M11"><mrow><msup><mi>σ</mi><mn>2</mn></msup><mrow><mo>(</mo><mrow><msup><mi>Q</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup></mrow><mo>)</mo></mrow><mo>=</mo><mi>E</mi><mrow><mo>(</mo><mrow><msup><mrow><mrow><mo>(</mo><mrow><msup><mi>Q</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow><mo>)</mo></mrow><mo>−</mo><msup><mi>E</mi><mn>2</mn></msup><mrow><mo>(</mo><mrow><msup><mi>Q</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup></mrow><mo>)</mo></mrow></mrow></math> , then Eq. (3) can be simplified as <math display="block" id="M12"><mrow><mtable><mtr><mtd><mrow><msubsup><mi>σ</mi><mtext>a</mtext><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mi>N</mi></mfrac><msup><mi>σ</mi><mn>2</mn></msup><mrow><mo>(</mo><mrow><msup><mi>Q</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup></mrow><mo>)</mo></mrow><mo>−</mo><mfrac><mn>1</mn><mrow><msup><mi>N</mi><mn>2</mn></msup></mrow></mfrac><mfrac><mn>1</mn><mrow><mi>N</mi><mo>−</mo><mn>1</mn></mrow></mfrac><mstyle displaystyle="true"><munder><mo>∑</mo><mi>i</mi></munder><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>j</mi><mo>≠</mo><mi>i</mi></mrow></munder><mrow><mtext>cov</mtext></mrow></mstyle></mrow></mstyle><mrow><mo>(</mo><mrow><msup><mi>Q</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup><mo>,</mo><msup><mi>Q</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup></mrow><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mrow></math> The true variance <math id="M13"><mrow><msubsup><mi>σ</mi><mtext>r</mtext><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow></mrow></math> of <math id="M14"><mover accent="true"><mi>Q</mi><mo>¯</mo></mover></math> is expressed by Eq. (5). <math display="block" id="M15"><mrow><mtable><mtr><mtd><mrow><msubsup><mi>σ</mi><mtext>r</mtext><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow><mo>=</mo><mi>E</mi><mrow><mo>(</mo><mrow><msup><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mn>2</mn></msup></mrow><mo>)</mo></mrow><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mi>E</mi><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup><mo>=</mo><mfrac><mn>1</mn><mi>N</mi></mfrac><msup><mi>σ</mi><mn>2</mn></msup><mrow><mo>(</mo><mrow><msup><mi>Q</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup></mrow><mo>)</mo></mrow><mo>+</mo><mfrac><mn>1</mn><mrow><msup><mi>N</mi><mn>2</mn></msup></mrow></mfrac><mstyle displaystyle="true"><munder><mo>∑</mo><mi>i</mi></munder><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>j</mi><mo>≠</mo><mi>i</mi></mrow></munder><mrow><mtext>cov</mtext></mrow></mstyle></mrow></mstyle><mrow><mo>(</mo><mrow><msup><mi>Q</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup><mo>,</mo><msup><mi>Q</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup></mrow><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mrow></math> The difference between <math id="M16"><mrow><msubsup><mi>σ</mi><mtext>r</mtext><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow></mrow></math> and <math id="M17"><mrow><msubsup><mi>σ</mi><mtext>a</mtext><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow></mrow></math> is <math display="block" id="M18"><mrow><mtable><mtr><mtd><mrow><msubsup><mi>σ</mi><mtext>r</mtext><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow><mo>−</mo><msubsup><mi>σ</mi><mtext>a</mtext><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><mi>N</mi><mrow><mo>(</mo><mrow><mi>N</mi><mo>−</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></mfrac><mstyle displaystyle="true"><munder><mo>∑</mo><mi>i</mi></munder><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>j</mi><mo>≠</mo><mi>i</mi></mrow></munder><mrow><mtext>cov</mtext></mrow></mstyle></mrow></mstyle><mrow><mo>(</mo><mrow><msup><mi>Q</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup><mo>,</mo><msup><mi>Q</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup></mrow><mo>)</mo></mrow><mo>.</mo></mrow></mtd></mtr></mtable></mrow></math> Thus, the correlation between the tally results of every cycle is the direct cause of the variance underestimation problem. Furthermore, we used the batch method by combining a certain number of active cycles into batches and tallying the physical quantities at the end of each batch. Suppose that the batch length is <math id="M19"><mrow><msub><mi>l</mi><mtext>b</mtext></msub></mrow></math> and the number of batches is <math id="M20"><mrow><msub><mi>n</mi><mtext>b</mtext></msub></mrow></math> . The coefficients confront Eq. (7). <math display="block" id="M21"><mrow><mtable><mtr><mtd><mrow><msub><mi>l</mi><mtext>b</mtext></msub><mo>×</mo><msub><mi>n</mi><mtext>b</mtext></msub><mo>=</mo><mi>N</mi><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo></mrow></mtd></mtr></mtable></mrow></math> We defined <math id="M22"><mrow><msubsup><mi>Q</mi><mtext>b</mtext><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup></mrow></math> as the estimate in batch j and <math id="M23"><mrow><msubsup><mi>Q</mi><mtext>b</mtext><mrow><mi>j</mi><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup></mrow></math> as the estimate of cycle i in batch j . The estimated value <math id="M24"><mrow><msub><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mtext>b</mtext></msub></mrow></math> and sample variance <math id="M25"><mrow><msubsup><mi>σ</mi><mrow><mtext>bs</mtext></mrow><mn>2</mn></msubsup><mrow><mo>(</mo><mrow><msub><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mtext>b</mtext></msub></mrow><mo>)</mo></mrow></mrow></math> can be as Eq. (8-9). <math display="block" id="M26"><mrow><mtable><mtr><mtd><mrow><msub><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mtext>b</mtext></msub><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>n</mi><mtext>b</mtext></msub></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mrow><msub><mi>n</mi><mtext>b</mtext></msub></mrow></munderover><mrow><msubsup><mi>Q</mi><mtext>b</mtext><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup></mrow></mstyle><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>n</mi><mtext>b</mtext></msub></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mrow><msub><mi>n</mi><mtext>b</mtext></msub></mrow></munderover><mrow><mrow><mo>(</mo><mrow><mfrac><mn>1</mn><mrow><msub><mi>l</mi><mtext>b</mtext></msub></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mrow><msub><mi>l</mi><mtext>b</mtext></msub></mrow></munderover><mrow><msubsup><mi>Q</mi><mtext>b</mtext><mrow><mi>j</mi><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup></mrow></mstyle></mrow><mo>)</mo></mrow></mrow></mstyle><mo>=</mo><mfrac><mn>1</mn><mi>N</mi></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><msup><mi>Q</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup></mrow></mstyle><mo>=</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover></mrow></mtd></mtr></mtable></mrow></math> <math display="block" id="M27"><mrow><mtable><mtr><mtd><mrow><msubsup><mi>σ</mi><mrow><mtext>bs</mtext></mrow><mn>2</mn></msubsup><mrow><mo>(</mo><mrow><msub><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mi>b</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>n</mi><mtext>b</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>n</mi><mtext>b</mtext></msub><mo>−</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mrow><msub><mi>n</mi><mtext>b</mtext></msub></mrow></munderover><mrow><msup><mrow><mrow><mo>(</mo><mrow><msubsup><mi>Q</mi><mtext>b</mtext><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup><mo>−</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></mstyle></mrow></mtd></mtr></mtable></mrow></math> The apparent variance of <math id="M28"><mrow><msub><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mtext>b</mtext></msub></mrow></math> , <math id="M29"><mrow><msubsup><mi>σ</mi><mrow><mtext>ba</mtext></mrow><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow></mrow></math> is defined as the expected value of the sample variance. <math display="block" id="M30"><mrow><mtable><mtr><mtd><mrow><msubsup><mi>σ</mi><mrow><mtext>ba</mtext></mrow><mn>2</mn></msubsup><mrow><mo>(</mo><mrow><msub><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mi>b</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><mi>E</mi><mrow><mo>(</mo><mrow><msubsup><mi>σ</mi><mrow><mtext>bs</mtext></mrow><mn>2</mn></msubsup><mrow><mo>(</mo><mrow><msub><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mtext>b</mtext></msub></mrow><mo>)</mo></mrow></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>n</mi><mtext>b</mtext></msub></mrow></mfrac><msup><mi>σ</mi><mn>2</mn></msup><mrow><mo>(</mo><mrow><msubsup><mi>Q</mi><mi>b</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup></mrow><mo>)</mo></mrow><mo>+</mo><mfrac><mn>1</mn><mrow><msub><mi>n</mi><mtext>b</mtext></msub><mo>−</mo><mn>1</mn></mrow></mfrac><mfrac><mn>1</mn><mrow><msubsup><mi>n</mi><mtext>b</mtext><mn>2</mn></msubsup></mrow></mfrac><mstyle displaystyle="true"><munder><mo>∑</mo><mi>i</mi></munder><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>j</mi><mo>≠</mo><mi>i</mi></mrow></munder><mrow><mtext>cov</mtext></mrow></mstyle></mrow></mstyle><mrow><mo>(</mo><mrow><msubsup><mi>Q</mi><mtext>b</mtext><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup><mo>,</mo><msubsup><mi>Q</mi><mtext>b</mtext><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup></mrow><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mrow></math> where <math id="M31"><mrow><msup><mi>σ</mi><mn>2</mn></msup><mrow><mo>(</mo><mrow><msubsup><mi>Q</mi><mtext>b</mtext><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup></mrow><mo>)</mo></mrow><mo>=</mo><mi>E</mi><mrow><mo>(</mo><mrow><msup><mrow><mrow><mo>(</mo><mrow><msubsup><mi>Q</mi><mtext>b</mtext><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow><mo>)</mo></mrow><mo>−</mo><msup><mi>E</mi><mn>2</mn></msup><mrow><mo>(</mo><mrow><msubsup><mi>Q</mi><mtext>b</mtext><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup></mrow><mo>)</mo></mrow></mrow></math> . Then, the difference between <math id="M32"><mrow><msubsup><mi>σ</mi><mrow><mtext>br</mtext></mrow><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow></mrow></math> and <math id="M33"><mrow><msubsup><mi>σ</mi><mrow><mtext>ba</mtext></mrow><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow></mrow></math> is <math display="block" id="M34"><mrow><mtable><mtr><mtd><mrow><msubsup><mi>σ</mi><mrow><mtext>br</mtext></mrow><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow><mo>−</mo><msubsup><mi>σ</mi><mrow><mtext>ba</mtext></mrow><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>n</mi><mtext>b</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>n</mi><mtext>b</mtext></msub><mo>−</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></mfrac><mstyle displaystyle="true"><munder><mo>∑</mo><mi>i</mi></munder><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>j</mi><mo>≠</mo><mi>i</mi></mrow></munder><mrow><mtext>cov</mtext></mrow></mstyle></mrow></mstyle><mrow><mo>(</mo><mrow><msubsup><mi>Q</mi><mtext>b</mtext><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup><mo>,</mo><msubsup><mi>Q</mi><mtext>b</mtext><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup></mrow><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mrow></math> The correlation among batch tally results is generally less significant than that among cycle tally results in the PI process, owing to the weaker FSD inter-batch correlation. Thus, the difference between <math id="M35"><mrow><msubsup><mi>σ</mi><mrow><mtext>br</mtext></mrow><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow></mrow></math> and <math id="M36"><mrow><msubsup><mi>σ</mi><mrow><mtext>ba</mtext></mrow><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow></mrow></math> was smaller than that between <math id="M37"><mrow><msubsup><mi>σ</mi><mtext>r</mtext><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow></mrow></math> and <math id="M38"><mrow><msubsup><mi>σ</mi><mtext>a</mtext><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow></mrow></math> . However, setting an appropriate batch length is challenging. A large batch length can almost eliminate underestimation, whereas the computational cost can be too high because it requires a minimum number of batches in common practice. A small batch length is appropriate for practical use, but the magnitude of variance underestimation elimination is unknown. Determining the optimum batch length was of interest in this study. Sliced Wasserstein distance algorithm 2 The Sliced Wasserstein (SW) distance was developed from Optimum Transport theory [ 21 21 ]. Its initial definition is the minimum "cost" to transform the 1D distribution <math id="M39"><mstyle mathvariant="bold-italic"><mi>p</mi></mstyle></math> into <math id="M40"><mstyle mathvariant="bold-italic"><mi>q</mi></mstyle></math> , as shown in Fig. 1 Fig. 1 . It can be efficiently calculated by <math display="block" id="M41"><mrow><msub><mrow><mtext>SW</mtext></mrow><mi>k</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>p</mi></mstyle><mo>,</mo><mstyle mathvariant="bold-italic"><mi>q</mi></mstyle></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mi>N</mi></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mtext>d</mtext></mstyle><mrow><mo>(</mo><mrow><msub><mi>x</mi><mi>i</mi></msub><mo>,</mo><mtext> </mtext><msub><mi>y</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where k is the X / Y / Z direction, <math id="M42"><mrow><mstyle mathvariant="bold-italic"><mi>p</mi></mstyle><mo>,</mo><mstyle mathvariant="bold-italic"><mi>q</mi></mstyle></mrow></math> are two 1D distributions, <math id="M43"><mrow><msub><mi>x</mi><mi>i</mi></msub><mo>∈</mo><mstyle mathvariant="bold-italic"><mi>p</mi></mstyle></mrow></math> , <math id="M44"><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>∈</mo><mstyle mathvariant="bold-italic"><mi>q</mi></mstyle></mrow></math> are the coordinates after sorting the neutron points in the X / Y / Z direction, and N is the number of particles per cycle in the MC simulations. In this paper, the 1-norm distance is used to calculate <math id="M45"><mrow><mi>d</mi><mrow><mo>(</mo><mrow><mi>x</mi><mo>,</mo><mi>y</mi></mrow><mo>)</mo></mrow></mrow></math> : <math display="block" id="M46"><mrow><mtable><mtr><mtd><mrow><mtext>d</mtext><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>x</mi></mstyle><mo>,</mo><mstyle mathvariant="bold-italic"><mi>y</mi></mstyle></mrow><mo>)</mo></mrow><mo>=</mo><mrow><mo>‖</mo><mrow><mstyle mathvariant="bold-italic"><mi>x</mi></mstyle><mo>−</mo><msub><mstyle mathvariant="bold-italic"><mi>y</mi></mstyle><mn>1</mn></msub></mrow><mo>‖</mo></mrow><mo>=</mo><mo>∑</mo><mrow><mo>|</mo><mrow><msub><mi>x</mi><mi>i</mi></msub><mo>−</mo><msub><mi>y</mi><mi>i</mi></msub></mrow><mo>|</mo></mrow></mrow></mtd></mtr></mtable><mo>.</mo></mrow></math> The SW distance algorithm is a commonly used mesh-free algorithm for measuring the difference between discrete vectors. Guo et al. (2021) introduced this algorithm into the scope of MC source convergence diagnosis [ 22 22 ]. Shen et al. (2022) adopted this algorithm to measure fission source error in the MC PI process [ 23 23 ]. In the calculation, the sorting procedure uses the fast sorting algorithm costs <math id="M47"><mrow><mi>O</mi><mrow><mo>(</mo><mrow><mi>N</mi><mtext>log</mtext><mrow><mo>(</mo><mi>N</mi><mo>)</mo></mrow></mrow><mo>)</mo></mrow></mrow></math> operations, and the SW distance calculation costs <math id="M48"><mrow><mi>O</mi><mrow><mo>(</mo><mi>N</mi><mo>)</mo></mrow></mrow></math> operations. Thus, the total cost is <math id="M49"><mrow><mi>O</mi><mrow><mo>(</mo><mrow><mi>N</mi><mtext>log</mtext><mrow><mo>(</mo><mi>N</mi><mo>)</mo></mrow></mrow><mo>)</mo></mrow></mrow></math> . The efficiency penalty is approximately 2% in MC criticality calculations for the PWR full-core model [ 22 22 ]. SW-based FSD inter-cycle correlation observation method 2 The PI method in the MC algorithm is expressed by Eq. (14) in the discretization. <math display="block" id="M52"><mrow><mtable><mtr><mtd><mrow><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msup><mo>=</mo><mfrac><mrow><mstyle mathvariant="bold-italic"><mi>H</mi></mstyle><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup></mrow><mrow><msup><mi>k</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup></mrow></mfrac><mo>+</mo><msup><mstyle mathvariant="bold-italic"><mi>ϵ</mi></mstyle><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup></mrow></mtd></mtr></mtable><mo>,</mo></mrow></math> where <math id="M53"><mrow><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup></mrow></math> is the fission source vector in cycle I, <math id="M54"><mstyle mathvariant="bold-italic"><mi>H</mi></mstyle></math> is the fission matrix, <math id="M55"><mrow><msup><mstyle mathvariant="bold-italic"><mi>ϵ</mi></mstyle><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup></mrow></math> is the stochastic error component, and <math id="M56"><mrow><msup><mi>k</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup></mrow></math> is the estimation of fundamental mode eigenvalue in cycle i . The number of particles per cycle is defined as m , the fundamental mode eigenvalue, and the normalized eigenvector of fission matrix <math id="M57"><mstyle mathvariant="bold-italic"><mi>H</mi></mstyle></math> is defined as <math id="M58"><mrow><msub><mi>k</mi><mn>0</mn></msub></mrow></math> and <math id="M59"><mrow><msub><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mn>0</mn></msub></mrow></math> . Eq. (15) can be expressed as <math display="block" id="M60"><mrow><mtable><mtr><mtd><mrow><mi>H</mi><msub><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mn>0</mn></msub><mo>=</mo><msub><mi>k</mi><mn>0</mn></msub><msub><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mn>0</mn></msub></mrow></mtd></mtr></mtable><mo>.</mo></mrow></math> The fission source error in cycle i is introduced as <math display="block" id="M61"><mrow><mtable><mtr><mtd><mrow><msup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup><mo>=</mo><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup><mo>−</mo><mi>m</mi><msub><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mn>0</mn></msub><mo>.</mo></mrow></mtd></mtr></mtable></mrow></math> By substituting and considering the magnitude of the error terms [ 10 10 ], we disregarded some terms and obtained <math display="block" id="M62"><mrow><mtable><mtr><mtd><mrow><msup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup><mo>≅</mo><mstyle mathvariant="bold-italic"><mi>A</mi></mstyle><msup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msup><mo>−</mo><mfrac><mrow><mstyle mathvariant="bold-italic"><mi>τ</mi></mstyle><mstyle mathvariant="bold-italic"><mi>H</mi></mstyle></mrow><mrow><msub><mi>k</mi><mn>0</mn></msub><mi>m</mi></mrow></mfrac><msup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msup><mstyle mathvariant="bold-italic"><mi>A</mi></mstyle><msup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msup><mo>+</mo><msup><mstyle mathvariant="bold-italic"><mi>ϵ</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msup><mo> </mo><mo>,</mo></mrow></mtd></mtr></mtable></mrow></math> where <math id="M63"><mstyle mathvariant="bold-italic"><mi>τ</mi></mstyle></math> is (1, 1, …, 1), and <math id="M64"><mstyle mathvariant="bold-italic"><mi>A</mi></mstyle></math> is the noise propagation matrix, written as Eq. (18). <math display="block" id="M65"><mrow><mtable><mtr><mtd><mrow><mstyle mathvariant="bold-italic"><mi>A</mi></mstyle><mo>=</mo><mfrac><mrow><mstyle mathvariant="bold-italic"><mi>I</mi></mstyle><mo>−</mo><msub><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mn>0</mn></msub><mstyle mathvariant="bold-italic"><mi>τ</mi></mstyle></mrow><mrow><msub><mi>k</mi><mn>0</mn></msub></mrow></mfrac><mstyle mathvariant="bold-italic"><mi>H</mi></mstyle></mrow></mtd></mtr></mtable></mrow></math> We iterated Eq. (17) i times and obtained Eq. (19). <math display="block" id="M66"><mrow><mtable><mtr><mtd><mrow><msup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup><mo>≅</mo><msup><mstyle mathvariant="bold-italic"><mi>A</mi></mstyle><mi>i</mi></msup><msup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mrow><mrow><mo>(</mo><mn>0</mn><mo>)</mo></mrow></mrow></msup><mo>+</mo><mi>O</mi><mrow><mo>(</mo><mrow><msup><mi>m</mi><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow><mo>)</mo></mrow><mo>+</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></munderover><msup><mstyle mathvariant="bold-italic"><mi>A</mi></mstyle><mi>j</mi></msup><msup><mstyle mathvariant="bold-italic"><mi>ϵ</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>−</mo><mi>j</mi></mrow><mo>)</mo></mrow></mrow></msup><mo> </mo></mrow></mtd></mtr></mtable></mrow></math> The fission source error consists of three components: the initial error, bias error, and random error terms. The FSD random error term is the dominant component after the source convergence, as expressed in Eq. (20) [ 24 24 ]. <math display="block" id="M67"><mrow><mtable><mtr><mtd><mrow><msubsup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mtext>R</mtext><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup><mo>≅</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><msup><mstyle mathvariant="bold-italic"><mi>A</mi></mstyle><mi>j</mi></msup><msup><mstyle mathvariant="bold-italic"><mi>ϵ</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>−</mo><mi>j</mi></mrow><mo>)</mo></mrow></mrow></msup></mrow></mstyle><mo> </mo></mrow></mtd></mtr></mtable><mo>,</mo></mrow></math> where <math id="M68"><mrow><msubsup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mi>R</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup></mrow></math> is the fission source random error term of cycle i , <math id="M69"><mstyle mathvariant="bold-italic"><mi>A</mi></mstyle></math> is the error propagation matrix, and <math id="M70"><mrow><msup><mstyle mathvariant="bold-italic"><mi>ϵ</mi></mstyle><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup></mrow></math> is the stochastic error component induced in cycle j . The FSD random error term is inter-cycle correlated because it is composed of stochastic error components induced in adjacent cycles. This is the theoretical basis of the FSD inter-cycle correlation. The fission error difference between adjacent cycles can be written as <math display="block" id="M71"><mrow><mtable><mtr><mtd><mrow><msubsup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mtext>R</mtext><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup><mo>−</mo><msubsup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mtext>R</mtext><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msubsup><mo>=</mo><msup><mstyle mathvariant="bold-italic"><mi>ϵ</mi></mstyle><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>A</mi></mstyle><mo>−</mo><mstyle mathvariant="bold-italic"><mi>I</mi></mstyle></mrow><mo>)</mo></mrow></mrow></mstyle><msup><mstyle mathvariant="bold-italic"><mi>A</mi></mstyle><mrow><mi>j</mi><mo>−</mo><mn>1</mn></mrow></msup><msup><mstyle mathvariant="bold-italic"><mi>ϵ</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>−</mo><mi>j</mi></mrow><mo>)</mo></mrow></mrow></msup><mo>,</mo></mrow></mtd></mtr></mtable></mrow></math> where <math id="M72"><mstyle mathvariant="bold-italic"><mi>I</mi></mstyle></math> is the identity matrix. In fission systems, noise propagation matrix <math id="M73"><mstyle mathvariant="bold-italic"><mi>A</mi></mstyle></math> is similar to the identity matrix. By comparing the scalar error differences, we obtained <math display="block" id="M74"><mrow><mtable><mtr><mtd><mrow><mrow><mo>‖</mo><mrow><msubsup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mtext>R</mtext><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup></mrow><mo>‖</mo></mrow><mo>></mo><mrow><mo>‖</mo><mrow><msubsup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mtext>R</mtext><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup><mo>−</mo><msubsup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mtext>R</mtext><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msubsup></mrow><mo>‖</mo></mrow></mrow></mtd></mtr></mtable><mo>,</mo></mrow></math> where <math id="M75"><mrow><msup><mi>ε</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup></mrow></math> is the scalar fission source error in cycle i . Furthermore, the fission error difference between cycles can be written as <math display="block" id="M76"><mrow><mtable><mtr><mtd><mrow><msubsup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mtext>R</mtext><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup><mo>−</mo><msubsup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mtext>R</mtext><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>−</mo><mi>j</mi></mrow><mo>)</mo></mrow></mrow></msubsup><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>j</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><msup><mstyle mathvariant="bold-italic"><mi>A</mi></mstyle><mi>k</mi></msup><msup><mstyle mathvariant="bold-italic"><mi>ϵ</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>−</mo><mi>k</mi></mrow><mo>)</mo></mrow></mrow></msup></mrow></mstyle><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>t</mi><mo>=</mo><mi>j</mi></mrow><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><mrow><mo>(</mo><mrow><msup><mstyle mathvariant="bold-italic"><mi>A</mi></mstyle><mi>t</mi></msup><mo>−</mo><msup><mstyle mathvariant="bold-italic"><mi>A</mi></mstyle><mrow><mi>t</mi><mo>−</mo><mi>j</mi></mrow></msup></mrow><mo>)</mo></mrow></mrow></mstyle><msup><mstyle mathvariant="bold-italic"><mi>ϵ</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>−</mo><mi>t</mi></mrow><mo>)</mo></mrow></mrow></msup><mo>.</mo></mrow></mtd></mtr></mtable></mrow></math> Similarly, the scalar error differences are expected to be <math display="block" id="M77"><mrow><mrow><mo>‖</mo><mrow><msubsup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mtext>R</mtext><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup></mrow><mo>‖</mo></mrow><mtable><mtr><mtd><mrow><mo>></mo><mrow><mo>‖</mo><mrow><msubsup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mtext>R</mtext><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup><mo>−</mo><msubsup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mtext>R</mtext><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>−</mo><mi>j</mi></mrow><mo>)</mo></mrow></mrow></msubsup></mrow><mo>‖</mo></mrow><mo>></mo><mo>…</mo><mo>></mo><mrow><mo>‖</mo><mrow><msubsup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mtext>R</mtext><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup><mo>−</mo><msubsup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mtext>R</mtext><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msubsup></mrow><mo>‖</mo></mrow><mo>,</mo><mtext> </mtext><mi>j</mi><mo>></mo><mn>1</mn></mrow></mtd></mtr></mtable></mrow></math> The scalar error differences increase with the interval j . When <math id="M78"><mrow><mi>E</mi><mrow><mo>(</mo><mrow><mrow><mo>‖</mo><mrow><msubsup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mi>R</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup></mrow><mo>‖</mo></mrow></mrow><mo>)</mo></mrow><mo>≈</mo><mi>E</mi><mrow><mo>(</mo><mrow><mrow><mo>‖</mo><mrow><msubsup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mi>R</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup><mo>−</mo><msubsup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mi>R</mi><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>−</mo><mi>n</mi></mrow><mo>)</mo></mrow></mrow></msubsup></mrow><mo>‖</mo></mrow></mrow><mo>)</mo></mrow></mrow></math> , the FSD inter-cycle correlation can be ignored after n cycles. The SW-based FSD inter-cycle correlation observation method uses the combined Sliced Wasserstein distance in three directions to measure the scalar fission source error <math id="M79"><mrow><msup><mi>ε</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup></mrow></math> , as expressed in Eq. (25). <math display="block" id="M80"><mrow><mtable><mtr><mtd><mrow><mtext>SW</mtext><mrow><mo>(</mo><mrow><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup><mo>,</mo><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup></mrow><mo>)</mo></mrow><mtext>   is used to measure </mtext><mrow><mo>‖</mo><mrow><msubsup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mi>R</mi><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup><mo>−</mo><msubsup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mi>R</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup></mrow><mo>‖</mo></mrow><mo>,</mo></mrow></mtd></mtr></mtable></mrow></math> where <math id="M81"><mrow><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup></mrow></math> is the fission source distribution in cycle i. The combined SW distance is calculated in Eq. (26). <math display="block" id="M82"><mrow><mtable><mtr><mtd><mrow><mtext>SW</mtext><mrow><mo>(</mo><mrow><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup><mo>,</mo><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup></mrow><mo>)</mo></mrow><mo>=</mo><msub><mrow><mtext>SW</mtext></mrow><mi>x</mi></msub><mrow><mo>(</mo><mrow><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup><mo>,</mo><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup></mrow><mo>)</mo></mrow><mo>×</mo><msub><mrow><mtext>SW</mtext></mrow><mi>y</mi></msub><mrow><mo>(</mo><mrow><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup><mo>,</mo><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup></mrow><mo>)</mo></mrow><mo>×</mo><msub><mrow><mtext>SW</mtext></mrow><mi>z</mi></msub><mrow><mo>(</mo><mrow><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup><mo>,</mo><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup></mrow><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mrow></math> We combined Eqs. 24 and 25, and the average SW distance between the FSD of adjacent cycles increased along with the adjacent cycle number, as expressed in Eq. (27). <math display="block" id="M83"><mrow><mtable><mtr><mtd><mrow><mo>∀</mo><mo> </mo><mi>i</mi><mo>></mo><mi>k</mi><mo>></mo><mi>j</mi><mo>,</mo><mtext>  </mtext><mi>E</mi><mrow><mo>(</mo><mrow><mtext>SW</mtext><mrow><mo>(</mo><mrow><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup><mo>,</mo><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup></mrow><mo>)</mo></mrow></mrow><mo>)</mo></mrow><mo>≥</mo><mi>E</mi><mrow><mo>(</mo><mrow><mtext>SW</mtext><mrow><mo>(</mo><mrow><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mi>k</mi><mo>)</mo></mrow></mrow></msup><mo>,</mo><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup></mrow><mo>)</mo></mrow></mrow><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mrow></math> The FSD inter-cycle correlation can be ignored when the average FSD SW distance fluctuates. Inter-cycle correlation length and Optimum batch length 2 We defined the interval cycle length <math id="M84"><mrow><msub><mi>n</mi><mn>0</mn></msub></mrow></math> where the average FSD SW distance starts to fluctuate as the FSD inter-cycle correlation length. The interval length <math id="M85"><mrow><msub><mi>n</mi><mn>0</mn></msub></mrow></math> conforms to Eq. (28). <math display="block" id="M86"><mrow><mtable><mtr><mtd><mrow><mo>∀</mo><mrow><mo>(</mo><mrow><mo> </mo><mi>j</mi><mo>≥</mo><msub><mi>n</mi><mn>0</mn></msub><mo> </mo></mrow><mo>)</mo></mrow><mo>,</mo><mo> </mo><mi>E</mi><mrow><mo>(</mo><mrow><mover accent="true"><mrow><mtext>SW</mtext><mrow><mo>(</mo><mrow><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup><mo>,</mo><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>−</mo><msub><mi>n</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow></mrow></msup></mrow><mo>)</mo></mrow></mrow><mo stretchy="true">¯</mo></mover></mrow><mo>)</mo></mrow><mo>≅</mo><mi>E</mi><mrow><mo>(</mo><mrow><mover accent="true"><mrow><mtext>SW</mtext><mrow><mo>(</mo><mrow><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup><mo>,</mo><msup><mstyle mathvariant="bold-italic"><mi>s</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>−</mo><mi>j</mi></mrow><mo>)</mo></mrow></mrow></msup></mrow><mo>)</mo></mrow></mrow><mo stretchy="true">¯</mo></mover></mrow><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mrow></math> According to Eq. (25), the FSD SW distance was used to measure the scalar fission source error. Thus, Eq. (29) can be expressed as <math display="block" id="M87"><mrow><mtable><mtr><mtd><mrow><mo>∀</mo><mrow><mo>(</mo><mrow><mo> </mo><mi>j</mi><mo>≥</mo><msub><mi>n</mi><mn>0</mn></msub><mo> </mo></mrow><mo>)</mo></mrow><mo>,</mo><mtext>​</mtext><mo> </mo><mi>E</mi><mrow><mo>(</mo><mrow><mover accent="true"><mrow><msup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup><mo>−</mo><msup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>−</mo><msub><mi>n</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow></mrow></msup></mrow><mo stretchy="true">¯</mo></mover></mrow><mo>)</mo></mrow><mo>≅</mo><mi>E</mi><mrow><mo>(</mo><mrow><mover accent="true"><mrow><msup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msup><mo>−</mo><msup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>−</mo><mi>j</mi></mrow><mo>)</mo></mrow></mrow></msup></mrow><mo stretchy="true">¯</mo></mover></mrow><mo>)</mo></mrow><mo>.</mo></mrow></mtd></mtr></mtable></mrow></math> If the interval cycle length j is greater than <math id="M88"><mrow><msub><mi>n</mi><mn>0</mn></msub></mrow></math> , the correlation between the FSDs of cycle i and cycle <math id="M89"><mrow><mi>i</mi><mo>−</mo><mi>j</mi></mrow></math> declines statistically to negligible. According to the Law of Large Numbers, we can obtain <math display="block" id="M90"><mrow><mtable><mtr><mtd><mrow><mo>∀</mo><mrow><mo>(</mo><mrow><mi>j</mi><mo>≥</mo><msub><mi>n</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><mo>,</mo><mtext> </mtext><mi>E</mi><mrow><mo>(</mo><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>0</mn></mrow><mi>l</mi></munderover><mrow><mrow><mo>‖</mo><mrow><msup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>+</mo><mi>k</mi></mrow><mo>)</mo></mrow></mrow></msup><mo>−</mo><msup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>+</mo><mi>k</mi><mo>−</mo><msub><mi>n</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow></mrow></msup></mrow><mo>‖</mo></mrow></mrow></mstyle></mrow><mo>)</mo></mrow><mo>≅</mo><mi>E</mi><mrow><mo>(</mo><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>0</mn></mrow><mi>l</mi></munderover><mrow><mrow><mo>‖</mo><mrow><msup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>+</mo><mi>k</mi></mrow><mo>)</mo></mrow></mrow></msup><mo>−</mo><msup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>+</mo><mi>k</mi><mo>−</mo><mi>j</mi></mrow><mo>)</mo></mrow></mrow></msup></mrow><mo>‖</mo></mrow></mrow></mstyle></mrow><mo>)</mo></mrow><mo>,</mo></mrow></mtd></mtr></mtable></mrow></math> where l is group length. We assumed the cumulated scalar fission error of l cycles to be uncorrelated after <math id="M91"><mrow><msub><mi>n</mi><mn>0</mn></msub></mrow></math> cycles. Specifically, when group length l is set as the FSD inter-cycle correlation length <math id="M92"><mrow><msub><mi>n</mi><mn>0</mn></msub></mrow></math> , Eq. (31) can be expressed as <math display="block" id="M93"><mrow><mtable><mtr><mtd><mrow><mo>∀</mo><mrow><mo>(</mo><mrow><mi>j</mi><mo>≥</mo><msub><mi>n</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><mo>,</mo><mtext> </mtext><mi>E</mi><mrow><mo>(</mo><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>0</mn></mrow><mrow><msub><mi>n</mi><mn>0</mn></msub></mrow></munderover><mrow><mrow><mo>‖</mo><mrow><msup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>+</mo><mi>k</mi></mrow><mo>)</mo></mrow></mrow></msup><mo>−</mo><msup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>+</mo><mi>k</mi><mo>−</mo><msub><mi>n</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow></mrow></msup></mrow><mo>‖</mo></mrow></mrow></mstyle></mrow><mo>)</mo></mrow><mo>≅</mo><mi>E</mi><mrow><mo>(</mo><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>0</mn></mrow><mrow><msub><mi>n</mi><mn>0</mn></msub></mrow></munderover><mrow><mrow><mo>‖</mo><mrow><msup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>+</mo><mi>k</mi></mrow><mo>)</mo></mrow></mrow></msup><mo>−</mo><msup><mstyle mathvariant="bold-italic"><mi>e</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>i</mi><mo>+</mo><mi>k</mi><mo>−</mo><mi>j</mi></mrow><mo>)</mo></mrow></mrow></msup></mrow><mo>‖</mo></mrow></mrow></mstyle></mrow><mo>)</mo></mrow><mo>.</mo></mrow></mtd></mtr></mtable></mrow></math> The cumulative scalar fission error of <math id="M94"><mrow><msub><mi>n</mi><mn>0</mn></msub></mrow></math> cycles is statistically uncorrelated after <math id="M95"><mrow><msub><mi>n</mi><mn>0</mn></msub></mrow></math> cycles. The flux results <math id="M96"><mrow><mrow><mo>{</mo><mrow><msubsup><mi>Q</mi><mtext>b</mtext><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup></mrow><mo>}</mo></mrow></mrow></math> are then statistically uncorrelated using <math id="M97"><mrow><msub><mi>n</mi><mn>0</mn></msub></mrow></math> as the batch length. The variance underestimation problem can be solved consequently, as expressed in Eq. (32). <math display="block" id="M98"><mrow><mtable><mtr><mtd><mrow><msubsup><mi>σ</mi><mrow><mtext>br</mtext></mrow><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow><mo>−</mo><msubsup><mi>σ</mi><mrow><mtext>ba</mtext></mrow><mn>2</mn></msubsup><mrow><mo>(</mo><mover accent="true"><mi>Q</mi><mo>¯</mo></mover><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>n</mi><mtext>b</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>n</mi><mtext>b</mtext></msub><mo>−</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></mfrac><mstyle displaystyle="true"><munder><mo>∑</mo><mi>i</mi></munder><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>j</mi><mo>≠</mo><mi>i</mi></mrow></munder><mrow><mtext>cov</mtext></mrow></mstyle></mrow></mstyle><mrow><mo>(</mo><mrow><msubsup><mi>Q</mi><mtext>b</mtext><mrow><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></msubsup><mo>,</mo><msubsup><mi>Q</mi><mtext>b</mtext><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup></mrow><mo>)</mo></mrow><mo>≅</mo><mn>0</mn></mrow></mtd></mtr></mtable></mrow></math> In conclusion, the FSD inter-cycle correlation length is the optimum batch length. Semi-quantitative Variance Underestimation Elimination (SeVUE) method 1 We used the SeVUE method to observe the FSD inter-cycle correlation length and eliminate the variance underestimation problem. Sections 3.1 and 3.2 describe the workflow of the SeVUE method and the implementation scheme, respectively ( Fig. 2 Fig. 2 ). Workflow 2 The SeVUE method comprises two procedures: the Estimation procedure and the Elimination procedure. In the Estimation procedure, the SeVUE method requires an MC calculation to estimate the FSD inter-cycle correlation length of each model. This calculation should be in the criticality mode with no tallies. The number of particles per cycle can be set to a relatively small value because the FSD inter-cycle correlation is almost independent of it. For example, setting 10,000 neutrons per cycle is sufficient for a full-core model. By contrast, the number of active cycles should be relatively large. Thousands of active cycles are recommended for models with a high dominance ratio. The SW distances between FSDs can be calculated using external scripts ( Fig. 4 Fig. 4 ) or embedded codes. At the end of this procedure, the relationship between the average FSD SW distances and the interval cycle length can be plotted ( Fig. 3 Fig. 3 ). The FSD inter-cycle correlation length is the point at which the average FSD SW distance starts to fluctuate. In the Elimination procedure, the batch method was used to eliminate variance underestimation. According to Eq. (19-22), the FSD inter-cycle correlation length was set as the batch length. The number of active cycles should be at least 10 times the FSD inter-cycle correlation length. The calculated variances of the tally results were close to the actual variances. Implementation 2 The SeVUE method was implemented using the MC code RMC [ 20 20 ]. RMC is a neutron photon-electron transport code developed by the Reactor Engineering Analysis Lab (REAL) of Tsinghua University. RMC has been applied for nuclear reactor criticality, burnup, shielding, and neutronics-thermal coupling calculations [ 25 25 - 28 28 25-28 ]. In the Estimation procedure, an external script was used to calculate the SW distances between FSDs of different cycles. In the RMC code, the fission source distribution is the output for each cycle. The external script can run independently and be easily coupled with other MC codes. The coupled workflow is illustrated in Fig. 3 Fig. 3 . In the Elimination procedure, the batch method is implemented in the RMC code. Tally results were processed at the end of each batch. Validation and discussion 1 To validate the SeVUE method and study its characteristics, we calculated the flux distribution of the sphere array model from the OECD/NEA source convergence benchmark and the power distribution of the BEAVRS full-core model from the BEAVRS benchmark [ 29 29 - 30 30 29-30 ]. The platform used was an AMD Ryzen 9 3990X. The base RMC version was V3.5.0. Flux distribution of sphere array model 2 The sphere array model was obtained from the OECD/NEA source convergence benchmark, with a dominance ratio of approximately 0.886. It comprised of 25 fissile metal spheres in a <math id="M99"><mrow><mn>5</mn><mo>×</mo><mn>5</mn><mo>×</mo><mn>1</mn></mrow></math> array ( Fig. 5 Fig. 5 ). The metal spheres were filled with highly enriched uranium. The center sphere was larger than the other eight spheres, with a radius of 10 cm. The flux distribution in this model was of interest in this study. The number of particles per cycle of 5,000, 50,000, and 5,000,000 was analyzed to observe the FSD inter-cycle correlation lengths. Each calculation had 1,000 inactive and 20,000 active cycles. Reference flux distribution and variance underestimation phenomenon 3 The reference flux distribution and relative standard deviation (RSD) distributions were obtained using 100 independent calculations. The fluxes in each sphere were tallied. The reference flux results and RSDs are shown in Fig. 6 Fig. 6 . According to the Law of large numbers, the RSD of the sample variance was approximately <math id="M100"><mrow><msqrt><mrow><mfrac><mn>2</mn><mrow><mi>N</mi><mo>−</mo><mn>1</mn></mrow></mfrac></mrow></msqrt></mrow></math> (where N is the number of samples). Thus, the reference standard deviation has a <math id="M101"><mrow><mo>±</mo><mn>10</mn><mtext>%</mtext></mrow></math> fluctuation. The variance underestimation ratio distribution in Fig. 7 Fig. 7 was obtained by comparing the relative standard deviations in a single simulation with the reference results. The Variance Underestimation Ratio ( VUR ) is defined in Eq. (33). <math display="block" id="M102"><mrow><mtable><mtr><mtd><mrow><mtext>VUR</mtext><mo>=</mo><mfrac><mrow><msub><mi>σ</mi><mtext>r</mtext></msub><mo>−</mo><msub><mi>σ</mi><mtext>a</mtext></msub></mrow><mrow><msub><mi>σ</mi><mtext>r</mtext></msub></mrow></mfrac><mo>×</mo><mn>100</mn><mi>%</mi><mo>,</mo></mrow></mtd></mtr></mtable></mrow></math> where <math id="M103"><mrow><msub><mi>σ</mi><mi>a</mi></msub></mrow></math> is the calculated standard deviation, and <math id="M104"><mrow><msub><mi>σ</mi><mi>r</mi></msub></mrow></math> is the true standard deviation. As shown in Fig. 7 Fig. 7 , the VURs in each sphere varied from 24% to 47%. The average VUR was 37.3%. Estimation procedure 3 The number of particles per cycle of 5,000, 50,000, and 5,000,000 was analyzed to observe the FSD inter-cycle correlation lengths in this procedure. Fig. 8 Fig. 8 shows that the FSD average SW distance evolves with the interval cycle length. A visual check of the curves showed that the FSD inter-cycle correlation lengths were approximately 50 for all three batch sizes. Thus, the FSD inter-cycle correlation was almost independent of batch size. Semi-quantitative information can be obtained from Fig. 8 Fig. 8 . The initial FSD average SW distance was approximately 40% of the stationary SW distance. This "similarity coefficient" can be further used to measure the variance underestimation ratio of the total flux in this model. Elimination procedure 3 We verified Eqs. (28-32) by setting different batch lengths. The batch lengths are listed in Table. 1 Table. 1 . Notably, in the SeVUE method, only the case with <math id="M105"><mrow><msub><mi>l</mi><mtext>b</mtext></msub><mo>=</mo><mn>50</mn></mrow></math> should be calculated. 1 1 Serial number 1 1 Batch length/cycle 1 1 Number of batches 1 1 Number of active cycles 1 1 1 1 1 1 1 1 20,000 1 1 20,000 1 1 2 1 1 2 1 1 10,000 1 1 20,000 1 1 3 1 1 4 1 1 5000 1 1 20,000 1 1 4 1 1 5 1 1 4000 1 1 20,000 1 1 5 1 1 8 1 1 2500 1 1 20,000 1 1 6 1 1 10 1 1 2000 1 1 20,000 1 1 7 1 1 16 1 1 1250 1 1 20,000 1 1 8 1 1 20 1 1 1000 1 1 20,000 1 1 9 1 1 25 1 1 800 1 1 20,000 1 1 10 1 1 32 1 1 625 1 1 20,000 1 1 11 1 1 40 1 1 500 1 1 20,000 1 1 12 1 1 50 1 1 400 1 1 20,000 Figure 9 Figure 9 shows that the average VURs evolve with the batch length. The average VURs of different batch sizes are approximately 37% when the batch length is one cycle. With the batch length increases, the average VURs gradually declines. When batch length is 50, the average VURs are approximately <math id="M106"><mrow><mo>±</mo><mn>1</mn><mi>%</mi></mrow></math> . Thus, 50 is the proper batch length for this model. This value is in satisfactory agreement with the FSD inter-cycle correlation length estimated in the former procedure. The results of a single calculation for <math id="M107"><mrow><msub><mi>l</mi><mtext>b</mtext></msub><mo>=</mo><mn>50</mn></mrow></math> are shown in Fig. 10 Fig. 10 . The VUR of each sphere varied from -5.0% to 2.6%, with an average VUR of -1.9%. Considering the fluctuations in the reference standard deviations, these differences were acceptable. Power distribution of BEAVRS 3D full-core model 2 The BEAVRS full-core model is a detailed reactor core model used to test indicators in real-world PWRs. The geometric configuration is illustrated in Fig. 11 Fig. 11 . It consists of 193 fuel assemblies in a 15 <math id="M109"><mo>×</mo></math> 15 <math id="M110"><mo>×</mo></math> 1 array with a dominance ratio of 0.989. The power distribution in this model was of interest in this study. The number of particles per 50,000 cycles was analyzed to observe the FSD inter-cycle correlation length. Each calculation had 1,000 inactive and 20,000 active cycles. The total number of neutron histories in the active cycles was 1 billion, which is practical for actual calculations. Reference power distribution and variance underestimation phenomenon 3 The reference power distribution and standard deviation distribution were obtained using 30 independent calculations, similar to the sphere array model. The power of each assembly was tallied. The reference results and RSDs are shown in Fig. 12 Fig. 12 . The variance underestimation ratio distribution in Fig. 13 Fig. 13 was obtained by comparing the RSDs in a single MC simulation with the reference results. The VURs in each assembly varied from 72% to 91%. The VURs of the assemblies close to the center were lower than the average. By contrast, the outer assemblies had a VUR above 90%. The average VUR was 87.5%. Estimation procedure 3 A batch size of 50,000 was analyzed to observe the FSD inter-cycle correlation length. Fig. 14 Fig. 14 shows that the FSD average SW distance evolves with the interval cycle length. A visual check of the curves showed that the FSD inter-cycle correlation length of this model was approximately 500. The initial FSD average SW distance was approximately 1/10 of the stationary SW distance. Thus, the variance underestimation problem is more severe in this model than in the sphere array model. This conclusion is consistent with the results presented in Figs. 7 Figs. 7 and 13 13 . Elimination procedure 3 We set different batch lengths to verify that the FSD inter-cycle correlation length was the optimum batch length. The batch lengths are listed in Table 2 Table 2 . In the SeVUE method, only the case where <math id="M111"><mrow><msub><mi>l</mi><mtext>b</mtext></msub><mo>=</mo><mn>500</mn></mrow></math> should be calculated. 1 1 Serial number 1 1 Batch length/cycle 1 1 Number of batches 1 1 Number of active cycles 1 1 1 1 1 1 1 1 20,000 1 1 20,000 1 1 2 1 1 2 1 1 10,000 1 1 20,000 1 1 3 1 1 10 1 1 2000 1 1 20,000 1 1 4 1 1 20 1 1 1000 1 1 20,000 1 1 5 1 1 40 1 1 500 1 1 20,000 1 1 6 1 1 80 1 1 250 1 1 20,000 1 1 7 1 1 100 1 1 200 1 1 20,000 1 1 8 1 1 125 1 1 160 1 1 20,000 1 1 9 1 1 200 1 1 100 1 1 20,000 1 1 10 1 1 250 1 1 80 1 1 20,000 1 1 11 1 1 400 1 1 50 1 1 20,000 1 1 12 1 1 500 1 1 40 1 1 20,000 1 1 13 1 1 625 1 1 32 1 1 20,000 1 1 14 1 1 800 1 1 25 1 1 20,000 1 1 15 1 1 1000 1 1 20 1 1 20,000 Figure 15 Figure 15 shows that the average VURs evolve with the batch length. The average VUR was approximately 87% when the batch length was one cycle. It gradually declined with the batch length increases. When the batch length exceeded 500, the average VURs were within <math id="M112"><mrow><mo>±</mo><mn>5</mn><mi>%</mi><mo>,</mo></mrow></math> and the curve fluctuated. Thus, 500 is the optimum batch length for this model. The results of the single MC simulation with <math id="M114"><mrow><msub><mi>l</mi><mtext>b</mtext></msub><mo>=</mo><mn>500</mn></mrow></math> are shown in Fig. 16 Fig. 16 . The VUR of each assembly varies from -10.5% to 13.4%, with an average VUR of -1.3%. The differences were acceptable, considering the fluctuations in the reference standard deviations and batch method results. Discussion 2 Relatively large tallies were used to exemplify the efficiency of the SeVUE method, and this method effectively eliminated the variance underestimation problem in MC criticality simulations. The FSD inter-cycle correlation length was validated to be independent of the batch size. The ratio between the initial FSD average SW distance and stationary SW distance in the Estimation procedure can be a potential reference for the variance underestimation of large local tallies. The penalty is that an additional Estimation calculation must be conducted to estimate the FSD inter-cycle correlation length of every model. In addition, the FSD inter-cycle correlation length may be hundreds of cycles in the reactor full-core model. Thousands of active cycles and billions of neutron histories are required to eliminate the variance underestimation problem. Moreover, the standard deviation of the tally results may have <math id="M115"><mrow><mo>±</mo><mn>10</mn><mtext>%</mtext></mrow></math> fluctuations in the VUR when a proper batch length is set. Conclusion 1 This study provides a mesh-free SeVUE method to obtain a credible confidence interval for the tallied results. This method comprises two procedures: Estimation and Elimination. The SW-based variance underestimation estimation method was used to estimate the FSD inter-cycle correlation length. The batch method was used in the Elimination procedure. The FSD inter-cycle correlation length is the interval cycle length <math id="M116"><mrow><msub><mi>n</mi><mn>0</mn></msub></mrow></math> where the average FSD SW distance starts to fluctuate. It was proved to be the optimum batch length to eliminate the variance underestimation problem. We exemplified this method by using the OECD sphere array model and 3D PWR BEAVRS model. The sphere flux and assembly power were analyzed. The average VURs of the sphere flux and assembly power decreased from 37% and 87% to within <math id="M117"><mo>±</mo></math> 5%, respectively. The calculated variances were successfully compared with variance estimates from independent simulations. One-step variance underestimation elimination methods are under development. Reference D.B. MacMillan ,  Monte Carlo confidence limits for iterated-source calculations. Nucl . Sci. Eng.  50 ,   73 - 75 ( 1973 ).  doi: 10.13182/NSE73-A22590 http://doi.org/10.13182/NSE73-A22590  Monte Carlo confidence limits for iterated-source calculations. Nucl P. Turinski , Few-group neutron diffusion equation solver utilizing the nodal expansion method for eigenvalue, adjoint, fixed-source steadystate and transient problems. North Carolina State Univ. Report . ( 1994 ). Few-group neutron diffusion equation solver utilizing the nodal expansion method for eigenvalue, adjoint, fixed-source steadystate and transient problems. North Carolina State Univ. Report H. J. Shim , C. H. Kim , Real variance estimation using an intercycle fission source correlation for Monte Carlo eigenvalue calculations . Nuclear science and engineering , 162 ( 1 ), 98 - 108 ( 2009 ). doi: 10.13182/NSE09-2 http://doi.org/10.13182/NSE09-2 Real variance estimation using an intercycle fission source correlation for Monte Carlo eigenvalue calculations H. J. Park , H. C. Lee , H. J. Shim , J. Y. Cho , Real variance analysis of Monte Carlo eigenvalue calculation by McCARD for BEAVRS benchmark . Annals of Nuclear Energy 90 , 205 - 211 ( 2016 ). doi: 10.1016/j.anucene.2015.12.009 http://doi.org/10.1016/j.anucene.2015.12.009 Real variance analysis of Monte Carlo eigenvalue calculation by McCARD for BEAVRS benchmark L.J. Pan , R.H. Wang , S. Jiang , Monte Carlo fission matrix acceleration method with adaptive mesh . Nucl. Sci. Tech . 26 , 050603 ( 2015 ). doi: 10.13538/j.1001-8042/nst.26.050603 http://doi.org/10.13538/j.1001-8042/nst.26.050603 Monte Carlo fission matrix acceleration method with adaptive mesh Q.Q. Pan , T.F. Zhang , X.J. Liu , et al. , SP3-coupled global variance reduction method based on RMC code . Nucl. Sci. Tech . 32 , 122 ( 2021 ). doi: 10.1007/s41365-021-00973-0 http://doi.org/10.1007/s41365-021-00973-0 SP3-coupled global variance reduction method based on RMC code X.-L. Pan , J.-Q. Wang , R. Yuan et al. , Biasing transition rate method based on direct MC simulation for probabilistic safety assessment . Nucl. Sci. Tech . 28 , 91 ( 2017 ). doi: 10.1007/s41365-017-0255-2 http://doi.org/10.1007/s41365-017-0255-2 Biasing transition rate method based on direct MC simulation for probabilistic safety assessment X.-C. Nie , J. Li , S.-L. Liu et al. , Global variance reduction method for global Monte Carlo particle transport simulations of CFETR . Nucl. Sci. Tech . 28 , 115 ( 2017 ). doi: 10.1007/s41365-017-0270-3 http://doi.org/10.1007/s41365-017-0270-3 Global variance reduction method for global Monte Carlo particle transport simulations of CFETR Z.G. Li , K. Wang , Y.C. Guo et al. , Forced propagation method for Monte Carlo fission source convergence acceleration in the RMC . Nucl. Sci. Tech . 32 , 27 ( 2021 ). doi: 10.1007/s41365-021-00868-0 http://doi.org/10.1007/s41365-021-00868-0 Forced propagation method for Monte Carlo fission source convergence acceleration in the RMC R.J. Brissenden , A.R. Garlick , Biases in the estimation of Keff and its error by Monte Carlo methods . Ann. Nucl. Energy 13 , 63 - 83 ( 1986 ). doi: 10.1016/0306-4549(86)90095-2 http://doi.org/10.1016/0306-4549(86)90095-2 Biases in the estimation of Keff and its error by Monte Carlo methods B. Herman , B. Forget , K. Smith , Analysis of tally correlation in large light water reactors . PHYSOR 2014 , Kyoto, Japan Analysis of tally correlation in large light water reactors J. Park , P. Zhang , H.C. Lee et al. , Performance evaluation of CMFD on inter-cycle correlation reduction of Monte Carlo simulation . Computer Phys. Communications 235 , 111 - 119 ( 2019 ). doi: 10.1016/j.cpc.2018.09.014 http://doi.org/10.1016/j.cpc.2018.09.014 Performance evaluation of CMFD on inter-cycle correlation reduction of Monte Carlo simulation E. M. Gelbard , R. Prael , Computation of standard deviations in Eigenvalue calculations . Prog. Nucl. Energy 24 ( 1 ), 237 - 241 ( 1990 ). doi: 10.1016/0149-1970(90)90041-3 http://doi.org/10.1016/0149-1970(90)90041-3 Computation of standard deviations in Eigenvalue calculations L. Deng , G. Li , B.Y. Zhang et al. , A high fidelity general purpose 3-D Monte Carlo particle transport program JMCT3.0 . Nucl. Sci. Tech . 33 , 108 ( 2022 ). doi: 10.1007/s41365-022-01092-0 http://doi.org/10.1007/s41365-022-01092-0 A high fidelity general purpose 3-D Monte Carlo particle transport program JMCT3.0 T. Ueki , Standard deviation of local tallies in global monte carlo calculation of nuclear reactor core . Journal of nuclear science and technology 47 ( 8 ), 739 - 753 ( 2010 ). doi: 10.3327/jnst.47.739 http://doi.org/10.3327/jnst.47.739 Standard deviation of local tallies in global monte carlo calculation of nuclear reactor core T.M. Sutton , Application of a discretized phase space approach to the analysis of Monte Carlo uncertainties . In: ANS MC2015 – Joint International Conference Mathematics and Computation (M&C), Supercomputing in Nuclear Applications (SNA) and the Monte Carlo (MC) Method, no. Mc , Nashville, TN , 2015 , pp. 1 – 17 . Application of a discretized phase space approach to the analysis of Monte Carlo uncertainties T. Ueki , Fractal dimension analysis for run length diagnosis of Monte Carlo criticality calculation . Journal of Nuclear Science and Technology 53 ( 3 ), 312 - 322 ( 2016 ). doi: 10.1080/00223131.2015.1039620 http://doi.org/10.1080/00223131.2015.1039620 Fractal dimension analysis for run length diagnosis of Monte Carlo criticality calculation J. Miao , B. Forget , K. Smith , Predicting correlation coefficients for Monte Carlo eigenvalue simulations with multitype branching process . Annals of Nuclear Energy 112 , 307 - 321 ( 2018 ). doi: 10.1016/j.anucene.2017.10.014 http://doi.org/10.1016/j.anucene.2017.10.014 Predicting correlation coefficients for Monte Carlo eigenvalue simulations with multitype branching process J. Miao , B. Forget , K. Smith , Correlation diagnosis method for heterogeneous Monte Carlo eigenvalue simulations based on a diffusion approximation . Annals of Nuclear Energy 130 , 301 - 318 ( 2019 ). doi: 10.1016/j.anucene.2019.02.048 http://doi.org/10.1016/j.anucene.2019.02.048 Correlation diagnosis method for heterogeneous Monte Carlo eigenvalue simulations based on a diffusion approximation K. Wang , Z.G. Li , D. She et al. , RMC–a monte carlo code for reactor core analysis . Ann. Nucl. Energy 82 , 121 - 129 ( 2015 ). doi: 10.1016/j.anucene.2014.08.048 http://doi.org/10.1016/j.anucene.2014.08.048 RMC–a monte carlo code for reactor core analysis M. Arjovsky , S. Chintala , L. Bottou , Wasserstein gan . 2017 . arXiv:1701.07875 Wasserstein gan X. Guo , Z. Li , S. Huang , et al. , Convergence diagnostics for Monte Carlo fission source distributions using the Wasserstein distance measure . Nuclear Engineering and Design 389 , 111675 ( 2022 ) Convergence diagnostics for Monte Carlo fission source distributions using the Wasserstein distance measure P. Shen , S. Jiang , Y. Hu et al. , Research on the correlations of fission source distribution in the Monte Carlo algorithm using mathematical tools . Proceedings of 29th International Conference on Nuclear Engineering(ICONE 29) ( 2022 ) Research on the correlations of fission source distribution in the Monte Carlo algorithm using mathematical tools P. Shen , X. Huo , S. Huang et al. , Capturing and utilizing the random feature in Monte Carlo fission source distributions . Ann. Nucl. Energy 180 , 109468 ( 2022 ). Capturing and utilizing the random feature in Monte Carlo fission source distributions Q. Pan , J. Rao , S. Huang et al. , Improved adaptive variance reduction algorithm based on RMC code for deep penetration problems . Ann. Nucl. Energy 137 , 107113 ( 2020 ). doi: 10.1016/j.anucene.2019.107113 http://doi.org/10.1016/j.anucene.2019.107113 Improved adaptive variance reduction algorithm based on RMC code for deep penetration problems P. Shen , X. Guo , K. Li et al. , Research on global neighbor list method in Monte Carlo code RMC . Ann. Nucl. Energy 167 , 108861 ( 2022 ). doi: 10.1016/j.anucene.2021.108861 http://doi.org/10.1016/j.anucene.2021.108861 Research on global neighbor list method in Monte Carlo code RMC K. Wang , S. Liu , Z. Li et al. , Analysis of beavrs two-cycle benchmark using rmc based on full core detailed model . Prog. Nucl. Energy 98 , 301 - 312 ( 2017 ). doi: 10.1016/j.pnucene.2017.04.009 http://doi.org/10.1016/j.pnucene.2017.04.009 Analysis of beavrs two-cycle benchmark using rmc based on full core detailed model Y. Guo , Z. Li , S. Huang et al. , A new neutronics-thermal-mechanics multi-physics coupling method for heat pipe cooled reactor based on RMC and OpenFOAM . Prog. Nucl. Energy 139 , 103842 ( 2021 ). doi: 10.1016/j.pnucene.2021.103842 http://doi.org/10.1016/j.pnucene.2021.103842 A new neutronics-thermal-mechanics multi-physics coupling method for heat pipe cooled reactor based on RMC and OpenFOAM R. Blomquist , A. Nouri , The OECD/NEA source convergence benchmark program . Transactions of the American Nuclear Society  87 , 143 - 145 ( 2002 ). The OECD/NEA source convergence benchmark program N. Horelik , B. Herman , B. Forget et al. , Benchmark for evaluation and validation of reactor simulations (beavrs) , v 1.0.1 , in: Proc. Int. Conf. Mathematics and Computational Methods Applied to Nuc. Sci. & Eng , pp. 63 – 68 ( 2013 ). Benchmark for evaluation and validation of reactor simulations (beavrs) All authors contributed to the study conception and design. Material preparation, data collection, and analysis were performed as described by Peng-Fei Shen. Programming and tests were strongly supported by Ze-Guang Li, Hai-Feng Yang, Zeng Shao and Peng Zhang. Xiaodong Huo and Kan Wang supervised the study. The first draft of the manuscript was written by Peng-Fei Shen, and all the authors commented on the previous versions of the manuscript. All authors read and approved the final manuscript. 10.1007/s41365-022-01156-1 This work was supported by China Nuclear Power Engineering Co., Ltd. Scientific Research Project (No. KY22104) and the fellowship of China Postdoctoral Science Foundation (No. 2022M721793). 2022-08-03 2022-11-03 2022-11-04 2023-01-30 2023-02-20 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2023 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-022-01156-1 http://dx.doi.org/10.1007/s41365-022-01156-1 . Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law. 2023 Nuclear Science and Techniques 01 34