Geant4 simulation of the moderating neutrons spectrum

In this paper we deal with the problem of predicting a steady-state neutron spectrum in media of arbitrary composition and geometry. The analytical calculations of such spectrum are often too complex, if at all possible. We describe a method of Geant4-based Monte Carlo calculation of the steady-state neutron spectrum in a medium containing a fixed neutron source. In addition to the steady-state spectrum, we obtain the snapshots of the neutron spectrum evolution in time, which may be thought of as the non-equilibrium neutron spectra, and their form is of considerable interest for further studies.


Introduction
The detailed knowledge of the neutron spectrum is crucial for numerous applications such as the nuclear reactor operation [1][2][3], the traveling wave reactor (TWR) development [4][5][6], including the search of the neutron energy ranges suitable for the wave nuclear burning [7,8], the search and prediction of the so-called "blowup modes" in neutron-multiplying media [9,10], the verification of neutron moderation theories and so on.
However, obtaining the neutron spectrum for a particular combination of geometry, medium composition, set of nuclear reactions etc. is a non-trivial and rather complex task by itself. Therefore, in practice, various assumptions and approximations are usually made. For example, in nuclear reactor physics the neutron spectrum is often composed of separate parts -the fission spectrum in high-energy range (fast neutrons), the Maxwellian spectrum in low-energy range (thermal neutrons), and the Fermi spectrum in between. These three parts are then joined together using some scaling coefficients [3]. For the same reasons, e.g. in [7], a well-known fission spectrum was used for the fast TWR calculations. In [8] the experimental WWER 1 spectrum was used for the epithermal TWR calculations. So there is currently no straightforward way to obtain a complete form of the neutron spectrum within a single theory. It is not surprising, because such theory is not at all easy to develop, and the well known and most cited theories of neutron moderation were largely developed with nuclear reactors in mind ("...raison d'être of slowing-down calculations (or experiments) is to provide information for criticality calculations." [1]). Nowadays, there are much more sophisticated and advanced theories, each working more or less satisfactorily in its area, but the practical considerations most often prevail [11,12].
In the absence of a convenient theory of neutron moderation, another possibility to deal with complex geometries, different medium compositions, and temperatures, is to use some Monte Carlo simulations which allow to tune all these parameters and see how they influence on the shape of the resulting steady-state neutron spectrum.
In Monte Carlo simulations of the neutron moderation by some medium the following approach is often used. The initial neutron spectrum is set, and the neutrons are directed to some moderating layers (optionally of variable thickness). The passage of neutrons through this medium is tracked, and the spectrum of neutrons coming at the opposite side is measured (Fig. 1). By comparing the resultant spectrum to the initial one for different moderator thickness, some conclusions can be made on how the moderator layers influence the initial neutron spectrum [13][14][15][16][17].
However, such approach is rather limited, and is disadvantageous for our task. The problem is that the Monte Carlo codes such as Geant4 treat the injection of each primary particle as a separate Event, and simulate its further interactions independently of all the other primary particles. It means that each primary particle is injected into the system at "zero" time (i.e. the particle enters a "fresh" system with no signs of other particles' presence or previous interactions), and is then tracked down to its disappearance either by absorption or by simply escaping the volume of interest. This is of course accompanied by the tracking of all the secondary particles if there are some.
Repeating this scenario N times for a large number of primary particles is actually equivalent to injecting a single pulse of N primary particles into the system. And the spectrum thus obtained is therefore a spectrum of a single "generation" of neutrons which were born (injected) at the same time, then passed through some moderating layers, and eventually formed some final distribution at the output.
But how would the form of the final spectrum change if we used a fixed neutron source, and the new particles were born and added to the system all the time? The new ("younger") particles would obviously possess their initial energy, and they would be added to the "older" particles which had already slowed down to lower energies. The total spectrum of such system will strongly depend on the time that the neutrons possess a certain energy (or the time the neutrons spend in a certain energy bin), as this determines the appearance of the spectral maxima and minima. For example, the maximum is expected to appear in the energy range in which the neutrons spend most of their time, and the fixed neutron source will make them gradually accumulate there (older ones + younger ones + even younger ones etc.).
For this reason, in problems involving the action of a fixed neutron source, it is necessary to work out some other way of calculating the neutron spectrum in Geant4 with due account taken of time.
In this paper we focus on the problem of predicting the steady-state neutron spectrum which sets in some medium under the action of a fixed neutron source, using the philosophy and capabilities of the Geant4 library [18, Chapter V. Tracking and physics], and restricting the consideration to the subcritical cases only (simulation of the super-critical conditions in Geant4 requires some special treatment which is not the subject of the present paper). In particular, in Sections 2-3 we describe our approach to the steady-state neutron spectrum calculation in Geant4. In Section 4 we perform the basic test of our algorithm for the simplest case of neutron moderation in pure hydrogen. Next, in Section 5 we discuss the advantages of the use of the time slices with the logarithmic temporal step, and describe a method of their implementation in Geant4. In Section 6 we also describe an alternative way of building the steady-state neutron spectrum without the time slices, and discuss its pros and cons. Finally, in Section 7 we discuss some technical details of the possible simulation optimizations.

The concept of time slices of neutron spectrum
In the framework of Monte-Carlo method it is possible to track the change in energy of a single neutron with time in certain medium. This procedure may be repeated many times from the start, imitating a simultaneous injection of a large number of neutrons into the studied system. This also allows to set the initial neutron spectrum of any form.
Let us conventionally call such a single portion of neutrons simultaneously injected into the system, a single "generation of neutrons".
If we know the neutron energies at any time, it is possible to build the "snapshots" of their spectra at different moments in time. We shall call them the time slices of the neutron spectrum. Such time slices will show the temporal evolution of the neutron spectrum in some medium with chosen parameters (Fig. 2).
Within the framework of Monte-Carlo method, it is reasonable to assume that the following neutron generations will demonstrate a similar spectrum evolution, provided that all of them are produced under similar conditions. This means that starting with zero time, after a period of ∆t, the spectrum of the first generation of neutrons will be n(E, ∆t). During this period a second generation of neutrons will be injected into the system from the neutron source with initial spectrum n(E, t = 0). After another period ∆t the first generation of neutrons will be moderated to n(E, 2∆t), and the second one -to n(E, ∆t), and also the third generation of neutrons will be injected with their initial spectrum. And so on and so forth.
After some time T max , when the first generation of neutrons is totally thermalized and absorbed (or escapes from the system), a stationary state will settle -each leaving generation of neutrons will be replaced by a younger one. The total neutron spectrum will stop changing, and may be considered a non-equilibrium stationary spectrum for the considered combination of the chosen medium and neutron source. This is of course true for the case of a constant neutron source only.
To determine such non-equilibrium stationary neutron spectrum is the main goal of the present research.
So in the volume which contains all of the mentioned neutron generations at once, the total spectrum is composed of the spectra of all these generations: where T max is the neutron lifetime in the studied system, n(E, t) is the spectrum of neutrons which spent the time t in the system. Or in the discrete case: where n i (E) is the spectrum of neutrons in i th time slice, ∆t is the time bin width. The time bins are not necessarily of equal width, so in general case, where (∆t) i is the width of the i th time bin. So, in order to apply all these considerations to some Monte Carlo simulation, it is necessary to find a way to calculate the mentioned spectrum time slices at chosen moments in time.

Time slices in Geant4
The Geant4 Monte-Carlo engine [19][20][21] generates the primary particles only one at a time. The particle's trajectory is being tracked from the moment of its appearance in the system -till the moment of its vanishing or complete stop [18]. After that another primary particle is generated and its trajectory is tracked. The trajectory thus obtained is a polygonal line consisting of consecutive steps. Each step is limited by its PreStepPoint and PostStepPoint.
Each StepPoint stores the information on: • particle position • particle momentum direction • particle kinetic energy • global time • and some other useful quantities Since each step is generated randomly, it is impossible to predict (or preset) the moments t at which the particle appears in PreStepPoint or PostStepPoint (these points are chosen randomly according to the physical processes applied). Therefore it is also impossible to obtain the neutron spectrum "snapshot" at arbitrary moment in time, as was described in Section 2. However, due to the fact that Geant4 stores the GlobalTime for each StepPoint of the particle track, it is possible to overcome this limitation.
One can divide the time axis into intervals -time bins -and using the values of energy and GlobalTime at each StepPoint, put the particle into the corresponding time bin (Fig. 3).
For the algorithm sketched in Fig. 3 a simple 2D array may be used. The dimensions of this array are the number of bins in the energy spectrum -numEnergyBins (often referred to as "channels" in spectrometry), and the number of time intervals we divided the time scale into -numTimeBins.
The GlobalTime of the current StepPoint as well as the current kinetic energy of the tracked particle may be accessed in the UserSteppingAction function of the G4UserSteppingAction class. The corresponding time and energy bins may be calculated as follows: 1 G4double currentTime = step->GetPreStepPoint()->GetGlobalTime(); 2 G4double currentKinE = step->GetPreStepPoint()->GetKineticEnergy();   1. The neutron is registered for the first time with zero GlobalTime. In this case it is considered to be the initial neutron, and its energy lets us to fill in the corresponding cell in the 2D time-spectrum array. 2. The neutron is registered once again with different energy and at different time. In this case we assume its energy did not change during all the intermediate time bins, and we fill in those array cells as well (Fig. 4a).
Here isFirstStep is true for the first fire of a neutron, lastTimeBin is the previous time bin, in which the neutron had been registered earlier, timeBin is the current time bin, in which the neutron is registered for the second time, lastEnergyBin is the previous energy bin, in which the neutron had been registered, the numTimeBins is the total number of time bins, and the pTimeSpectrum is the pointer to the 2D array containing the spectra. 3. The neutron hits the same time bin for the second time, but with different energy. In this case we have to subtract it from the previous energy bin and add it to a new energy bin (Fig. 4b).
4. The neutron is registered for the first time with non-zero GlobalTime. Such neutron is considered a secondary particle, so there is no need to fill in all the earlier time bins (Fig. 4c). This was actually the reason to check the isFirstStep value in the 2 nd point above.

Hydrogen test
As the analytical neutron moderation theory is best established for the hydrogen nuclei (namely, the protium), it seems reasonable to start the verification of the algorithm from this material.
Any simulation in Geant4 starts with the construction of the system's geometry and definition of all the substances and materials used. Here we use the hydrogen density at STP: By this test, the neutron moderation theory [22] taking into account the medium temperature may also be supplemented with the Monte Carlo simulations. For this purpose we are currently interested in simulating the neutron moderation in an infinite homogeneous medium first.
To simulate an infinite medium, we use a "Detecor" in a form of a cube with the side comparable to the size of the visible Universe (∼ 10 26 m).
1 G4Box *world_box = new G4Box("world", 1e26*m, 1e26*m, 1e26*m); All this cube will act as a moderator and a sensitive volume at the same time. It means that the entire cube is filled uniformly with hydrogen atoms, and the neutrons are being tracked throughout its entire volume.
Next it is necessary to define the physical processes. We choose the following set of processes described in Geant4 library: • electromagnetic interactions [21], • elastic hadron interactions [23], • QGSP BERT HP ("Quark-gluon String, Precompound, Bertini, High Precision") which uses the quark-gluon model for the hadron energies above 12 GeV [24,25] (such energies are not required for our current purpose, but makes the program more flexible and adjustable to future applications), and the high-precision data on neutron scattering, capture and fission from the G4NDL database [26]. RegisterPhysics(new G4HadronPhysicsQGSP_BERT_HP(0)); Finally, a PrimaryGeneratorAction class must be defined. This class describes the source of initial particles. For our test it may be a simple point source located at the center of the "Detector" cube:   Note that under such conditions (an extremely large moderator volume and the neutron source located at its center), there is no realistic way for a neutron to escape from the moderator. Therefore, the only possible way for a neutron to disappear from the system (and stop being tracked by Geant4 engine) is to be absorbed by the moderator nuclei.
An example spectrum of neutrons obtained this way is shown in Fig. 5 in three different forms. As can be seen from this figure, the most convenient visual form of the spectra is the log-log scale (Fig. 5c), which requires the use of the logarithmic energy bins at the stage of calculations. Such scale will be used in most cases below.
The corresponding energy bin should be calculated in the following way: where numEnergyBins is the number of energy bins, currentKinE is the current kinetic energy, minOrder is the minimum required magnitude order for the energy, maxEnergy is the maximum possible energy. This way we obtain a form of the spectrum as in Fig. 5c. In the context of nuclear reactor applications, a Watt fission spectrum may be used for the initial neutrons: where a = 1.036, b = 2.29 and c = 0.4527 for 235 U fission spectrum. Now everything is ready to run the simulation. We generate a sufficiently large number of primary neutrons whose initial energies are distributed according to a fission spectrum of 235 U , given above, and then apply the algorithm described in Section 3, to each of them. Finally, we get the time slices of the neutron spectrum of a single generation of neutrons (remember, all the primary neutrons were generated under exactly the same conditions and tracked independently, which is equivalent to injecting them into the moderator at zero time, as was discussed above). These time slices are shown in Fig. 6. Figure 7 shows the respective stationary neutron spectrum calculated from all the intermediate ones (the time slices) using Eq. (2). As can be seen in this figure, the stationary neutron spectrum in hydrogen moderator has exactly the expected form. Indeed, the comparison with the typical neutron spectrum in thermal reactor (inset in Fig. 7) shows that there is an extremely close correspondence, aside, of course, from the fact that our spectrum was simulated in the pure hydrogen medium, which explains the absence of the resonance area and a substantially lower level of fission spectrum.
This basic test of the neutron moderation in pure hydrogen, although may seem oversimplified, is still very important, as it validates our approach and proves that the algorithm described in Section 3 is able to produce correct results. Since this algorithm is rather generic, and does not impose any restrictions on the moderator's geometry and material, as well as on the form of the primary particle source and initial energy, it can be easily applied to arbitrary medium of any form and composition, and to any source of primary particles.
For example, in order to change the geometry of the studied medium, instead of using a simple cube (G4Box class) like we did in our test above, one can create instances of any classes describing solids -either the basic shapes (e.g. G4Tubs, G4Cons, G4Sphere, and so on), or even more complex tesselated solids. Geant4 allows also to put one solid into another (hierarchical structure), or produce the unions, intersections or subtractions of solid volumes [18, Section 4.1.2 Solids].
In our test we used a standard hydrogen from the NIST materials database (code G4 H). However, each of the created volumes may have its own material (described by the G4Material class) including its state, temperature, density, pressure, and of course, the isotope composition, which is very important in our case. In fact, the simulated medium may even include some fissile materials (acting as neutron multiplicators), as long as the sub-criticality condition is met. In this case, the multiplication of neutrons due to fission reactions will be  [27]. Note that the main image shows the neutron population (not flux) distribution over the logarithmic bins (bin width is proportional to energy), so the 1/E slope is consistent between the main image and the inset. treated correctly by the Geant4 engine, since the number of neutrons always remains finite (due to sub-criticality condition). A very handy G4NistManager class gives quick access to a rich collection of predefined widely used materials.
In order to change the parameters of the primary particle source (the particle type, its initial energy and momentum, and the position of its first appearance in the simulated volume, one can either set them for the corresponding G4ParticleGun instance while constructing a necessary G4VUserPrimaryGeneratorAction-derived class, or even use a universal G4GeneralParticleSource class to allow the user to load those parameters from an external macro file at run-time without the need of re-compilation.
Such flexibility of Geant4 allows us to apply our algorithm described in Section 3, to a variety of tasks beyond reactor physics, ranging from e.g. the radiation shielding and protection to medical imaging applications and radiation therapy.
Turning back to the analysis of the results of our hydrogen test, it is interesting to note the evolution of the spectrum shape with time, shown in Fig. 6 (please note that the spectra shown in this figure represent the "snapshots" of the spectrum of a single generation of neutrons, or alternatively, a single pulse of neutrons, taken at different times). For a different visual representation of the spectrum evolution it is possible to build a kind of pseudo-3D (or even a 3D) plot (Fig. 8). As can be seen in Fig. 8, the distribution of neutron energies shifts to the thermal area very quickly, and then its position remains the same for quite a long time, and it just gradually fades out as the neutrons get absorbed by the medium. This picture served as a motivation to the use of a logarithmic time scale described in the next section. It is hard to see any details of the spectrum evolution in such view. Therefore, by analogy to what we have done with the energy distribution earlier, we can apply the logarithmic scale for the time axis here. The basic idea is to make the earlier time bins short, and the later ones -long. As long as the spectrum changes noticeably at first stages only, and the later time slices do not differ drastically from each other, this procedure will not introduce any serious error.

Logarithmic time scale
The use of the logarithmic time scale has also some technical advantages. In order to catch the details of rapidly changing neutron spectrum at early stages, the time bin width ∆t must be chosen rather small. On the other hand, in order to track the neutrons till the very end of their lifetime in a weakly-absorbing medium, the maximum simulated time T max must be chosen rather large. Putting these two conditions together leads to a large overall number of time bins required. It can make the 2D TimeSpectrum array too large and inconvenient to operate.
In case of the energy scale it was not necessary to start with zero, so it was possible to choose some minimum energy and simply calculate the logarithm from there on. As to the time scale, it must start from zero, so applying a logarithm requires some special care. Here we suggest to choose the maximum time T max and the desired number of time bins m, and then by imposing a condition that the first and the second bins must be of equal width, find the start of the logarithmic time scale T min (see Fig. 9) as follows: Given the current time t, the corresponding time bin may be calculated as follows: bin = (m − 1) · log 10 (t/T min ) log 10 (T max /T min ) + 1 For the times less than T min the logarithm is negative, so all such values should go into the 1 st time bin (number zero in terms of C/C++ syntax).
In fact, 2 in Eq. (5) represents the ratio of the two adjacent bins, so in general case, it may be substituted by some other desirable value depending on the task.
Since the time bins are different in length here, the total stationary spectrum should be calculated using Eq. (3).
The time slices of neutron spectrum in hydrogen taken at logarithmic time steps (∆t) i are shown in Fig. 10 as 2D time slices, and in Fig. 11 as pseudo-3D graphs. These images, and especially Fig. 10, demonstrate that, in addition to the calculation of a steady-state spectrum, the suggested approach gives us access to the spectrum shapes at different moments in time (or rather, at different stages of moderation) in a wide range, starting from the particle entrance to the moderator medium, and on. It could be an interesting task to analyze the statistical parameters of those distributions to study the process of neutron moderation in more detail. The authors have made some investigations of the spectrum variance and entropy for several materials in [28].
The total stationary neutron spectrum calculated on the basis of the time slices thus taken is shown in Fig. 12. It is apparently somewhat more accurate than the one in Fig. 7 because of a) the more detailed account for the spectrum evolution at early stages, and b) the larger time scales available (10 ms versus 100 µs in Fig. 7).
One remark should be made here. For certain media, a production of secondary neutrons is possible (the case of neutron multiplying media mentioned earlier). Such neutrons appear for the first time in some later time bins (corresponding to non-zero time). Since the later time bins are wider than the earlier ones, the contribution of such neutrons may be overestimated when multiplied by the corresponding bin width (∆t) i . One option to overcome this drawback   may be to use a 2D array of real numbers instead of integers, and instead of incrementing its values by a unit per each neutron, increment them by a fraction of time bin during which the neutron actually existed.

No-timeslice approach
Finally, a couple of words may be said on the case when there is no need to study the intermediate stages of the neutron spectrum (its time slices), but there is still interest in the form of a stationary spectrum. In this case the algorithm described in Section 3 may be further simplified as follows.
Since each of the PreStepPoint and PostStepPoint provide the information on particle's GlobalTime t and kinetic energy E, it is easy to calculate the particle's contribution to the corresponding energy bin in the total stationary spectrum: This method does not require the introduction of the concept of time slices described above, and is therefore free of its accompanying complications. While it is still usable for calculating the final stationary spectrum, there are many cases when it will not suffice, because the stages of spectrum evolution may be of considerable interest.
An example of stationary neturon spectrum calculated this way is shown in Fig. 13 compared to the one obtained with time slices approach (Fig. 12). Figure 13: Stationary neutron spectra calculated through the time slices (red dots) and widthout the time slices (violet line). As can be seen, with sufficiently high number of primary particles, these curves coincide very closely. The random fluctuations in the low energy range here are more prominent than in Figs. 7 and 12 because of a much smaller number of primaries simulated, while the bin width in the lower energy range is orders of magnitude smaller than that in the higher energy range (logarithmic bin width described in Section 3, which is why the adjacent bins are filled very unevenly).
As can be seen in this figure, the steady state spectra calculated using both methods (with time slices and without them) are practically the same even for a rather small number of primary events generated (orders of magnitude smaller than those used to produce Figs. 7 and 12). This simple method without time slices, described in this section, is also roughly equivalent to the use of the built-in CellFlux scorer in Geant4. Indeed, to retrieve a neutron spectrum, one can create a number of identical CellFlux scorers and attach corresponding particle and kinetic energy filters to each of them (a scorer per each energy bin) [29]. Then, after the simulation, one can collect the values accumulated by all scorers to obtain a complete spectrum. Aside from the minor inconvenience of creating and operating the numerous scorers, the CellFlux scorer does not yet support attaching to real-world logical volumes (as of Geant4 version 11.1), so in order to use it, one has also to create a dedicated scoring mesh to fit the geometry of their detector. Recent versions of Geant4 (starting with 10.7) feature some classes for easier accumulation of histograms. In particular, with the new G4TScoreHistFiller class, a primitive scorer can directly fill a 1-D histogram. This spares a user the necessity to deal with multiple scorers and collect all the values at the end of simulation, but also comes with its own limitations. For example, one can only use this histogram filler with the "real-world" or "probe" scorers. Thus, the only possible way of utilizing the CellFlux scorer together with ScoreHistFiller at the moment is to attach them to a "probe" (which is essentially a cube). This was exactly the way we obtained the data shown with light orange shaded area in Fig. 13 for comparison. It should be noted that we had to limit the maximum length of a single step to 1 mm for neutrons in order to achieve the accurate results for the hydrogen at normal conditions.
The comparison shows that the form of the steady-state spectrum is exactly the same, but the downside of the approach with built-in CellFlux scorer is the impossibility to attach it directly to a logical volume of arbitrary form, which would obviously complicate things in more realistic use cases. Additionally, both the "timeslices" and "no timeslices" curves in Fig. 13 were calculated much faster than the one obtained with the built-in CellFlux scorer and a histogram filler on the same test machine.
And again, all the methods described in this section give absolutely no information about the intermediate shapes of the neutron spectrum, as opposed to our first, "time-slices approach" described in Section 3 above. Meanwhile, these intermediate "time slices" or "snapshots" of the neutron spectrum carry a lot of useful information for future studies (see e.g. [28]).

Parallelization 7.1.1. Shared memory parallelization
In order to utilize the shared memory parallelism in Geant4, one should first instantiate the G4MTRunManager instead of G4RunManager. This run manager creates several instances of the G4Run class and divides the task between them. The SetNumberOfThreads() method may be used to set the desired number of threads. Each instance of G4Run will collect the data independently, so it is convenient to make the pTimeSpectrum array a member of its subclass. At the end of the job the data collected by all instances of G4Run must be merged together. For this purpose the G4MTRunManager calls the master run's Merge() method with each of them. The spectrum arrays may be simply added like this:  After that the RunAction::EndOfRunAction() method can output the resulting arrays for the master run, which may be checked with the isMaster variable.

Distributed memory parallelization
The Geant4 simulation can also be parallelized over the distributed memory system, e.g. several computers with network connection. A G4MPImanager from Geant4's extended examples can be used for this 2 . First, it is necessary to open the MPI session just before instantiating the G4RunManager or G4MTRunManager: 1 G4MPImanager *g4MPI = new G4MPImanager(argc,argv); 2 G4MPIsession *session = g4MPI->GetMPIsession(); 3 4 //G4RunManager *runManager = new G4RunManager; 5 G4MTRunManager *runManager = new G4MTRunManager; In this case the MPI workload manager or the directly invoked mpirun utility will launch a separate Geant4 application on each machine with its own run manager.
At the end of the job the collected data must be sent over the network to the MPI process with rank 0. For this purpose, the mentioned Geant4 example provides the G4VUserMPIrunMerger class, which has the Pack() and UnPack() methods to send the collected data over network as a single 1D array (hence the names). The inherited RunMerger class must be instantiated and activated in the RunAction::EndOfRunAction() method: and stored in some ExternalXSDataset class inheriting from the Geant4's G4VCrossSectionDataSet class and implementing the IsIsoApplicable() and GetIsoCrossSection() methods.
These datasets must be instantiated in the PhysicsList::ConstructProcess() and attached to each process of interest: It is then the user's responsibility to supply the correct cross-sections for the specified temperature.

Conclusions
We considered the case of neutron moderation under the action of a fixed neutron source. The knowledge of the non-equilibrium stationary neutron spectrum is essential to a number of neutron-related physical problems.
We developed an algorithm of the stationary spectrum calculation based on the concept of spectrum "time slices", which provides the view of the spectrum evolution stages for a single "generation" of neutrons, as well as for the whole set of neutrons including all new generations coming from the constantly acting neutron source. We also showed a way to implement this algorithm in the custom software built on top of the Geant4 Monte Carlo simulation toolkit.
The developed software algorithm was tested against the known neutron moderation spectrum in hydrogen and demonstrated a good agreement. We also illustrated the advantages of the use of logarithmic time scale in such problems and suggested a simple and quick way of its implementation.
Since the described algorithm and its implementation are rather universal, this software may be used for various particular tasks with only minor modifications (as small as change the system geometry and the neutron source spectrum). Despite the fact that we tested the software for the infinite homogeneous medium only, it may as well be applied for any kind of heterogeneous structures.