# 5. Neutron Physics¶

There are limited differences between physics treatments used in the continuous-energy and multi-group modes. If distinctions are necessary, each of the following sections will provide an explanation of the differences. Otherwise, replacing any references of the particle’s energy (E) with references to the particle’s energy group (g) will suffice.

## 5.1. Sampling Distance to Next Collision¶

As a particle travels through a homogeneous material, the probability distribution function for the distance to its next collision \(\ell\) is

where \(\Sigma_t\) is the total macroscopic cross section of the material. Equation (1) tells us that the further the distance is to the next collision, the less likely the particle will travel that distance. In order to sample the probability distribution function, we first need to convert it to a cumulative distribution function

By setting the cumulative distribution function equal to \(\xi\), a random number on the unit interval, and solving for the distance \(\ell\), we obtain a formula for sampling the distance to next collision:

Since \(\xi\) is uniformly distributed on \([0,1)\), this implies that \(1 - \xi\) is also uniformly distributed on \([0,1)\) as well. Thus, the formula usually used to calculate the distance to next collision is

## 5.2. \((n,\gamma)\) and Other Disappearance Reactions¶

All absorption reactions other than fission do not produce any secondary neutrons. As a result, these are the easiest type of reactions to handle. When a collision occurs, the first step is to sample a nuclide within a material. Once the nuclide has been sampled, then a specific reaction for that nuclide is sampled. Since the total absorption cross section is pre-calculated at the beginning of a simulation, the first step in sampling a reaction is to determine whether a “disappearance” reaction occurs where no secondary neutrons are produced. This is done by sampling a random number \(\xi\) on the interval \([0,1)\) and checking whether

where \(\sigma_t\) is the total cross section, \(\sigma_a\) is the absorption cross section (this includes fission), and \(\sigma_f\) is the total fission cross section. If this condition is met, then the neutron is killed and we proceed to simulate the next neutron from the source bank.

Note that photons arising from \((n,\gamma)\) and other neutron reactions are not produced in a microscopically correct manner. Instead, photons are sampled probabilistically at each neutron collision, regardless of what reaction actually takes place. This is described in more detail in Photon Production.

## 5.3. Elastic Scattering¶

Note that the multi-group mode makes no distinction between elastic or inelastic scattering reactions. The specific multi-group scattering implementation is discussed in the Multi-Group Scattering section.

Elastic scattering refers to the process by which a neutron scatters off a nucleus and does not leave it in an excited state. It is referred to as “elastic” because in the center-of-mass system, the neutron does not actually lose energy. However, in lab coordinates, the neutron does indeed lose energy. Elastic scattering can be treated exactly in a Monte Carlo code thanks to its simplicity.

Let us discuss how OpenMC handles two-body elastic scattering kinematics. The first step is to determine whether the target nucleus has any associated motion. Above a certain energy threshold (400 kT by default), all scattering is assumed to take place with the target at rest. Below this threshold though, we must account for the thermal motion of the target nucleus. Methods to sample the velocity of the target nucleus are described later in section Effect of Thermal Motion on Cross Sections. For the time being, let us assume that we have sampled the target velocity \(\mathbf{v}_t\). The velocity of the center-of-mass system is calculated as

where \(\mathbf{v}_n\) is the velocity of the neutron and \(A\) is the
atomic mass of the target nucleus measured in neutron masses (commonly referred
to as the *atomic weight ratio*). With the velocity of the center-of-mass
calculated, we can then determine the neutron’s velocity in the center-of-mass
system:

where we have used uppercase \(\mathbf{V}\) to denote the center-of-mass system. The direction of the neutron in the center-of-mass system is

At low energies, elastic scattering will be isotropic in the center-of-mass system, but for higher energies, there may be p-wave and higher order scattering that leads to anisotropic scattering. Thus, in general, we need to sample a cosine of the scattering angle which we will refer to as \(\mu\). For elastic scattering, the secondary angle distribution is always given in the center-of-mass system and is sampled according to the procedure outlined in Sampling Angular Distributions. After the cosine of the angle of scattering has been sampled, we need to determine the neutron’s new direction \(\mathbf{\Omega}'_n\) in the center-of-mass system. This is done with the procedure in Transforming a Particle’s Coordinates. The new direction is multiplied by the speed of the neutron in the center-of-mass system to obtain the new velocity vector in the center-of-mass:

Finally, we transform the velocity in the center-of-mass system back to lab coordinates:

In OpenMC, the angle and energy of the neutron are stored rather than the velocity vector itself, so the post-collision angle and energy can be inferred from the post-collision velocity of the neutron in the lab system.

For tallies that require the scattering cosine, it is important to store the scattering cosine in the lab system. If we know the scattering cosine in the center-of-mass, the scattering cosine in the lab system can be calculated as

However, equation (11) is only valid if the target was at rest. When the target nucleus does have thermal motion, the cosine of the scattering angle can be determined by simply taking the dot product of the neutron’s initial and final direction in the lab system.

## 5.4. Inelastic Scattering¶

Note that the multi-group mode makes no distinction between elastic or inelastic scattering reactions. The specific multi-group scattering implementation is discussed in the Multi-Group Scattering section.

The major algorithms for inelastic scattering were described in previous sections. First, a scattering cosine is sampled using the algorithms in Sampling Angular Distributions. Then an outgoing energy is sampled using the algorithms in Sampling Energy Distributions. If the outgoing energy and scattering cosine were given in the center-of-mass system, they are transformed to laboratory coordinates using the algorithm described in Transforming a Particle’s Coordinates. Finally, the direction of the particle is changed also using the procedure in Transforming a Particle’s Coordinates.

Although inelastic scattering leaves the target nucleus in an excited state, no secondary photons from nuclear de-excitation are tracked in OpenMC.

## 5.5. \((n,xn)\) Reactions¶

Note that the multi-group mode makes no distinction between elastic or inelastic scattering reactions. The specific multi-group scattering implementation is discussed in the Multi-Group Scattering section.

These types of reactions are just treated as inelastic scattering and as such are subject to the same procedure as described in Inelastic Scattering. For reactions with integral multiplicity, e.g., \((n,2n)\), an appropriate number of secondary neutrons are created. For reactions that have a multiplicity given as a function of the incoming neutron energy (which occasionally occurs for MT=5), the weight of the outgoing neutron is multiplied by the multiplicity.

## 5.6. Multi-Group Scattering¶

In multi-group mode, a scattering collision requires that the outgoing energy group of the simulated particle be selected from a probability distribution, the change-in-angle selected from a probability distribution according to the outgoing energy group, and finally the particle’s weight adjusted again according to the outgoing energy group.

The first step in selecting an outgoing energy group for a particle in a given incoming energy group is to select a random number (\(\xi\)) between 0 and 1. This number is then compared to the cumulative distribution function produced from the outgoing group (g’) data for the given incoming group (g):

If the scattering data is represented as a Legendre expansion, then the value of \(\Sigma_{s,g \rightarrow g'}\) above is the 0th order for the given group transfer. If the data is provided as tabular or histogram data, then \(\Sigma_{s,g \rightarrow g'}\) is the sum of all bins of data for a given g and g’ pair.

Now that the outgoing energy is known the change-in-angle, \(\mu\) can be determined. If the data is provided as a Legendre expansion, this is done by rejection sampling of the probability distribution represented by the Legendre series. For efficiency, the selected values of the PDF (\(f(\mu)\)) are chosen to be between 0 and the maximum value of \(f(\mu)\) in the domain of -1 to 1. Note that this sampling scheme automatically forces negative values of the \(f(\mu)\) probability distribution function to be treated as zero probabilities.

If the angular data is instead provided as a tabular representation, then the value of \(\mu\) is selected as described in the Tabular Angular Distribution section with a linear-linear interpolation scheme.

If the angular data is provided as a histogram representation, then the value of \(\mu\) is selected in a similar fashion to that described for the selection of the outgoing energy (since the energy group representation is simply a histogram representation) except the CDF is composed of the angular bins and not the energy groups. However, since we are interested in a specific value of \(\mu\) instead of a group, then an angle is selected from a uniform distribution within from the chosen angular bin.

The final step in the scattering treatment is to adjust the weight of the neutron to account for any production of neutrons due to \((n,xn)\) reactions. This data is obtained from the multiplicity data provided in the multi-group cross section library for the material of interest. The scaled value will default to 1.0 if no value is provided in the library.

## 5.7. Fission¶

While fission is normally considered an absorption reaction, as far as it concerns a Monte Carlo simulation it actually bears more similarities to inelastic scattering since fission results in secondary neutrons in the exit channel. Other absorption reactions like \((n,\gamma)\) or \((n,\alpha)\), on the contrary, produce no neutrons. There are a few other idiosyncrasies in treating fission. In an eigenvalue calculation, secondary neutrons from fission are only “banked” for use in the next generation rather than being tracked as secondary neutrons from elastic and inelastic scattering would be. On top of this, fission is sometimes broken into first-chance fission, second-chance fission, etc. The nuclear data file either lists the partial fission reactions with secondary energy distributions for each one, or a total fission reaction with a single secondary energy distribution.

When a fission reaction is sampled in OpenMC (either total fission or, if data exists, first- or second-chance fission), the following algorithm is used to create and store fission sites for the following generation. First, the average number of prompt and delayed neutrons must be determined to decide whether the secondary neutrons will be prompt or delayed. This is important because delayed neutrons have a markedly different spectrum from prompt neutrons, one that has a lower average energy of emission. The total number of neutrons emitted \(\nu_t\) is given as a function of incident energy in the ENDF format. Two representations exist for \(\nu_t\). The first is a polynomial of order \(N\) with coefficients \(c_0,c_1,\dots,c_N\). If \(\nu_t\) has this format, we can evaluate it at incoming energy \(E\) by using the equation

The other representation is just a tabulated function with a specified interpolation law. The number of prompt neutrons released per fission event \(\nu_p\) is also given as a function of incident energy and can be specified in a polynomial or tabular format. The number of delayed neutrons released per fission event \(\nu_d\) can only be specified in a tabular format. In practice, we only need to determine \(nu_t\) and \(nu_d\). Once these have been determined, we can calculated the delayed neutron fraction

We then need to determine how many total neutrons should be emitted from fission. If no survival biasing is being used, then the number of neutrons emitted is

where \(w\) is the statistical weight and \(k_{eff}\) is the effective multiplication factor from the previous generation. The number of neutrons produced is biased in this manner so that the expected number of fission neutrons produced is the number of source particles that we started with in the generation. Since \(\nu\) is not an integer, we use the following procedure to obtain an integral number of fission neutrons to produce. If \(\xi > \nu - \lfloor \nu \rfloor\), then we produce \(\lfloor \nu \rfloor\) neutrons. Otherwise, we produce \(\lfloor \nu \rfloor + 1\) neutrons. Then, for each fission site produced, we sample the outgoing angle and energy according to the algorithms given in Sampling Angular Distributions and Sampling Energy Distributions respectively. If the neutron is to be born delayed, then there is an extra step of sampling a delayed neutron precursor group since they each have an associated secondary energy distribution.

The sampled outgoing angle and energy of fission neutrons along with the position of the collision site are stored in an array called the fission bank. In a subsequent generation, these fission bank sites are used as starting source sites.

The above description is similar for the multi-group mode except the data are provided as group-wise data instead of in a continuous-energy format. In this case, the outgoing energy of the fission neutrons are represented as histograms by way of either the nu-fission matrix or chi vector.

## 5.8. Secondary Angle-Energy Distributions¶

Note that this section is specific to continuous-energy mode since the multi-group scattering process has already been described including the secondary energy and angle sampling.

For a reaction with secondary products, it is necessary to determine the outgoing angle and energy of the products. For any reaction other than elastic and level inelastic scattering, the outgoing energy must be determined based on tabulated or parameterized data. The ENDF-6 Format specifies a variety of ways that the secondary energy distribution can be represented. ENDF File 5 contains uncorrelated energy distribution whereas ENDF File 6 contains correlated energy-angle distributions. The ACE format specifies its own representations based loosely on the formats given in ENDF-6. OpenMC’s HDF5 nuclear data files use a combination of ENDF and ACE distributions; in this section, we will describe how the outgoing angle and energy of secondary particles are sampled.

One of the subtleties in the nuclear data format is the fact that a single reaction product can have multiple angle-energy distributions. This is mainly useful for reactions with multiple products of the same type in the exit channel such as \((n,2n)\) or \((n,3n)\). In these types of reactions, each neutron is emitted corresponding to a different excitation level of the compound nucleus, and thus in general the neutrons will originate from different energy distributions. If multiple angle-energy distributions are present, they are assigned incoming-energy-dependent probabilities that can then be used to randomly select one.

Once a distribution has been selected, the procedure for determining the outgoing angle and energy will depend on the type of the distribution.

### 5.8.2. Product Angle-Energy Distributions¶

If the secondary distribution for a product was given in file 6 in ENDF, the angle and energy are correlated with one another and cannot be sampled separately. Several representations exist in ENDF/ACE for correlated angle-energy distributions.

#### 5.8.2.3. N-Body Phase Space Distribution¶

Reactions in which there are more than two products of similar masses are sometimes best treated by using what’s known as an N-body phase distribution. This distribution has the following probability density function for outgoing energy and angle of the \(i\)-th particle in the center-of-mass system:

where \(n\) is the number of outgoing particles, \(C_n\) is a normalization constant, \(E_i^{max}\) is the maximum center-of-mass energy for particle \(i\), and \(E'\) is the outgoing energy. We see in equation (47) that the angle is simply isotropic in the center-of-mass system. The algorithm for sampling the outgoing energy is based on algorithms R28, C45, and C64 in the Monte Carlo Sampler. First we calculate the maximum energy in the center-of-mass using the following equation:

where \(A_p\) is the total mass of the outgoing particles in neutron masses, \(A\) is the mass of the original target nucleus in neutron masses, and \(Q\) is the Q-value of the reaction. Next we sample a value \(x\) from a Maxwell distribution with a nuclear temperature of one using the algorithm outlined in Maxwell Fission Spectrum. We then need to determine a value \(y\) that will depend on how many outgoing particles there are. For \(n = 3\), we simply sample another Maxwell distribution with unity nuclear temperature. For \(n = 4\), we use the equation

where \(\xi_i\) are random numbers sampled on the interval \([0,1)\). For \(n = 5\), we use the equation

After \(x\) and \(y\) have been determined, the outgoing energy is then calculated as

There are two important notes to make regarding the N-body phase space distribution. First, the documentation (and code) for MCNP5-1.60 has a mistake in the algorithm for \(n = 4\). That being said, there are no existing nuclear data evaluations which use an N-body phase space distribution with \(n = 4\), so the error would not affect any calculations. In the ENDF/B-VII.1 nuclear data evaluation, only one reaction uses an N-body phase space distribution at all, the \((n,2n)\) reaction with H-2.

## 5.9. Transforming a Particle’s Coordinates¶

Since all the multi-group data exists in the laboratory frame of reference, this section does not apply to the multi-group mode.

Once the cosine of the scattering angle \(\mu\) has been sampled either from a angle distribution or a correlated angle-energy distribution, we are still left with the task of transforming the particle’s coordinates. If the outgoing energy and scattering cosine were given in the center-of-mass system, then we first need to transform these into the laboratory system. The relationship between the outgoing energy in center-of-mass and laboratory is

where \(E'_{cm}\) is the outgoing energy in the center-of-mass system, \(\mu_{cm}\) is the scattering cosine in the center-of-mass system, \(E'\) is the outgoing energy in the laboratory system, and \(E\) is the incident neutron energy. The relationship between the scattering cosine in center-of-mass and laboratory is

where \(\mu\) is the scattering cosine in the laboratory system. The scattering cosine still only tells us the cosine of the angle between the original direction of the particle and the new direction of the particle. If we express the pre-collision direction of the particle as \(\mathbf{\Omega} = (u,v,w)\) and the post-collision direction of the particle as \(\mathbf{\Omega}' = (u',v',w')\), it is possible to relate the pre- and post-collision components. We first need to uniformly sample an azimuthal angle \(\phi\) in \([0, 2\pi)\). After the azimuthal angle has been sampled, the post-collision direction is calculated as

## 5.10. Effect of Thermal Motion on Cross Sections¶

Since all the multi-group data should be generated with thermal scattering treatments already, this section does not apply to the multi-group mode.

When a neutron scatters off of a nucleus, it may often be assumed that the target nucleus is at rest. However, the target nucleus will have motion associated with its thermal vibration, even at absolute zero (This is due to the zero-point energy arising from quantum mechanical considerations). Thus, the velocity of the neutron relative to the target nucleus is in general not the same as the velocity of the neutron entering the collision.

The effect of the thermal motion on the interaction probability can be written as

where \(v_n\) is the magnitude of the velocity of the neutron, \(\bar{\sigma}\) is an effective cross section, \(T\) is the temperature of the target material, \(\mathbf{v}_T\) is the velocity of the target nucleus, \(v_r = || \mathbf{v}_n - \mathbf{v}_T ||\) is the magnitude of the relative velocity, \(\sigma\) is the cross section at 0 K, and \(M (\mathbf{v}_T)\) is the probability distribution for the target nucleus velocity at temperature \(T\) (a Maxwellian). In a Monte Carlo code, one must account for the effect of the thermal motion on both the integrated cross section as well as secondary angle and energy distributions. For integrated cross sections, it is possible to calculate thermally-averaged cross sections by applying a kernel Doppler broadening algorithm to data at 0 K (or some temperature lower than the desired temperature). The most ubiquitous algorithm for this purpose is the SIGMA1 method developed by Red Cullen and subsequently refined by others. This method is used in the NJOY and PREPRO data processing codes.

The effect of thermal motion on secondary angle and energy distributions can be accounted for on-the-fly in a Monte Carlo simulation. We must first qualify where it is actually used however. All threshold reactions are treated as being independent of temperature, and therefore they are not Doppler broadened in NJOY and no special procedure is used to adjust the secondary angle and energy distributions. The only non-threshold reactions with secondary neutrons are elastic scattering and fission. For fission, it is assumed that the neutrons are emitted isotropically (this is not strictly true, but is nevertheless a good approximation). This leaves only elastic scattering that needs a special thermal treatment for secondary distributions.

Fortunately, it is possible to directly sample the velocity of the target nuclide and then use it directly in the kinematic calculations. However, this calculation is a bit more nuanced than it might seem at first glance. One might be tempted to simply sample a Maxwellian distribution for the velocity of the target nuclide. Careful inspection of equation (55) however tells us that target velocities that produce relative velocities which correspond to high cross sections will have a greater contribution to the effective reaction rate. This is most important when the velocity of the incoming neutron is close to a resonance. For example, if the neutron’s velocity corresponds to a trough in a resonance elastic scattering cross section, a very small target velocity can cause the relative velocity to correspond to the peak of the resonance, thus making a disproportionate contribution to the reaction rate. The conclusion is that if we are to sample a target velocity in the Monte Carlo code, it must be done in such a way that preserves the thermally-averaged reaction rate as per equation (55).

The method by which most Monte Carlo codes sample the target velocity for use in elastic scattering kinematics is outlined in detail by [Gelbard]. The derivation here largely follows that of Gelbard. Let us first write the reaction rate as a function of the velocity of the target nucleus:

where \(R\) is the reaction rate. Note that this is just the right-hand side of equation (55). Based on the discussion above, we want to construct a probability distribution function for sampling the target velocity to preserve the reaction rate – this is different from the overall probability distribution function for the target velocity, \(M ( \mathbf{v}_T )\). This probability distribution function can be found by integrating equation (56) to obtain a normalization factor:

Let us call the normalization factor in the denominator of equation (57) \(C\).

### 5.10.1. Constant Cross Section Model¶

It is often assumed that \(\sigma (v_r)\) is constant over the range of relative velocities of interest. This is a good assumption for almost all cases since the elastic scattering cross section varies slowly with velocity for light nuclei, and for heavy nuclei where large variations can occur due to resonance scattering, the moderating effect is rather small. Nonetheless, this assumption may cause incorrect answers in systems with low-lying resonances that can cause a significant amount of up-scatter that would be ignored by this assumption (e.g. U-238 in commercial light-water reactors). We will revisit this assumption later in Energy-Dependent Cross Section Model. For now, continuing with the assumption, we write \(\sigma (v_r) = \sigma_s\) which simplifies (57) to

The Maxwellian distribution in velocity is

where \(m\) is the mass of the target nucleus and \(k\) is Boltzmann’s constant. Notice here that the term in the exponential is dependent only on the speed of the target, not on the actual direction. Thus, we can change the Maxwellian into a distribution for speed rather than velocity. The differential element of velocity is

Let us define the Maxwellian distribution in speed as

To simplify things a bit, we’ll define a parameter

Substituting equation (62) into equation (61), we obtain

Now, changing variables in equation (58) by using the result from equation (61), our new probability distribution function is

Again, the Maxwellian distribution for the speed of the target nucleus has no dependence on the angle between the neutron and target velocity vectors. Thus, only the term \(|| \mathbf{v}_n - \mathbf{v}_T ||\) imposes any constraint on the allowed angle. Our last task is to take that term and write it in terms of magnitudes of the velocity vectors and the angle rather than the vectors themselves. We can establish this relation based on the law of cosines which tells us that

Thus, we can infer that

Inserting equation (66) into (64), we obtain

This expression is still quite formidable and does not lend itself to any natural sampling scheme. We can divide this probability distribution into two parts as such:

In general, any probability distribution function of the form \(p(x) = f_1(x) f_2(x)\) with \(f_1(x)\) bounded can be sampled by sampling \(x'\) from the distribution

and accepting it with probability

The reason for dividing and multiplying the terms by \(v_n + v_T\) is to ensure that the first term is bounded. In general, \(|| \mathbf{v}_n - \mathbf{v}_T ||\) can take on arbitrarily large values, but if we divide it by its maximum value \(v_n + v_T\), then it ensures that the function will be bounded. We now must come up with a sampling scheme for equation (69). To determine \(q(v_T)\), we need to integrate \(f_2\) in equation (68). Doing so we find that

Thus, we need to sample the probability distribution function

Now, let us do a change of variables with the following definitions

Substituting equation (73) into equation (72) along with \(dx = \beta dv_T\) and doing some crafty rearranging of terms yields

It’s important to make note of the following two facts. First, the terms outside the parentheses are properly normalized probability distribution functions that can be sampled directly. Secondly, the terms inside the parentheses are always less than unity. Thus, the sampling scheme for \(q(x)\) is as follows. We sample a random number \(\xi_1\) on the interval \([0,1)\) and if

then we sample the probability distribution \(2x^3 e^{-x^2}\) for \(x\) using rule C49 in the Monte Carlo Sampler which we can then use to determine the speed of the target nucleus \(v_T\) from equation (73). Otherwise, we sample the probability distribution \(\frac{4}{\sqrt{\pi}} x^2 e^{-x^2}\) for \(x\) using rule C61 in the Monte Carlo Sampler.

With a target speed sampled, we must then decide whether to accept it based on the probability in equation (70). The cosine can be sampled isotropically as \(\mu = 2\xi_2 - 1\) where \(\xi_2\) is a random number on the unit interval. Since the maximum value of \(f_1(v_T, \mu)\) is \(4\sigma_s / \sqrt{\pi} C'\), we then sample another random number \(\xi_3\) and accept the sampled target speed and cosine if

If is not accepted, then we repeat the process and resample a target speed and cosine until a combination is found that satisfies equation (76).

### 5.10.2. Energy-Dependent Cross Section Model¶

As was noted earlier, assuming that the elastic scattering cross section is constant in (56) is not strictly correct, especially when low-lying resonances are present in the cross sections for heavy nuclides. To correctly account for energy dependence of the scattering cross section entails performing another rejection step. The most common method is to sample \(\mu\) and \(v_T\) as in the constant cross section approximation and then perform a rejection on the ratio of the 0 K elastic scattering cross section at the relative velocity to the maximum 0 K elastic scattering cross section over the range of velocities considered:

where it should be noted that the maximum is taken over the range \([v_n - 4/\beta, 4_n + 4\beta]\). This method is known as Doppler broadening rejection correction (DBRC) and was first introduced by Becker et al.. OpenMC has an implementation of DBRC as well as an accelerated sampling method that samples the relative velocity directly.

## 5.11. S(\(\alpha,\beta,T\)) Tables¶

Note that S(\(\alpha,\beta,T\)) tables are only applicable to continuous-energy transport.

For neutrons with thermal energies, generally less than 4 eV, the kinematics of scattering can be affected by chemical binding and crystalline effects of the target molecule. If these effects are not accounted for in a simulation, the reported results may be highly inaccurate. There is no general analytic treatment for the scattering kinematics at low energies, and thus when nuclear data is processed for use in a Monte Carlo code, special tables are created that give cross sections and secondary angle/energy distributions for thermal scattering that account for thermal binding effects. These tables are mainly used for moderating materials such as light or heavy water, graphite, hydrogen in ZrH, beryllium, etc.

The theory behind S(\(\alpha,\beta,T\)) is rooted in quantum mechanics and is quite complex. Those interested in first principles derivations for formulae relating to S(\(\alpha,\beta,T\)) tables should be referred to the excellent books by [Williams] and [Squires]. For our purposes here, we will focus only on the use of already processed data as it appears in the ACE format.

Each S(\(\alpha,\beta,T\)) table can contain the following:

Thermal inelastic scattering cross section;

Thermal elastic scattering cross section;

Correlated energy-angle distributions for thermal inelastic and elastic scattering.

Note that when we refer to “inelastic” and “elastic” scattering now, we are
actually using these terms with respect to the *scattering system*. Thermal
inelastic scattering means that the scattering system is left in an excited
state; no particular nucleus is left in an excited state as would be the case
for inelastic level scattering. In a crystalline material, the excitation of the
scattering could correspond to the production of phonons. In a molecule, it
could correspond to the excitation of rotational or vibrational modes.

Both thermal elastic and thermal inelastic scattering are generally divided into
incoherent and coherent parts. Coherent elastic scattering refers to scattering
in crystalline solids like graphite or beryllium. These cross sections are
characterized by the presence of *Bragg edges* that relate to the crystal
structure of the scattering material. Incoherent elastic scattering refers to
scattering in hydrogenous solids such as polyethylene. As it occurs in ACE data,
thermal inelastic scattering includes both coherent and incoherent effects and
is dominant for most other materials including hydrogen in water.

### 5.11.1. Calculating Integrated Cross Sections¶

The first aspect of using S(\(\alpha,\beta,T\)) tables is calculating cross sections to replace the data that would normally appear on the incident neutron data, which do not account for thermal binding effects. For incoherent inelastic scattering, the cross section is stored as a linearly interpolable function on a specified energy grid. For coherent elastic data, the cross section can be expressed as

where \(E_i\) are the energies of the Bragg edges and \(s_i\) are related to crystallographic structure factors. Since the functional form of the cross section is just 1/E and the proportionality constant changes only at Bragg edges, the proportionality constants are stored and then the cross section can be calculated analytically based on equation (78). For incoherent elastic data, the cross section can be expressed as

where \(\sigma_b\) is the characteristic bound cross section and \(W'\) is the Debye-Waller integral divided by the atomic mass.

### 5.11.2. Outgoing Angle for Coherent Elastic Scattering¶

Another aspect of using S(\(\alpha,\beta,T\)) tables is determining the outgoing energy and angle of the neutron after scattering. For incoherent and coherent elastic scattering, the energy of the neutron does not actually change, but the angle does change. For coherent elastic scattering, the angle will depend on which Bragg edge scattered the neutron. The probability that edge \(i\) will scatter then neutron is given by

After a Bragg edge has been sampled, the cosine of the angle of scattering is given analytically by

where \(E_i\) is the energy of the Bragg edge that scattered the neutron.

### 5.11.3. Outgoing Angle for Incoherent Elastic Scattering¶

For incoherent elastic scattering, OpenMC has two methods for calculating the cosine of the angle of scattering. The first method uses the Debye-Waller integral, \(W'\), and the characteristic bound cross section as given directly in an ENDF-6 formatted file. In this case, the cosine of the angle of scattering can be sampled by inverting equation 7.4 from the ENDF-6 Format:

where \(\xi\) is a random number sampled on unit interval and \(c = 2EW'\). In the second method, the probability distribution for the cosine of the angle of scattering is represented as a series of equally-likely discrete cosines \(\mu_{i,j}\) for each incoming energy \(E_i\) on the thermal elastic energy grid. First the outgoing angle bin \(j\) is sampled. Then, if the incoming energy of the neutron satisfies \(E_i < E < E_{i+1}\) the cosine of the angle of scattering is

where the interpolation factor is defined as

To better represent the true, continuous nature of the cosine distribution, the sampled value of \(mu'\) is then “smeared” based on the neighboring values. First, values of \(\mu\) are calculated for outgoing angle bins \(j-1\) and \(j+1\):

Then, a final cosine is calculated as:

where \(\xi\) is again a random number sampled on the unit interval. Care must be taken to ensure that \(\mu\) does not fall outside the interval \([-1,1]\).

### 5.11.4. Outgoing Energy and Angle for Inelastic Scattering¶

Each S(\(\alpha,\beta,T\)) table provides a correlated angle-energy secondary distribution for neutron thermal inelastic scattering. There are three representations used in the ACE thermal scattering data: equiprobable discrete outgoing energies, non-uniform yet still discrete outgoing energies, and continuous outgoing energies with corresponding probability and cumulative distribution functions provided in tabular format. These three representations all represent the angular distribution in a common format, using a series of discrete equiprobable outgoing cosines.

#### 5.11.4.1. Equi-Probable Outgoing Energies¶

If the thermal data was processed with \(iwt = 1\) in NJOY, then the outgoing energy spectra is represented in the ACE data as a set of discrete and equiprobable outgoing energies. The procedure to determine the outgoing energy and angle is as such. First, the interpolation factor is determined from equation (84). Then, an outgoing energy bin is sampled from a uniform distribution and then interpolated between values corresponding to neighboring incoming energies:

where \(E_{i,j}\) is the j-th outgoing energy corresponding to the i-th incoming energy. For each combination of incoming and outgoing energies, there is a series equiprobable outgoing cosines. An outgoing cosine bin is sampled uniformly and then the final cosine is interpolated on the incoming energy grid:

where \(\mu_{i,j,k}\) is the k-th outgoing cosine corresponding to the j-th outgoing energy and the i-th incoming energy.

#### 5.11.4.2. Skewed Equi-Probable Outgoing Energies¶

If the thermal data was processed with \(iwt=0\) in NJOY, then the outgoing energy spectra is represented in the ACE data according to the following: the first and last outgoing energies have a relative probability of 1, the second and second-to-last energies have a relative probability of 4, and all other energies have a relative probability of 10. The procedure to determine the outgoing energy and angle is similar to the method discussed above, except that the sampled probability distribution is now skewed accordingly.

#### 5.11.4.3. Continuous Outgoing Energies¶

If the thermal data was processed with \(iwt=2\) in NJOY, then the outgoing energy spectra is represented by a continuous outgoing energy spectra in tabular form with linear-linear interpolation. The sampling of the outgoing energy portion of this format is very similar to Correlated Energy and Angle Distribution, but the sampling of the correlated angle is performed as it was in the other two representations discussed in this sub-section. In the Law 61 algorithm, we found an interpolation factor \(f\), statistically sampled an incoming energy bin \(\ell\), and sampled an outgoing energy bin \(j\) based on the tabulated cumulative distribution function. Once the outgoing energy has been determined with equation (34), we then need to decide which angular distribution data to use. Like the linear-linear interpolation case in Law 61, the angular distribution closest to the sampled value of the cumulative distribution function for the outgoing energy is utilized. The actual algorithm utilized to sample the outgoing angle is shown in equation (88). As in the case of incoherent elastic scattering with discrete cosine bins, the sampled cosine is smeared over neighboring angle bins to better approximate a continuous distribution.

## 5.12. Unresolved Resonance Region Probability Tables¶

Note that unresolved resonance treatments are only applicable to continuous-energy transport.

In the unresolved resonance energy range, resonances may be so closely spaced that it is not possible for experimental measurements to resolve all resonances. To properly account for self-shielding in this energy range, OpenMC uses the probability table method. For most thermal reactors, the use of probability tables will not significantly affect problem results. However, for some fast reactors and other problems with an appreciable flux spectrum in the unresolved resonance range, not using probability tables may lead to incorrect results.

Probability tables in the ACE format are generated from the UNRESR module in NJOY following the method of Levitt. A similar method employed for the RACER and MC21 Monte Carlo codes is described in a paper by Sutton and Brown. For the discussion here, we will focus only on use of the probability table table as it appears in the ACE format.

Each probability table for a nuclide contains the following information at a number of incoming energies within the unresolved resonance range:

Cumulative probabilities for cross section bands;

Total cross section (or factor) in each band;

Elastic scattering cross section (or factor) in each band;

Fission cross section (or factor) in each band;

\((n,\gamma)\) cross section (or factor) in each band; and

Neutron heating number (or factor) in each band.

It should be noted that unresolved resonance probability tables affect only integrated cross sections and no extra data need be given for secondary angle/energy distributions. Secondary distributions for elastic and inelastic scattering would be specified whether or not probability tables were present.

The procedure for determining cross sections in the unresolved range using probability tables is as follows. First, the bounding incoming energies are determined, i.e. find \(i\) such that \(E_i < E < E_{i+1}\). We then sample a cross section band \(j\) using the cumulative probabilities for table \(i\). This allows us to then calculate the elastic, fission, and capture cross sections from the probability tables interpolating between neighboring incoming energies. If interpolation is specified, then the cross sections are calculated as

where \(\sigma_{i,j}\) is the j-th band cross section corresponding to the i-th incoming neutron energy and \(f\) is the interpolation factor defined in the same manner as (84). If logarithmic interpolation is specified, the cross sections are calculated as

where the interpolation factor is now defined as

A flag is also present in the probability table that specifies whether an inelastic cross section should be calculated. If so, this is done from a normal reaction cross section (either MT=51 or a special MT). Finally, if the cross sections defined are above are specified to be factors and not true cross sections, they are multiplied by the underlying smooth cross section in the unresolved range to get the actual cross sections. Lastly, the total cross section is calculated as the sum of the elastic, fission, capture, and inelastic cross sections.

## 5.13. Variance Reduction Techniques¶

### 5.13.1. Survival Biasing¶

In problems with highly absorbing materials, a large fraction of neutrons may be
killed through absorption reactions, thus leading to tallies with very few
scoring events. To remedy this situation, an algorithm known as *survival
biasing* or *implicit absorption* (or sometimes *implicit capture*, even though
this is a misnomer) is commonly used.

In survival biasing, absorption reactions are prohibited from occurring and instead, at every collision, the weight of neutron is reduced by probability of absorption occurring, i.e.

where \(w'\) is the weight of the neutron after adjustment and \(w\) is the weight of the neutron before adjustment. A few other things need to be handled differently if survival biasing is turned on. Although fission reactions never actually occur with survival biasing, we still need to create fission sites to serve as source sites for the next generation in the method of successive generations. The algorithm for sampling fission sites is the same as that described in Fission. The only difference is in equation (14). We now need to produce

fission sites, where \(w\) is the weight of the neutron before being
adjusted. One should note this is just the expected number of neutrons produced
*per collision* rather than the expected number of neutrons produced given that
fission has already occurred.

Additionally, since survival biasing can reduce the weight of the neutron to very low values, it is always used in conjunction with a weight cutoff and Russian rouletting. Two user adjustable parameters \(w_c\) and \(w_s\) are given which are the weight below which neutrons should undergo Russian roulette and the weight should they survive Russian roulette. The algorithm for Russian rouletting is as follows. After a collision if \(w < w_c\), then the neutron is killed with probability \(1 - w/w_s\). If it survives, the weight is set equal to \(w_s\). One can confirm that the average weight following Russian roulette is simply \(w\), so the game can be considered “fair”. By default, the cutoff weight in OpenMC is \(w_c = 0.25\) and the survival weight is \(w_s = 1.0\). These parameters vary from one Monte Carlo code to another.

### 5.13.2. Weight Windows¶

In fixed source problems, it can often be difficult to obtain sufficiently low variance on tallies in regions that are far from the source. The weight window method was developed to increase the population of particles in important spatial regions and energy ranges by controlling particle weights. Each spatial region and particle energy range is assigned upper and lower weight bounds, \(w_u\) and \(w_\ell\), respectively. When a particle is in a given spatial region / energy range, its weight, \(w\), is compared to the lower and upper bounds. If the weight of the particle is above the upper weight bound, the particle is split into \(N\) particles, where

and \(N_{max}\) is a user-defined maximum number of splits. To ensure a fair game, each of the \(N\) particles is assigned a weight \(w/N\). If the weight is below \(w_\ell\), it is Russian rouletted as described in Survival Biasing with a survival weight \(w_s\) that is set equal to

where \(f_s\) is a user-defined survival weight ratio greater than one.

On top of the standard weight window method described above, OpenMC implements two additional checks intended to mitigate problems with long histories. First, particles with a weight that falls below some very small cutoff (defaults to \(10^{-38}\)) are killed with no Russian rouletting. Additionally, the total number of splits experienced by a particle is tracked and if it reaches some maximum value, it is prohibited from splitting further.

At present, OpenMC allows weight windows to be defined on all supported mesh types.

References

- Gelbard
Ely M. Gelbard, “Epithermal Scattering in VIM,” FRA-TM-123, Argonne National Laboratory (1979).

- Squires
G. L. Squires,

*Introduction to the Theory of Thermal Neutron Scattering*, Cambridge University Press (1978).- Williams
M. M. R. Williams,

*The Slowing Down and Thermalization of Neutrons*, North-Holland Publishing Co., Amsterdam (1966).**Note:**This book can be obtained for free from the OECD.