{'fraction_non_alphanumeric': 0.036677164144711574, 'fraction_numerical': 0.007508257162608495, 'mean_word_length': 4.5990753682399745, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analysis of Population Control Techniques for Time-Dependent and Eigenvalue Monte Carlo Neutron Transport Calculations

An extensive study of population control techniques (PCTs) for time-dependent and eigenvalue Monte Carlo (MC) neutron transport calculations is presented. We define PCT as a technique that takes a censused population and returns a controlled, unbiased population. A new perspective based on an abstraction of particle census and population control is explored, paving the way to improved understanding and application of the concepts. Five distinct PCTs identified from the literature are reviewed: Simple Sampling (SS), Splitting-Roulette (SR), Combing (CO), modified Combing (COX), and Duplicate-Discard (DD). A theoretical analysis of how much uncertainty is introduced to a population by each PCT is presented. Parallel algorithms for the PCTs applicable for both time-dependent and eigenvalue MC simulations are proposed. The relative performances of the PCTs based on runtime and tally mean error or standard deviation are assessed by solving time-dependent and eigenvalue test problems. It is found that SR and CO are equally the most performant techniques, closely followed by DD.


I. Introduction
The Monte Carlo (MC) method is indispensable in neutron transport calculations due to its ability to perform high-fidelity, continuous-energy transport simulations with minimal approximation. MC, however, suffers from stochastic uncertainties requiring an expensive computation of a large number of neutron source samples or histories.
Nevertheless, thanks to the advancement of high-performance parallel computing, the inherently parallel features of MC can be effectively exploited to a very large extentwhich can significantly reduce run time to solution, particularly for the computationally expensive time-dependent neutron transport simulations [1,2,3,4,5,6,7,8,9,10].
During a time-dependent MC simulation, particle population size can monotonically grow or decay depending on the criticality of the system. This monotonic evolution of population makes time-dependent MC simulation particularly challenging in two different ways. First, in a supercritical system, particle population size can quickly grow beyond the limited computational resources. Additionally, some MC implementations and variance reduction techniques-such as the precursor forced decay technique in [2] and time-dependent adaptation of the hybrid source iteration methods in [11,12]may promote monotonic population growth, which raises the same issue on the limited computational memory. Second, in a subcritical system without a significant persisting external source-such as in pulsed-reactor and shut-down experiments-particle population size can quickly decay to zero, which leads to a lack of samples and yields statistically noisy tally results at later times of the simulation.
One typically uses a Population Control Technique (PCT) to address the monotonic population growth and collapse issues discussed above. PCT essentially controls the size of a particle population to be near a desired value while preserving certain statistical expectations to achieve unbiased MC simulation. In the implementation of PCT, time census is employed to limit the population growth/collapse. The census introduces a time boundary that stops particles whenever they are about to cross it.
When all particles have already hit the time boundary, the time census is completed, and PCT can be performed on the censused particles. More recent applications of PCT include the use of random particle duplication or discard [3] in Serpent 2 [13], splitting and Russian-Roulette technique [4] in MCATK [14], particle combing technique [15] in TRIPOLI-4 [16,17] and GUARDYAN [8], and a modified combing technique which is most recently introduced in [10].
An innovative approach to performing time-dependent MC is proposed by [9]. The central idea is to re-purpose the generally available k-eigenvalue MC simulation infrastructure to perform time-dependent simulations. This approach works because there is a built-in population control in k-eigenvalue MC simulation. Besides the introduction of the 1/k factor on the fission operator, which is essential in achieving a steady-state configuration, simple sampling is typically performed to ensure that a certain number of particles are sampled from the fission bank and then used as the particle source for the simulation of the next fission generation. Observing the significance of that connection between the k-eigenvalue and time-dependent MC simulations offers an improved understanding of PCT. Such a study has been done to an extent by Cullen et al.
in [18]. Nevertheless, one can take advantage of this connection further by exploring potential benefits from and for both of the simulation modes.
Despite the multiple distinct PCTs proposed in the literature [2,3,4,9,10], documented studies in characterizing and assessing relative performances of all the identified PCTs are still very limited. A more recent effort found in [4] specifically compares the splitting and Russian-Roulette technique [4] to the particle combing technique [15]-hereafter referred to as Splitting-Roulette (SR) and Combing (CO), respectively.
Sweezy et al. [4] propose a normalized SR as an alternative to CO which may suffer from unwanted behavior due to possible correlations in the particle order. On the other hand, Faucher et al. [17] and Legrady et al. [19] prefer the use of CO instead of SR due to the inherent bias in the normalized SR [4] and suggest that the unwanted behavior of CO is unlikely to occur in practice. This support for CO, or, if you will, ctenophilia, is strengthened by an assertion [19] stating that the SR technique described in [4] is at least 2-3 times less efficient than CO. This assertion is based on a comparative study of CO and a "Russian roulette and splitting" technique reported in [20]. However, the Russian roulette and splitting technique in [20] seems to be different from the SR technique described in [4]; and for the record, the study [20] does not claim that the technique compared refers to [4]. Consistent implementation and fair comparison of the two techniques, SR [4] and CO [15], would shed light on their actual relative performances.
In this paper, we present an extensive study on PCT. In Sec. II, we start by making an abstraction of related concepts-i.e., particle census and population controlfollowed by reviewing PCTs identified from the literature. In Sec. III, we perform an analysis to reveal the theoretical uncertainty introduced by each of the PCTs, which directly affects the performance of the technique; these theoretical uncertainties are then verified numerically. Sec. IV presents a parallel PCT algorithm that exploits the abstraction established in Sec. II and adapts the nearest-neighbor parallel fission bank algorithm proposed in [21]. In Secs. V and VI, we implement and test the PCTs on time-dependent and eigenvalue MC neutron transport problems, respectively. Finally, Sec. VII summarizes the takeaways of the study. It is worth mentioning that while this paper is focused on PCT application for MC neutron transport, discussions and analyses presented in this paper are also applicable to PCT application in other transport simulations, such as the Implicit Monte Carlo thermal radiative transfer [22].

II. Population Control Technique (PCT)
Population control can be loosely defined as any MC technique that involves altering the number of particles being simulated; this includes many variance reduction techniques (e.g., cell importance and weight window) and even the introduction of 1/k factor in eigenvalue simulations [18]. However, in this paper, we specifically define population control as a technique that controls populations of censused particles.
In this section, we present an abstraction of particle census and population control (their definitions and how they are characterized) and then discuss distinct techniques identified from the literature.

II.A. Particle Census
Census is a process where we (1) stop particles, (2) remove them from the current simulation, and then (3) store them into a census bank. Census can be performed at arbitrary steps during simulation; however, there are several triggering events that physically make sense to perform the census.
Perhaps the most obvious one is time-grid crossing. In this time census, we stop particles whenever they are about to cross a predetermined time grid; these censused particles are then removed from the current simulation and stored into a time bank (the census bank).
Another useful triggering event is fission emission. In this fission census, neutrons emitted from fission reactions are removed from the current simulation and stored into a fission bank. One can see that this is actually a standard practice that has been long used in k-eigenvalue MC transport simulations. We can take a step further and census not only the fission neutrons but also the scattering neutrons-this results to collision census, which is typically used in the c-eigenvalue MC calculations [23].
There are several reasons to perform particle census. One is to limit particle population growth so that population control (discussed in more detail next) can be performed. Another reason is to allow the system (the MC model) to change-which can be geometry, composition, or parameter changes due to multi-physics feedback. Additionally, one can also see census as a manifestation of an iterative scheme performed to solve an equation-e.g., power iteration in k-eigenvalue problem.
It is worth noting that census time grid for population control does not necessarily need to be identical to other time grids possibly used in MC simulation. These other time grids include the one for tally scoring (also known as tally filters in some MC codes, such as OpenMC [24]), time grid for variance reduction techniques (e.g., weight window and forced precursor decay [2]), and census time grid for model change or multi-physics feedback. Given an initial population of size N , the objective of population control is to return a controlled final population with a size around, or exactly at, a predetermined value M , as illustrated in Fig. 1. In a supercritical system, typically N > M ; while in a subcritical one, N < M . The final population is then used as the source bank for the successive census-enabled transport simulation, during which a census bank is populated by a certain census mechanism (e.g., time census or fission census, as discussed in Sec. II.A). Once the transport is completed (i.e., both source and secondary particle banks are exhausted), the census bank becomes the initial population to be controlled by a PCT of choice. It is evident that population control does not care about what kind of transport simulation is being performed, whether it is a time-dependent fixedsource or an eigenvalue one. This also implies that any PCT can be used in any kind of transport simulation; as a particular example, one can use the particle combing technique [15] in k-eigenvalue simulation.

II.B. Population Control
The final population basically consists of copies of the initial particles, but how many times a particle gets copied will differ between particles, and some particles may not get copied at all. The procedure on determining how many times each initial particle get copied to the final population is the essence of PCT and has to be done in a way such that the MC simulation is not biased-i.e., the expectations of the population actions, and thus the expectations of simulation tally results, are preserved.
The only requirement for a PCT to be unbiased is to preserve the expected weight of each particle in the initial population. That is, for initial particle i having weight w i : where E[·] denotes the expectation of a random variable argument, d i is the number of copies of particle i in the final population, w i is the controlled weight assigned to the copies of particle i, and C i is the total weight represented by the copies of particle i in the final population. Now that we have described the minimum requirements-i.e., controlling population size from N to around M , while ensuring that Eq. (1) holds-we next point out two desirable characteristics of PCT.
The first is that we wish to have a low uncertainty in C i . In particular, we want its standard deviation, σ[C i ], to be as low as possible. In the absence of PCT, each initial particle i is copied once to the final population, and we will have C i = w i and σ[C i ] = 0; however, if a PCT is being used, σ[C i ] ≥ 0. The value of σ[C i ] affects the variance of the actions of particle i in the MC simulation, and later we numerically demonstrate that it ultimately affects the simulation tally results.
The second desirable characteristic is that we would like our PCT to preserve the initial population total weight W as much as possible; in other words, if W is the final population total weight: and we would like W to be close or equal to W . Booth [15] suggests that such strict equality of W = W is generally unimportant for neutron and photon transport, but it may be very important in charged particle transport. Therefore, we consider it a desirable characteristic, not a requirement, of PCT.
As a remark, PCT is a technique that takes an initial population of size N and total weight W and returns a controlled final population that: (1) has a size equal or close to M , (2) preserves the expected total weight of each particle (i.e., satisfies Eq.  (3) and (4) are desirable characteristics.

II.C.1. Combing (CO)
Perhaps the most standardized 1 PCT is the particle combing technique (CO). Per our classification, the "Simple Comb" proposed by Booth (Section II in [15]) is weightbased CO. CO techniques are best explained with graphical illustrations. Let us consider a population control problem with N =6 and M =4. Weight-based CO combs the population as shown in Fig. 2, where ξ is a random number (from 0 to 1) used to determine the offset of the initial tooth of the comb. Once the initial tooth location is set, the remaining teeth are spaced W/M apart. Per Fig. 2, Particles 1 and 5 are copied once, Particle 3 is copied twice, and Particles 2, 4, and 6 are not copied at all.
To ensure unbiased MC simulation (c.f. Eq. (1)), the copies of particle i are assigned with weight w i = W/M . Booth also proposes the "Importance-Weighted Comb" (Section III in [15]), which per our classification is importance-based CO. Importance-based CO is similar to the weight-based CO shown in Fig. 2, but instead of using w i for the particle axis, W/M for the distance between teeth, ξW/M for the offset of the comb, and final weight

II.C.2. Modified Combing (COX)
A modification of CO is recently proposed by Ajami et al. [10]. Different from the weight-based CO shown in Fig. 2, the weight-based COX combs the initial particle as shown in Fig. 4. In COX, instead of having uniformly-spaced teeth and sampling the offset of the whole comb, we allow the teeth to be non-uniformly spaced by offsetting each tooth with a different random number. The controlled weight w i assigned to the particle copies to ensure unbiased MC simulation (Eq. (1)) are identical to those of CO. Ajami et al. [10] provides limited discussion and demonstration on how COX compares to CO. In Sec. II.D.2, we discuss how COX may actually avoid a particular drawback of CO; and then in Sec. III, we discuss how that remedy comes at a significant expense.

II.C.3. Splitting-Roulette (SR)
Sweezy et al. [4] proposes the weight-based splitting-roulette (SR). In SR, we assign each initial particle i with splitting number s i . For uniform, weight-based, and importance-based SR, the values for s i are respectively M/N , w i /(W/M ), and u i /(U/M ). We split each particle i into s i + 1 copies, and then Russian-roulette the last copy with surviving probability s i − s i ; the function · denotes the floor function, which produces the greatest integer not greater than the variable argument.
Finally, to ensure unbiased MC simulation, the surviving particle copies are assigned with controlled weight w i , which happen to be identical to those of CO techniques.
SR techniques neither exactly produce a final population of size M nor exactly preserve the initial total weight W -however, they preserve the expectations. To exactly preserve the population's total weight W , Sweezy et al. suggest performing a weight normalization at the end of SR. This weight normalization can be applied to other PCTs that do not exactly preserve the population's total weight as well (e.g., uniform and importance-based CO). The significance of this PCT weight normalization is further discussed later in Sec. II.D.3.

II.C.4. Simple Sampling (SS)
Simple sampling (SS) is the typical PCT employed in k-eigenvalue MC simulations [21]. In SS, we simply sample M particles from the initial population to be the final population. For uniform SS, all particles have a uniform probability to be sampled at each draw; while for weight-based and importance-based SS, the probability for a particle to be sampled at each draw is proportional to its weight w i and the product of its weight and importance u i , respectively. Finally, to ensure unbiased MC simulation, the sampled particles are assigned with controlled weight w i which values happen to be identical to those of the other PCTs.

II.C.5. Duplicate-Discard (DD)
We identify the PCT proposed by Leppänen in [3] as the uniform duplicate-discard technique (DD), due to its mechanism of randomly duplicating or discarding particles to achieve the desired population size. In particle duplication (for N < M ), we first copy each initial particle once to the final population, and then, on top of that, randomly sample M − N particles from the initial population to be copied to the final population. In particle discard (for N > M ), we randomly sample N − M particles from the initial population; the sampled particles do not get copied to the final population, while the rest are copied once. Finally, the controlled weight w i that satisfies the unbiased MC simulation requirement is identical to that of the other uniform PCTs: One If the initial population has a uniform weight, the weight-based sampling is identical to the uniform sampling, since W = N w i . However, if the initial particles have varying weights, the weight-based sampling simultaneously functions as a variance reduction technique as well: particles having relatively large weights tend to be split into multiple copies, which leads to variance reduction; on the other hand, particles with relatively low weights tend to be Russian-rouletted, which may lead to more efficient computation by not spending time tracking small-weight particles. Nevertheless, particle weight does not necessarily indicate particle importance. If the initial particles are assigned with some importance values, the importance-based sampling offers more effective variance reduction than the weight-based.
One may argue that uniform sampling is the least optimal as it assigns all particles an identical splitting number or surviving probability regardless of their weights and importance. However, uniform sampling can be the most optimum choice in two cases.
The first is when the population has a uniform weight and unknown importance, which is the case in a fixed-source problem without any variance reduction technique and in the typical k-eigenvalue simulation where all the fission neutrons are emitted with uniform weight. The second case is when the MC simulation is already equipped with some variance reduction techniques, such as the weight window or the uniform fission site method [25], because particle distribution and weight profile of the population would be optimized already, such that particles can be treated equally by the PCT.
In other words, avoiding redundancy in variance reduction effort. In particular, in the application of an effective weight window or the uniform fission site method, the use of weight-based sampling may actually ruin the already optimized particle distribution.
The interplay between PCT and variance reduction technique briefly described above is outside the scope of this study. Furthermore, while the theoretical analysis performed in Sec. III is applicable to all sampling bases, only the uniform PCTs are implemented and tested in Secs. IV-VI.

II.D.2. Correlation Issue in CO
In Sec. II.C, it is interesting to observe that CO techniques only require one random number to perform the population control (as a comparison, SS and SR respectively require M and N random numbers); in other words, a single random number determines the fate of all particles in the population. This unfortunately yields correlation in the particle sampling. As an example, Particles 1 and 2 in Fig. 2 will never be sampled together. This correlation may produce unwanted behavior depending on how the initial particles are ordered. half. While such behavior may result in a larger tally variance, the expectation is still preserved. Outside this postulated PCT problem, some extent of physics is naturally embedded in the particle population order (e.g., adjacent particles may be originated from the same emission event). However, there has never been any observable effect of this correlation issue in the practical application of CO [2,17,20,8,19].
If one wishes to eliminate this possible correlation issue, the initial population order must be randomized before CO is applied. However, in massively parallel computation with reproducibility requirement, this pre-randomization process will require a large number of communications, which may ruin the parallel scalability of the simulation.
The modified combing technique COX proposed by Ajami et al. [10], to some extent, remedies this correlation issue as demonstrated in the lower part of Fig. 5. Nevertheless, this remedy comes at the expense of increasing σ[C i ], which is discussed later in Sec.

II.D.3. PCT Weight Normalization
Some PCTs-i.e, uniform and importance-based PCTs, and all SR techniques-do not exactly preserve the population total weight W . However, the expectation of the total weight is still preserved because where the first and the second equalities respectively use Eqs. (4) and (3).
To exactly preserve W , Sweezy et al. [4] suggest performing weight normalization after population control is performed. This is done by multiplying all of the final particles with the factor W/W , so that C PCT weight normalization introduces bias as Eq. (1) is now violated: where the inequality comes from Jensen's inequality [26,4] suggesting E[W/W ] ≥ 1.
Nevertheless, it can be seen that by using a large number of particles, the bias in the normalized PCTs can be minimized; however, it is also the case for the lack of exact total weight preservation in the non-normalized PCTs. In other words, PCT weight normalization suggested in [4] is only recommended if preserved total weight is more important than unbiased MC simulation.

II.D.4. More Advanced PCTs
The techniques considered in this work are those of basic PCTs. More advanced PCTs include the one proposed by Booth in Section IV of [15], which introduces the idea of partial population weight adjustment-an unbiased alternative to the weight normalization proposed by Sweezy et al. [4] (see Sec. II.D.3)-to exactly preserve the population's total weight W . This partial adjustment is technically more advanced than the weight normalization technique; it introduces tunable parameters (i.e., the adjusted partial population size and the number of recursive partial adjustments) and additional challenges for parallel computing implementation. While the proposed partial population weight adjustment is applied to the importance-based CO in [15], it basically can be applied to other PCTs that do not exactly preserve W as well. Other developments of advanced PCTs include the more recent study by Legrady et al. [19], which introduces several advanced CO techniques specifically improved for extensive variance reduction.

III.A. Theoretical Analysis
By determining the first and second moments of C i (the total weight of the copies of initial particle i in the final population), we can determine the variance introduced by a PCT: Another and perhaps more illustrative quantity is the relative uncertainty (standard deviation) introduced by the PCT to each particle i in the initial population: Unless normalized (as discussed in Sec. II.D.3), all of the identified PCTs (SS, SR, CO, COX, and DD) are unbiased, which means E[C i ] = w i . However, the second moments E[C 2 i ] of the PCTs may be different and thus become the key to determine the relative performance on how large uncertainty σ r [C i ] is introduced by the techniques.
In SR (described in Sec. II.C.3), each initial particle i is either copied s i + 1 times with a probability of s i − s i , or otherwise copied s i times. This suggests where we note that w i = w i /s i .
In CO (described in Sec. II.C.1), each initial particle i is either copied s i − 1 times with a probability of s i − s i , or otherwise copied s i times. The quantity s i (splitting number) used in this context happens to be identical to that of SR; and the function · denotes the ceiling function, which produces the smallest integer not smaller than the variable argument. Following the similar process to that of SR in the previous paragraph, we obtain In SS (described in Sec. II.C.4), each particle i can be copied multiple times up to M ; this means where we use the same definition of s i used in the other PCTs. Per binomial theorem, we can find that and thus where the approximation is due to the fact that typically s i M (or equivalently N 1 for uniform PCTs).
In uniform DD (described in Sec. II.C.5), we have two different cases. In the case of N > M , we uniformly discard N − M particles from the initial population. Therefore, particle i has to survive all of the discard draws to get copied once, otherwise it will not get copied at all. This means, for N > M we have where again we use the same definition of s i used in the other PCTs. On the other hand, in the case of N < M , DD keeps M/N copies of the initial population, and then uniformly draw a particle duplicate (M mod N ) times out of it. This process is similar to that of SS, except that we sample (M mod N ) particles instead of M particles and we pre-keep M/N copies of each initial particle. This gives where the approximation is again due to N 1.
In COX (described in Sec. II.C.2), things are more involved in that deriving the relative uncertainty σ r [C i ] is not as straightforward. First, let us observe how Fig. 4 of COX differs from Fig. 2 of CO. We can see that Particle 1 suffers from the same uncertainty in both methods; in this case, both CO and COX introduce identical uncertainty to Particle 1. However, it is not the case for the other particles. For example, in CO, Particle 2 has only two possibilities, get copied once or not at all; but in COX, Particle 2 has an additional possibility, which is to get copied twice. Due to this additional possibility, COX introduces higher uncertainty to Particle 2 than CO does. Similar findings can be observed for Particle 3. These observations indicate that CO , depending on how the particle i is located relative to the comb grid (the broken line in Fig. 4).   On the upper part of Fig. 6, we have s i ≤ 1 and 0 ≤ δ i ≤ 1 − s i ; in this case, COX and CO are identical. On the lower part of Fig. 6, we have s i ≤ 1 and 1 − s i < δ i ≤ 1; in this case, we have On the upper part of Fig. 7, we have s i ≥ 1 and 0 < δ i ≤ s i − s i ; in this case, we have Finally, on the lower part of Fig. 7, we have s i ≥ 1 and s i − s i < δ i ≤ 1; in this case, we have   16) and (17) for DD-are plotted in Fig. 9. Different to those of the other PCTs, σ r [C i ] of COX is dependent on δ i as shown in Fig. 8; thus, in Fig. 9, we plot its average value and shade the region (min to max) of its possible values. The x-axis is chosen to be 1/s i , which is equivalent to the ratio w i /w i -or N/M for the uniform PCTs. This x-axis effectively represents a measure of the system's population growth, which is dependent on the system criticality and the census frequency. Roughly speaking, one can say that N/M increases with the criticality of the system as illustrated with the arrows in the figure. The larger σ r [C i ], the larger the uncertainty introduced by the PCTs, which may lead to less precise (more statistical noise) results. From Fig. 9, it is evident that in a growing population regime ("Super"), the larger the ratio N/M , the larger the uncertainty introduced by the PCTs; this trend generally extends to the decaying population regime ("Sub"). However, some methods (SR, CO, and DD) take advantage of the pure-splitting scenario-in which M is a multiple of N -such that σ r [C i ] drops to zero. We note that given a reasonably accurate prediction of population decay rate, one can take advantage of this behavior to minimize uncertainty introduced by In the test problem, we perform population control to an initial population with a cosine statistical weight distribution: Each initial particle i is associated to tally bin i. All copies of particle i in the final population will score their controlled weight w i to the tally bin i; in other words, we are tallying C i .    Next, we would like to use the PCT test problem to verify the theoretical PCT uncertainties σ r [C i ] derived in this section. We set the target size M to be 1000 and consider multiple values of N such that N/M ranges from 0.75 to 1.25. In each case, the population control is repeated 100 times so that we can determine the relative standard deviation σ r [C i ] based on the accumulation of C i and C 2 i . Furthermore, we randomize the particle order in the population "stack" at each repetition. In uniform PCTs, σ r [C i ] is independent of i as it only depends on the value of N/M , as shown in Fig. 9. Therefore, in each case of N/M , we take the average of σ r [C i ] over all i as the final result. Finally, these numerical results from all cases of N/M are compared to the theoretical values, as shown in Fig. 12. The numerical results are denoted by the markers, and the lines are the theoretical values identical to those in Fig. 9. Excellent agreement is observed, even for COX with its ranging theoretical σ r [C i ] (the shaded area). This verifies not only the theoretical σ r [C i ] derived in Sec. III.A, but also the PCT implementations.

IV. Parallel PCT Algorithm
Romano and Forget [21] introduce an efficient, reproducible, parallel fission bank algorithm for k-eigenvalue MC simulation; in the paper, the typical uniform SS is described as the PCT. However, per our discussion in Secs. II.A and II.B, we can actually apply the algorithm not only to the k-eigenvalue MC simulation (fission census) but also to the time-dependent fixed-source with time census. This allows us to adapt the algorithm and design a common population control code routine for both simulation modes. Furthermore, the PCT of choice can be any of the five PCTs discussed in Sec.
II.C. In this section, we describe the adaptation of the parallel particle bank algorithm proposed in [21], present the resulting pseudo-codes of the different PCTs, and perform a weak scaling study to their implementations. Fig. 3 in [21], Fig. 13 illustrates an example of how particle banks are managed and population controlled using the proposed parallel algorithm. In the example, we consider 1000 source particles evenly distributed to 4 processors-each processor holds a Source Bank of size 250. The source particles are then transported in parallel. The transported particles are subject to a census mechanism, which can be a time census for time-dependent simulation or fission census for eigenvalue one.

Generalized from
Once the particle census is completed, population control is performed to the Census Bank using one of the PCTs (SS, SR, CO, COX, or DD). Finally, the resulting final population (Sampled Bank) is evenly redistributed to the processors via the nearestneighbor bank-passing algorithm, where each processor only needs to communicate (send or receive) with its adjacent neighbors as needed, without any global particle bank formation nor typical master-slave communication [21]. Two exclusive scans need to be performed in the proposed parallel algorithm. An exclusive scan to the Census Bank is required to determine the total size N and the position of the processor's local bank relative to the "global" bank, so that reproducible population control, regardless of the number of processors, can be achieved by consistently following the same random number sequence. The other scan is performed to the Sampled Bank so that we can determine local bank offsets required to perform the nearest-neighbor bank passing.
Algorithms 1 and 2 respectively show the pseudo-codes for bank-scanning and bankpassing processes, which are used in all of the PCT algorithms: Algs. 3-7. The PCT algorithms only take the minimum information required to perform the population control-the Census Bank (which can be either fission or time bank) and the target size M -and return the controlled, evenly distributed across processors, final bank. Therefore, the proposed parallel PCT algorithms are applicable for both time-dependent fixed-source and eigenvalue MC simulation modes. We also note that the algorithms are designed to start and return with the same random number seed across all processors, which is important for maintaining reproducibility. Next, we perform a weak scaling test to assess the relative parallel scalabilities of the different PCTs. The test is similar to the verification test problem in Sec. III.B, except that M is set to be 10 5 times the number of processors, N ∈ [0.5M, 1.5M ] is randomly picked in 50 repetitions, and the initial particles are randomly distributed to the processors. We note that there is no particle transport performed; we specifically measure the runtime required to perform population control and manage the particle bank, which includes communicating attributes and members of the bank throughout the processors. repetitions. It is found that SR, CO, and COX identically scale the best, followed by DD, and then SS. The sampling mechanisms of SR, CO, and COX are embarrassingly parallel, which make the techniques scale very well. On the other hand, SS, which scales the worst, needs to serially sample all of the M particles. The sampling mechanism of DD is also done in serial, but it only needs to sample as many as the difference between N and M ; this makes DD scale better than SS, and also explains why DD's runtime has a relatively larger standard deviation. Finally, it is worth mentioning that we purposely picked a modest number of particles per processor (10 5 ) to get a balanced demonstration on the significance of work (population control sampling) and communication. Should we significantly increase the number of particles per processor, such that the amount of work far outweighs the amount of communication, parallel scalabilities of SR, CO, and COX will improve, approaching the perfect scaling (horizontal line).

V. Time-Dependent Problems
In this section, we devise time-dependent MC test problems and then solve them with the PCTs to assess their relative performances. We adapt the homogeneous infinite 1D-slab medium problem of the analytical time-dependent benchmark suite AZURV1 [28]: which is subject to the following boundary and initial conditions: Note that particle position and time are respectively measured in mean-free-path where v is particle speed; and we also have the typical scattering parameter c and H(·) denotes the heaviside function.  The test problems are initiated by an isotropic neutron pulse at x = t = 0 (c.f. (22)).
In both cases, the scalar flux solution gradually diffuses throughout the medium. The difference is that the significant neutron absorption promotes population decay in the subcritical case. On the other hand, while the solution of the supercritical case initially behaves similarly to that of the subcritical, it eventually raises up due to the significant fission multiplication-at t = 20, the population size reaches exp(2) = 7.39 times of the initial value.

V.A. Verifying Time-Dependent Features of the MC code
To the authors' knowledge, there are three different time-dependent scalar flux quantities that can be calculated via MC simulation: (1) spatial-average time-average φ j,k , (2) spatial-average time-edge φ j (t), and (3) spatial-edge time-average φ k (x), where j and k respectively denote spatial and time mesh indices. The first tally uses the typical track-length estimator averaged over time mesh. The second uses a time-edge estimator, which accumulates the product of neutron speed and weight whenever a time-grid is crossed, averaged over spatial mesh. The third uses the typical spatialmesh-crossing estimator, which scores particle weight divided by absolute of normal product of particle direction and the surface, averaged over time mesh. The use of the track length estimator (spatial-average time-average φ j,k ) is typically desired because it generally gets more samples compared to "event-triggered" estimators like time-edge-crossing (φ j (t)) and spatial-mesh-crossing (φ k (x)). Nevertheless, it is worth mentioning that in some applications, time-edge quantities φ j (t) may be more desired than the time-average one φ j,k .
To simulate the supercritical (c = 1.1) and subcritical (c = 0.9) cases, we consider purely fission media with ν = c. The test problems are simulated using the research MC code, and we record the scalar flux using the three tally estimators that are Additionally, this also suggests that this set of test problems, the AZURV1 benchmark [28], serve as a good verification tool to test time-dependent features of MC codes.    Table I and Fig. 9). Table I shows census period and the expected ratio N/M associated with the simulated cases. By referring to Fig. 9, we can estimate the expected uncertainty σ r [C i ] introduced by a PCT at a given value of N/M . For convenience, plots showing the expected σ r [C i ] associated to the simulated cases are provided in Fig. 17. Note that this expected uncertainty is introduced every time the population control is performede.g., with 4 number of censuses, we perform census and population control, and introduce the associated uncertainty, once every 4 mean-free-times. This means, smaller σ r [C i ] due to larger census frequency does not necessarily lead to smaller uncertainty in the simulation result, because the more frequently we perform population control, the more frequently we introduce the uncertainty σ r [C i ] (even though small) to the population.
Two performance metrics are considered: (1) total runtime T and (2)  is calculated over the N r = 100 repetitions: where the superscript (i) denotes repetition index. We note that a repetition (or a realization) can be seen as a batch of source particles that makes a single, independent history. Finally, the averaged sample standard deviation of the scalar flux is calculated as follows:σ The resulting performance metrics are shown in Fig. 18. A figure of merit (FOM) based on the two performance metrics, is also shown in the figure. Finally, the analog (without PCT) solution, also run in 100 repetitions, is shown in the figure as well as a reference point. Figure 18      Finally, figure (c) of Fig. 18 shows that the FOMs of all PCTs are always lower than that of the analog simulation, and they monotonically decrease as we increase the census frequency-it seems that PCT is parasitic in this MC simulation. However, we should note that the main reason of applying PCT in a supercritical system is to limit population size being tracked in the simulation. Nevertheless, in general, one may find situations where the PCTs have larger FOMs than the analog one (for smaller census frequency) should the advantage of runtime reduction significantly outweighs the uncertainty introduced by the PCTs. Another important takeaway from the figure is that SR, CO, and DD are in the same ballpark as the best PCTs, which are followed by COX, and then SS.

V.B.2. Subcritical problem
The main motivation of population control in a subcritical problem is to maintain population size so that we have enough samples to yield more precise (less noisy) solution. However, this comes at the expense of increasing overall runtime as more neutrons need to be tracked. Fig. 18 shows that applying PCT in a subcritical problem increases overall runtime, and it increases further as we perform the population control more frequently. It is also worth to mention that DD has a similar runtime to those of SR, CO, and COX in higher number of censuses; this is because DD only needs to sample as many as |N − M |, which gets closer to zero as we increase census frequency. Finally, figure (c) of Fig. 18 shows that the PCTs offer improved FOMs relative to the analog. The FOM is improved further as we perform population control more frequently. However, it starts to consistently degrade as the effects of the increasing runtime and of the significant uncertainty introduced by the PCT start to dominate.

Figure (a) of
Note that this is similar to the typical trend of a variance reduction technique: it helps to improve FOM, but will degrade FOM if it is used too much. Another important takeaway from the figure (c) of Fig. 18 is that-similar to the supercritical case-SR, CO, and DD are in the same ballpark as the best PCT, followed by COX, and then SS.

VI. k-Eigenvalue Problem
In a MC calculation, the k-eigenvalue transport problem is typically solved via the method of successive generations. The MC simulation involves accumulation of fission neutrons in a fission bank. At the end of each generation (i.e., when the fission census is completed), the eigenvalue k is updated, and the generated fission bank is normalized such that its total weight is identical to the target population size M , which is the number of histories per generation. Finally, the normalized fission bank is set to be the source bank for the subsequent generation.
The procedure described above is considered to be the "analog" approach where there is no PCT being used. The eigenvalue update and weight normalization are useful for the MC simulation, because they direct the neutronics system into the steady-state configuration which helps maintaining the number of source particles It is worth mentioning that the "eigenfunction normalization" described in the previous paragraphs and the "PCT normalization" discussed in Sec. II.D.3 serve different purposes. The eigenfunction normalization is a necessary step to ensure that scores accumulated into simulation tallies are not arbitrary in magnitude. On the other hand, PCT normalization is an optional step to preserve the total weight of the initial population passed to the PCT (at the expense of introducing bias, as discussed in Sec. II.D.3). As another clear distinction, the eigenfunction normalization is performed before we apply PCT, while the optional PCT normalization is performed after.
Similar to the PCT application in time-dependent MC simulation, PCT application in k-eigenvalue simulation would also introduce the uncertainty σ r [C i ] to the population. However, it should be emphasized that this uncertainty introduced by PCT is not a bias. There is a well-known bias associated with the method of successive generations (which can be mitigated given sufficiently large number of histories per generation [29]). This bias is introduced when the fission bank is normalized, which is in effect regardless whether population control is applied. PCT neither enhances nor eliminates this bias-it merely introduces the additional uncertainty σ r [C i ] to the already-bias simulation.
In an eigenvalue simulation, how many times population control is performed is determined by the total number of generations, which is typically a very large number (in the order of 10 2 to 10 3 ). This means, the uncertainty σ r [C i ] would be introduced by the PCT to the population so many times, which, according to the findings in the previous section, may lead to highly noisy solutions. However, the effect of the uncertainty introduced by the PCT on an eigenvalue simulation is expected to be much less pronounced than that in a time-dependent one. This is because once the eigenvalue convergence is achieved, we essentially simulate a steady-state system, where the ratio N/M is expected to be around and close to unity, in which most PCTs introduce minimum uncertainties, as shown in Fig. 9. Nevertheless, some PCTs (SS and COX) still introduce considerably high uncertainties even with N/M ≈ 1. In this section we would like to investigate the relative performances of the identified PCTs in a keigenvalue transport calculation-particularly, we would like to see whether there is any discernible effect due to the different magnitudes of uncertainty σ r [C i ] introduced by the techniques.
We consider the k-eigenvalue transport problem of the mono-energetic two-region slab medium from [30]: Similar to Sec. V, all physical quantities will be presented in the unit of mean-freepath. The first and the second regions respectively occupy x ∈ [0, 1.5] and x ∈ [1.5, 2.5].
The cross-sections of the two regions are νΣ f,1 = 0.6, Σ s,1 = 0.9, νΣ f,2 = 0.3, and Σ s,2 = 0.2. Finally, the two-region slab is subject to vacuum boundaries. By using a deterministic transport method, Kornreich and Parsons [30] provide reference values for the fundamental k-eigenvalue, k = 1.28657, and the associated scalar fluxes at certain points (shown in Fig. 19).   follows the definition in Eq. (29). It is found that SS-the simplest, yet the most well-known technique for k-eigenvalue simulation-is the least performant PCT having FOM about 24% lower than the analog technique. On the other hand, SR, CO, and DD are the best PCTs with FOM about 2 to 4% lower. However, it is worth mentioning that the discernible decrease in FOMs of SR, CO, and DD is due to their higher runtimes (right figure of figure (b)), since their sample standard deviations ( figure   (a)) are about the same as the one with the analog technique. The higher runtimes of the PCTs can be hidden should we run the simulations with a sufficiently large number of particles per processor. In such a situation, SR, CO, and DD would be as performant as the analog technique.

VII. Conclusion
An extensive study on population control technique (PCT) for time-dependent and eigenvalue Monte Carlo (MC) neutron transport calculations is presented. We define PCT as a technique that takes a censused population and returns controlled, unbiased one. A new perspective based on an abstraction of particle census and population control is explored, paving the way to improved understanding and application of the concepts. We discuss how different kinds of census-e.g., time, fission, and collision censuses-are performed in time-dependent and eigenvalue problems. We also discuss the requirements and desirable characteristics of a PCT.
Identified from the literature, five distinct PCTs are reviewed: Simple Sampling (SS), Splitting-Roulette (SR) [4], Combing (CO) [15], modified Combing (COX) [10], and Duplicate-Discard (DD) [3]. While SS has been the typical procedure of choice in handling fission bank in MC eigenvalue simulations, the other four techniques have been almost exclusively applied for time-dependent simulation. The review encompasses the basic procedures of the techniques, significance of their sampling bases (uniform, weight-based, and importance-based), bias in PCT weight normalization, possible correlation issue in CO, thorough characterization of the recently introduced COX, and relation to the more advanced PCTs [15,19]. A short remark of the five PCTs, highlighting their respective caveats, are summarized in Table II. A theoretical analysis on the uncertainty introduced to population by each of the  (2) preserve total weight CO Subject to possible undesirable behavior due to correlation in particle order COX Avoid CO's issue, at the expense of increased uncertainty introduced PCTs is presented. The resulting theoretical uncertainties (shown in Fig. 9) are useful not only for theoretically assessing the relative performance of the PCTs, but also for numerically verifying whether the techniques are appropriately implemented. It is found that CO and SR are equally the most performant techniques based on this metric (smallest uncertainty introduced), followed by DD, then COX, while SS introduces the largest uncertainty. We hypothesized that this uncertainty would proportionally affect simulation tally results, which was later confirmed when we ran test problems using the different PCTs.
Parallel algorithms for the five PCTs are proposed. The algorithms are based on a generalized version of the parallel fission bank algorithm [21] and designed to be applicable for both eigenvalue and time-dependent simulations. The use of abstract base class for streamlined implementations of the five PCTs is also suggested. Weak scaling results of the PCTs are performed to demonstrate the parallel scalability of the techniques. It is found that SS and DD have limited scalabilities due to their respective significant serial sampling procedures.
Supercritical and subcritical time-dependent test problems based on the analytical benchmark AZURV1 [28] are devised; we found that these test problems serve as a good benchmark suite for verifying time-dependent features of a MC code. With the test problems we not only compare the relative performances of the PCTs, but also demonstrate typical behaviors of the PCTs in supercritical and subcritical problems as a function of census frequency. Two performance metrics are considered, total runtime and averaged sample standard deviation of the scalar flux, which are aggregated into figure of merit (FOM). Similar analysis is performed to a k-eigenvalue test problem.
We found that SR and CO are equally the most performant PCTs, closely followed by DD, and then COX.
The results of the time-dependent and k-eigenvalue test problems demonstrate the superiority of SR and CO. However, that does not mean that everyone would be confident in using one of those techniques in all of their MC simulations. This is particularly true since there is caveat for each of the PCTs, as summarized in II. The proposed generalized and streamlined PCT parallel algorithm offers MC code developers a minimally-invasive way to implement all of the PCTs into their code and allow the code users to pick the technique themselves. Finally, it is worth mentioning that the proposed PCT parallel algorithm can be generally applied to other MC eigenvalue simulations, such as c- [23] and α-eigenvalues [31].
Future work includes implementing the generalized parallel PCT algorithm into a production MC code and assessing the relative performances of the PCTs in simulating more practical, multi-dimensional, continuous-energy problems. Additionally, investigating the interplay between PCT and variance reduction technique (VRT), which is briefly discussed in Sec. II.D.1, can be a potential future work. Finally, a recent analysis of neutron clustering in k-eigenvalue transport calculation demonstrates that sampling source sites without replacement (which is a variation of DD) versus with replacement (which is SS) can greatly reduce the degree of clustering [32]. Studying how the different PCTs (particularly those who are not included in [32], i.e., SR, CO, and COX) affect the clustering phenomenon can be a potential future work as well.