A Classical and Semi-Analytical Analysis of Isotropic Scattering of Neutrons through a Protective Wall of Iron-Via Two Energy Groups

Introduction

In service to the field of health physics, the essential background topics in medical physics, and the reviews of radiological safety of reactors of nuclear reactor operations, various radiation related calculations have been conducted with reasonable caution for over 70 years. These calculations include: radiation shielding calculations, penetration assessments [1], and radiation dosimetry calculations [2]. Tying into the physics of interaction of radiation with matter, a beam which is comprised of either energetic neutrons or high-energy photons (i.e. gamma rays or hard X-rays) have conventionally been looked upon as candidates for interception via the standard X-section (i.e. cross-section) inspired models, for which one uses the typical expression

$$I(z) = I(0)e^{-\mu z} \quad (1)$$

where $z$ is the depth and $\mu$ is the attenuation parameter or 'coefficient'.

One could refer to this expression as the "simplistic" attenuation formula. Many shielding calculations done by medical physicists and health physicists over the years (since at least 1970) have been done using an extension of this "simplistic attenuation formula" coupled with a buildup factor [3]. The buildup factor is necessary if we are to use an attenuation formula as our principal tool of analysis to account for scattered neutrons or energetic photons, respectively [4]. Some of these "n's" and "γ's" are almost (or totally) elastically scattered, but some of these are down-scattered in energy. In keeping with the aspiration for excellence from the era of the “Space Age”, some computer/electrical engineers and nuclear engineers have enhanced these efforts by conducting Monte Carlo simulations of the transport and penetration of neutrons or photons through walls and various barriers with various respected codes (packages with X-section libraries) such as MCNP, EGS4, EGSNrc, and the versatile but highly tedious GEANT4 [5, 6, 7]. In this short paper, we introduce and present a deterministic and semi-analytical method for doing transport analysis on neutrons and isotropically scattering ‘hard’ photons which are placed in two energy and, with future ambitions, into 3 energy groups. There are advantages for doing such 2-group and higher multi-group analysis of radiative particles (i.e. neutrons and photons). These advantages are that we can more directly keep track of what percentage of radiative particles are close to the original high energy and how many are at significantly lower energy. An inspection of the profile of any build up function shows that the function is slightly larger than 1.0 at entry, then it rises to perhaps 2 or 3 within roughly one mean free path of the fast primary particles, and finally approaches the asymptote of 1.0 as the penetration depth gets progressively larger.

Review of Buildup and Discussion of Methods of Analysis

The attenuation formula in the introduction expresses the particulate intensity, not the energetic intensity in our convention. This corresponds to the choice of analysis of particulate flux rather than energetic flux of radiation. For those with a non-nuclear background: “Flux” is used by the medical physics community and health physicists with a definition considerably different from that of the flux of electric fields. Our flux has units of “particles” per cm$^2$ per second. See Frank Attix’s text [8] if this is unclear. Particulate intensity is less than or equal to the scalar flux of the particles. Very often intensity is defined as the magnitude of net current of transported particles per cm$^2$ of surface per sec. Indeed, in a case where equally many particles approach and penetrate a wall bidirectional the net current is zero. However, the scalar flux in such an example is much larger than zero. Admittedly, one can make some inferences on the approximate ratio of down-scattered particles at a given depth as a function of position by inspecting the buildup function, should it be available in published tables for a given shape and material. Here the intensity with a buildup coefficient can be expressed as

$$I(z) = B(E,z)I(0)e^{-\mu z}$$

However, the buildup coefficient (i.e. B(E, z)) does not make clear just what percentage of scattered $\gamma$-ray or X-ray radiation is scattered so as to retain most of its energy and how much has been “demoted” to photons with 50% less or more of energy per radiative particle. This reality (of radiative particles often undergoing elastic or nearly elastic scattering) holds for ‘free’ neutrons which scatter off of nuclei with an Atomic Number greater than 4. For example, when a ‘n’ with a kinetic energy of 1.0MeV collides with an Fe-56 nucleus, it has a 99% chance of undergoing elastic or quasielastic scattering (retaining most of its MeV). This same ‘n’ has a probability of less than 0.8% of down-scattering to a “low” energy neutron with less than 0.201MeV of kinetic energy. The chance of capture at 1MeV is less than one in a 1000. Fast neutrons are not easily captured, they usually are just scattered. This information is given on the basis of inspecting generally available public claims of two group data and, more importantly, by our careful conducting of MCNP [5] simulations in which we reproduce the conditions of scattering and so called “buildup” of scattered neutrons in rectangular walls of Iron (Fe-56). We are authorized to use and has extensive experience with the very versatile MCNP Monte Carlo code. Regarding neutrons, the buildup coefficient data either is not widely published for nuclear engineers or not readily available. Thus, the case of neutron shielding analysis offers a major service for performing multi-group energetic neutron flux and dose calculations. In this paper we stick with 2-group n’s (i.e. neutrons). The following formula is an overly simplistic & inadequate expression often used for metals with atomic number less than 84, away from the “uranic” family

$$\Phi(z) = \Phi(0)e^{-\sum [iotl,at,1MeV]^2}$$

However, Equation (3) gives a very incomplete story of the local neutron flux. Just as Equation (1) of I(z), which equals I(0)exp$^{-\mu z}$, gives a very incomplete story of local current/intensity of X-rays and $\gamma$-rays, when one fails to include the Buildup factor included in Equation (2), B(E, z) in Equation (2) is the buildup factor for the calculation of local ‘flux’ of $\gamma$-rays or photons at various depths (e.g. z) of penetration. It actually is more insightful to replace I(z) from Equation (1) with $\Phi(z)$ in Equation (3) where Flux is more appropriate than ‘n’ or photon current in the geometry of a wall or box since many of the neutrons or photons no longer travel straight forward along the z-axis after one or two collisions of scatter occur.

If the nuclei/atoms of a medium which is entered by the beam of n’s or photons is a pure absorber, then Equation (1) and Equation (3) are acceptable solutions for penetration and the differential equation which explains the transport of the particle(s) is given by:

$$\frac{d\left(I[2](z)\right)}{dz} + (\mu) \cdot I[2](z) = 0, \text{ for } \mu = \sigma_{aborb}^n \text{ rad}$$ and where $n_{\text{rad}}$ is the number density from radiation. If the medium is Boron-10 and the neutrons are at low energy, then Equation (4) would be a realistic equation for modeling transport of the neutrons, because B-10 is nearly a pure absorber. However, most materials are not pure absorbers.

If we presume that there are two energy levels for neutrons and photons (i.e. fast n's and slow n's), then it is appropriate to write a double energy-group Maxwell Boltzmann Transport Equation (MBTE) in order to express

$$\bar{n}(\theta) \cdot \nabla \Psi_{[2]}(\bar{r}, \theta) + (\Sigma_{ab,2} + \Sigma_{s,2} + \Sigma_{s[2,1]}) \cdot \Psi_{[2]}(\bar{r}, \theta) = \oint \Sigma_{s,2} \cdot \Psi_{[2]}(\bar{r}, \theta') \cdot (jacobi(\mathcal{L})) d\theta'$$

$$\bar{n}(\theta) \cdot \nabla \Psi_{[1]}(\bar{r}, \theta) + (\Sigma_{ab,1} + \Sigma_{s,1}) \cdot \Psi_{[1]}(\bar{r}, \theta) = \oint \Sigma_{s,1} \cdot \Psi_{[1]}(\bar{r}, \theta') + \Sigma_{s[2,1]} \cdot \Psi_{[2]}(\bar{r}, \theta')$$

Suppose that $\Sigma_{s[2, 1]}$ and $\Sigma_{s[2, 2]}$ are equal to zero, as we might imagine for a medium of "super-Boron". Then, Equation (5a) can be rewritten as:

$$\bar{n} \cdot \nabla \Psi_{[2]}(\bar{r}, _{-}) + (\Sigma_{ab,2} + 0) \cdot \Psi_{[2]}(\bar{r}, _{-}) = 0.$$

If our wall is very broad and if the distribution to approximation depends only on coordinate $z$, then this simple 1st order equation is completely equivalent to Equation (4) on attenuation on the 2nd page. In this paper, we presume that scattering is isotropic, which is often a good approximation for the scattering of neutrons. Therefore, we establish that $\Sigma_{s[2, 1]}$ and $\Sigma_{s[2, 2]}$ are completely independent of $\theta, \theta'$, and any angle. This 2-group version of the MBTE is an example of a pair of integrodifferential equations. Generally, it is easier to solve a purely integral equation, such as a Fredholm Int. Equation (12). Holding on to the presumption of isotropic scattering, Equations (5a) and (5b) can be subjected to a special integral transformation via Green's functions in order to re-express them as the following integral equations.

Equations (6a) & (6b) are challenging, but some solutions have been found. In the Russian Math academy of the 1950s and by a 'computations' group at Los Alamos in the 1950s some solutions have been found to Equation (6). The author(s) have found a method for numerically with arbitrary precision to iteratively solve the monoenergetic version of eqn.(6a) and the 2-group Equation (6) [13].

Equation 6 is more demanding and difficult to solve than it is to simplistically use Equation (2) with the buildup coefficient to calculate relative intensities. The information which we get out of Equations (6) for $\Phi_{[2]}$ and $\Phi_{[1]}$ is much richer than what we can get from the calculation of I(z) as a function of penetration into a wall Equation(2). This is especially apparent for of very broad slab of shielding (presuming breadth of slab is more than 6 times > than thickness of slab). Presuming a beam of neutrons enters from the left at the interface where $z=0$, the user would need to read, download, and interpolate a table of buildup coefficients values (or crude local formula) for a slab of the given material (such as iron, Pb, or concrete). These tables are almost non-existent for neutrons. Such tables do exist for $\gamma$-ray and X-ray photons for industrial and medical materials, but they are limited in their range of and diversity of examples. Thus, the dosimetrist who is borrowing or using the data often is stuck with having to interpolate from near fits of other examples most similar to the geometry which he or she has chosen to design or assess for predictions or dose verifications. As a reminder of geometrical concepts, the portion of "battered" particulate current density which escapes from the right-hand boundary of a rectangular slab of shielding from "mid-face" equals approximately $\frac{1}{2}$ or

0.6 times the "battered" scalar flux which is present on the boundary of escape from the rectangularly shaped shield. Battered here is the condition that a particle has never been coerced or forced to scatter. "Battering" a neutron changes its direction and can reduce its energy. For a very thin shield whose thickness is less than 1/4 of a mean free path of a neutron, the first author has verified that $I_{[battered]}(escape) = \Phi_{[battered]}(boundary)$ by analytical means of integrating the Green's function of the radiative source. For thicker samples MCNP simulations have verified with converging persuasiveness that the escape ratio is in the range of 1/2 through 0.65, depending on thickness of shielding. $I_{[battered]}$ refers to the particulate current density of 'n's which have undergone at least one collision before escaping from either the left-hand side (LHS) or right-hand side (RHS) of the slab. Many books casually refer to 'I' in $I_{[battered]}$ and $I_{[entering_beam]}$ as intensity. In this paper, 'I' shall mean current density of particles (mostly neutrons here) per cm$^2$ (or m$^2$) per second. If the beam and scattered particle are all mono-energetic, then, of course, energetic I, or energetic intensity, is given by

$$\Phi_{[2][n]}(r) = I(0) \cdot e^{-\left( \sum_{total,1MeV} \frac{z}{\cos(\theta)} \right)} + \Phi_{[2][1]}(z) + \Phi_{[2][2]}(z) + \Phi_{[2][3]}(z) + \dots$$

where $\Phi_{[2][n]}(z)$ is generated iteratively from the input of $\Phi_{[2][n-1]}(z)$ into the RHS of Equation(6). The writing out of the summation of $\Phi_{[2][n]}(z)$ with respect to iteration index $n$ is akin to the early methods of perturbative quantum electrodynamics done by Hans Bethe circa [15], but only without the need to subtract infinities out. Renormalization methods are not required to solve Eqs.(6), analytically or computationally. The first iteration, $\Phi_{[2][1]}(z)$, contains logarithms of $z$, Ei functions of $z$, and terms of exp($-\Sigma_z z$) for several factors. The "2" refers to the 1-MeV neutrons and can be efficiently approximated by a much more tractable sum of a polynomials and log expressions. $\Phi_{[2][1]}(z)$ is the flux of fast neutrons which have been scattered only once. We allow "batter-2" to refer to radiative particles which have undergone collisions (preferably) in highest group only once. Following thru on the iteration, $\Phi_{[n][3]}(z)$ is the flux of neutrons which have been battered once, where the $n$ equals either 1 or 2. Likewise, $\Phi_{[m][3]}(z)$ is the flux of fast neutrons which have been battered 2 times (i.e. 3-1 times). 'Fast' albedo is the ratio of the sum fast scattering cross sections divided by the total fast cross section of neutrons. So long as the albedo is less than or equal to 1, the series of $\Phi_{[m][m]}(z)$ is guaranteed to converge. Five to eight iterations have proven to be allow for sufficiently convergent solutions. We use the approximations that the current of escaping neutrons which penetrated = 'unbattered' flux + 0.5 (sum of fluxes comprised of scattered neutrons). 0.5 is the analytically and geometrically guaranteed minimum. Judiciously, we occasionally replace 0.5 with 0.52 or even 0.55. This is the unit-less factor of escape which relates surface $\Phi$ value to current density of scattered neutrons which depart a surface.

For the sake of space and the context of this paper, we will focus on predictions for scattering deterministically of the transmission of neutrons through a sample slabs iron at various thicknesses and on the prediction for the portion of neutrons which are returned backwards (due to back-scatter) from the slab. Our deterministic method distinguishes between the population of fast scattered neutrons and of slow neutrons (where the '1' neutrons are the designated down-scattered (i.e. slow) neutrons). We also conducted Monte Carlo simulations of neutrons from a beam which approach the same wall of Iron material. The two Monte Carlo (MC) codes used are: MCNP (which was developed and updated by LANL) and a homemade Monte Carlo code which was developed early in 2014 [13] and had proven to be valid for the modeling of isotropic scattering. In Table 1, just like in the paper of 2014, the name assigned to our homemade MC code is SMUSKE, where SMUSKE is designed to simulate isotropic scattering or radiative particles in arbitrarily chosen rectangular geometries.

Two examples are given below for a broad rectangular slab of homogeneous cast-iron where a beam of neutrons approaches the barrier at the normal angle. Our deterministic predictions, the predictions of MCNP, and the predictions of the MC code, SMUSKE, are included in Table 1. In the common style, MFP shall denote Mean-Free-Path of neutron in response to the total cross section for the fast (1MeV to .9MeV) neutrons. Our deterministic algorithm is named “IntegIterator”, which stands for Iterative Integral Equation Solver. “Deterministic Iterator” is a slightly lengthy alias for our abbreviation of IntegIterator. As a reminder, Equations (6a) & (6b) are Fredholm integral equations in terms of mathematical structure. “mfpm” in column one means the value of MFP/(Thickness_wall). “Grp#’ in the 2nd column of table refers to the energy group number of the neutrons. Grp# 2 includes all neutrons which are in the energy range of 0.9MeV through 1.0MeV. These are the designated fast neutrons. Grp# 1 denotes the slow neutrons, which include all of the free neutrons which have energy in the range of 0 through 0.160MeV. There were less than One in 5000 neutrons found in the energy range which lies in between the ‘slow’ Group and the ‘fast’ group of neutrons. Therefore, we do not include such an intermediate group of negligible population.

In the first half of Table (1), thickness of the rectangular slab is 1 MFP long, which is 2.884cm. According to SMUSKE, out of 10,000 incoming fast neutrons, 3132.8 fast n’s travel out or backward from the wall, and 106.5 slow n’s escape back toward the source of the beam. Accordingly, SMUSKE predicts that 66917 fast n’s escape forwards out through the iron wall, and 63.7 slow n’s escape forward. Now we point out the level of agreement between our Deterministic Iterator and SMUSKE’s predictions. We review the thick wall with 2.884cm first. For collective back scattered n’s (which are to escape out of wall backward), the % differences between the Deterministic Iterator and SMUSKE are: 25.3% for slow n’s (i.e. of grp.1) and 14.6% for fast n’s (i.e. of grp.2). For the rate of forward escape, or transmission, the % differences between IntegIterator and SMUSKE are: 63% for slow n’s (i.e. of grp.1) and -4.85% for fast n’s (i.e. of grp.2). We briefly consider the thin wall of 1.442cm now: For rate of forward escape, or transmission, the percentage of disagreement between our IntegIterator and SMUSKE are: -30.2% for slow n’s and 0.30% for fast n’s (i.e. of grp.2). Let us compare the #s of SMUSKE to those of MCNP: For rate of backward escape, returning to beam source, when the wall is one MFP thick, the percentage of disagreement between the predictions of MCNP and SMUSKE are: 133.45% for slow n’s and −120% for fast n’s. For rate of forward escape, when the wall is one MFP thick, the percentage of disagreement between the predictions of MCNP and SMUSKE are: 154.84% for slow n’s and 19.80% for fast n’s.

For rate of backward escape, returning to beam source, when the wall is ½ MFP thick, the percentage of disagreement between the predictions of MCNP and SMUSKE are: 154.3% for slow n’s and -85.27% for fast n’s. For rate of forward escape, when the wall is ½ MFP thick, the percentage of disagreement between the famous MCNP modeler and SMUSKE are: 158.1% for slow n’s and 3.74% for fast n’s. On the other hand, on a more impressive note, there is a -12.29% disagreement between the respective predictions of our Deterministic Iterator (i.e. IntegIterator) and those of MCNP for the number of forward transmitted fast neutrons for the of wall which has 1MFP (2.884cm) . Also, on a more impressive note, there is a -1.71879% disagreement between the respective predictions of IntegIterator and MCNP for the number of forward transmitted neutrons for the of wall with thickness of ½ MFP (1.442cm). One can see that our deterministic IntegIterator makes predictions which are closer in agreement to those of MCNP than SMUSKE makes. We observe this in spite of the fact that SMUSKE and IntegIterator are at the disadvantage by presuming that all of the ‘n’-scattering cross sections are isotropic. MCNP does consider angular probabilities to great detail for ‘n’- scattering for all of the well-known isotopes on the period table of nuclides.

Indeed, the compilation of the ‘scatter-kernel’ data of MCNP was a project which spanned more than a decade. Thus, it would be very difficult to summarize such a large amount angular data with approximations of the zeroeth, first, and second Legendre` polynomials of the cosine of scatter-angle with any tractable and manageable database which could function without incurring “data strangulation” of an analytical iterator (e.g. our IntegIterator) written in Maplesoft, high-level Python, or similar language/package. Some ask “why not just give the ‘jobs’ to GEANT4 to execute?”. With all due acknowledgement of the formidable abilities of GEANT4, the three smaller codes SMUSKE, MCNP, and IntegIterator are all faster and easier to work with than GEANT4 for designated rectangular walls of metal bombarded by n’s.

Conclusion

Our deterministic iterative algorithm, IntegIterator, agrees reasonably well in regard to the prediction of the energy distribution and the direction distribution to the corresponding distributions of energy and way from our simulations done via MCNP, in spite of the deficiency of Integiterator not being able to process anisotropic scattering yet. This can be seen from the results posted in Table I and in the summary of Table I above. On the other hand, neither our deterministic code nor our isotropically designed SMUSKE agree extremely well with the predictions of back- scatter of MCNP. All three methods of calculation agree well in the prediction of percent neutrons captured in the iron wall (which end ups being a very small percentage loss at the thickness of 1 or 2 cm), and in the overall statistic of neutron escape or “leakage”. However, the directionality and ratio of down-scattering is subject to disagreement. As suggested above, MCNP has extremely detailed X-section libraries. Iron-56 turns out to be one of these extremely an isotropically scattering isotopes. With the exception of the work of Chandrasekhar [16] and his use of H-functions for tracing intensities of scattered photons, there is almost zero analytical work recorded on predicting flux densities theoretically which are solutions to the MBTE in which the scattering X-sections are anisotropic. It is “semi-easy” to do analysis of solutions of the MBTE when the neutron scatterers (i.e. the nuclei) are isotropic. Many experts of reactor physics and shielding analysis approximate the multi-group MBTE as a diffusion equation of radiative particles. By its intrinsic nature, it is virtually impossible to do angle dependent “ray track tracing” of neutrons if one does modeling with a diffusion equation rather than the MBTE.

It would be convenient from a clinical radiological treatment planners’ point of view to carefully look up the data for B(E,z), where B(E,z) is the buildup factor included in Equation(2). However, in regard to neutrons, buildup coefficient data either is not widely published or is not available to the broad national/international communities of health physics and engineering. Moreover, a significant benefit of our deterministic method (and algorithm) of IntegIterator is the superior speed which it offers by its retention of the definitions of chosen geometries per wall to be bombarded and in its output extraction compared to the time and duties required at the conclusion of a corresponding run of an MCNP simulation. And GEANT4 is even more tedious than MCNP. Thus our 2-group IntegIterator algorithm and formulation is much faster than MCNP when predicting penetration ratios as well as the distribution of energy of penetrating neutrons.

The simulation of MCNP is sufficiently fast for modeling transmissions and down-scattering of neutrons through rectangular walls. However, the processing of the output data from the output files generated by the MCNP code require considerable data processing which is done best either in a UNIX console environment or a DOS console environment. Many of the younger physical engineers have neglected LINUX training and thus tend to rely on a Windows environment or a ‘Mac-Windows’ environment to process outputs of their chosen software modelers within Windows, ‘Mac-Windows’, or XWindows (if using LINUX version of Maplesoft). Maplesoft is the language which we have selected for IntegIterator in order to conduct or local Flux calculations. Our ‘Maple’ version of IntegiIterator can operate within the environments of Windows, XWindows of Mac OS, and Linux - as valid versions of Maplesoft can be placed in these OS’s. Much of our Maplesoft code can be translated into Matlab code, for the potential accommodation of electrical engineers (who gravitate to Matlab). Another, great benefit of our deterministic code is that one can generate a polynomial approximation of the local dose of neutron flux at any depth within the metal. With MCNP such a feat would require writing of an Input file which is more than ten-fold more elaborate than the input file for the IntegIterator code for the same slab of metallic material.

It is reasonable to anticipate a future effort of 3-group neutron flux distributions w.r.t. energy by us or radiation physics colleagues. However, for now in the 2020/21 academic year, we focus on constructing 2-group databases of neutron X-sections of various important materials besides just iron and boron and subsequently using the algorithm and code(s) of Integiterator to predict collective forward escape and backward escape of neutrons which initially enter slabs of the respective materials of interest.