Batson Iii committed Nov 17, 2020 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 

CENTRM: A Neutron Transport Code for Computing Continuous-Energy Spectra in General One-Dimensional Geometries and Two-Dimensional Lattice Cells

M. L. Williams

Abstract

CENTRM computes continuous-energy neutron spectra for infinite media, general one-dimensional (1D) systems, and two-dimensional (2D) unit cells in a lattice, by solving the Boltzmann transport equation using a combination of pointwise and multigroup nuclear data. Several calculational options are available, including a slowing-down computation for homogeneous infinite media, 1D discrete ordinates in slab, spherical, or cylindrical geometries; a simplified two-region solution; and 2D method of characteristics for a unit cell within a square-pitch lattice. In SCALE, CENTRM is used mainly to calculate problem-specific fluxes on a fine energy mesh (10,000–70,000 points), which may be used to generate self-shielded multigroup cross sections for subsequent radiation transport computations.

ACKNOWLEDGMENTS

Several current and former ORNL staff made valuable contributions to the CENTRM development. The author acknowledges the contributions of former ORNL staff members D. F. Hollenbach and N. M. Greene; as well as current staff member L. M. Petrie. Special thanks go to Kang-Seog Kim who developed the 2D method of characteristics option for CENTRM. Portions of the original code development were performed by M. Asgari as partial fulfillment of his PhD dissertation research at Louisiana State University (LSU); and Riyanto Raharjo from LSU made significant programming contributions for the inelastic scattering and thermal calculations.

Introduction

CENTRM (Continuous ENergy TRansport Module) computes “continuous-energy” neutron spectra using various deterministic approximations to the Boltzmann transport equation. Computational methods are available for infinite media, general one-dimensional (1-D) geometries, and two-dimensional (2D) unit cells in a square-pitch lattice. The purpose of the code is to provide fluxes and flux moments for applications that require a high resolution of the fine-structure variation in the neutron energy spectrum. The major function of CENTRM is to determine problem-specific fluxes for processing multigroup (MG) data with the XSProc self-shielding module (Introduction in XSProc chapter), which is executed by all SCALE MG sequences. XSProc calls an application program interface (API) to perform a CENTRM calculation for a representative model (e.g., a unit cell in a lattice), and then utilizes the spectrum as a problem-dependent weight function for MG averaging. The MG data processing is done in XSProc by calling an API for the PMC code, which uses the CENTRM continuous-energy (CE) flux spectra and cross-section data to calculate group-averaged cross sections over some specified energy range. The resulting application-specific library is used for MG neutron transport calculations within SCALE sequences. In this approach the CENTRM/PMC cross-section processing in XSProc becomes an active component in the overall transport analysis. CENTRM can also be executed as a standalone code, if the user provides all required input data and nuclear data libraries; but execution through XSProc is much simpler and less prone to error.

Description of problem solved

CENTRM uses a combination of MG and pointwise (PW) solution techniques to solve the neutron transport equation over the energy range ~0 to 20 MeV. The calculated CE spectrum consists of PW values for the flux per unit lethargy defined on a discrete energy mesh, for which a linear variation of the flux between energy points is assumed. Depending on the specified transport approximation, the flux spectrum may vary as a function of space and direction, in addition to energy. Spherical harmonic moments of the angular flux, which may be useful in processing MG matrices for higher order moments of the scattering cross section, can also be determined as a function of space and energy mesh.

CENTRM solves the fixed-source (inhomogeneous) form of the transport equation, with a user-specified fixed source term. The input source may correspond to MG histogram spectra for volumetric or surface sources or it may be a “fission source” which has a continuous-energy fission-spectrum distribution (computed internally) appropriate for each fissionable mixture. Note that eigenvalue calculations are not performed in CENTRM—these must be performed by downstream MG transport codes that utilize the self-shielded data processed with the CENTRM spectra.

Nuclear data required for CENTRM

A MG cross section library, a CE cross section library, and a CE thermal kernel [S(α, β)] library are required for the CENTRM transport calculation. During XSProc execution for a given unit cell in the CELL DATA block, the MG library specified in the input is processed by BONAMI prior to the CENTRM calculation, in order to provide self-shielded data based on the Bondarenko approximation for the MG component of the CENTRM solution. The shielded MG cross sections are also used in CENTRM to correct infinitely dilute CE data in the unresolved resoance range. The CRAWDAD module is executed by XSProc to generate the CENTRM CE cross section and thermal kernel libraries, respectively, by concatenating discrete PW data read from individual files for the nuclides in the unit cell mixtures. CE resonance profiles are based strictly on specifications in the nuclear data evaluations; e.g., Reich-Moore formalism is specified for most materials in ENDF/B-VII. PW data in the CENTRM library are processed such that values at any energy can be obtained by linear interpolation within some error tolerance specified during the library generation (usually ~0.1% or less). CRAWDAD also interpolates the CE cross section data and the Legendre moments of the thermal scattering kernels to the appropriate temperatures for the unit cell mixtures. The format of the CENTRM library is described in Description of the CENTRM CE cross section file.

Code assumptions and features

 Batson Iii committed Jan 19, 2021 141 

As shown in Fig. 206, the energy range of interest is divided into  Batson Iii committed Nov 17, 2020 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 three intervals called the Upper Multigroup Range (UMR), Pointwise Range (PW), and Lower Multigroup Range (LMR), respectively, which are defined by input. MG fluxes are computed using standard multigroup techniques for the UMR and LMR, and these values are then divided by the group lethargy width to obtain the average flux per lethargy within each group. This “pseudo-pointwise” flux is assigned to the midpoint lethargy of the group, so that there is one energy point per group in the UMR and LMR energy intervals. However, for each group in the PW range there are generally several, and possibly many, energy points for which CENTRM computes flux values. In this manner a problem-dependent spectrum is obtained over the entire energy range.

The default PW range goes from 0.001 eV to 20 keV, but the user can modify the PW limits. The energy range for the PW transport calculation is usually chosen to include the interval where the important absorber nuclides have resolved resonances, while MG calculations are performed where the cross sections characteristically have a smoother variation or where shielding effects are less important. In the SCALE libraries the thermal range is defined to be energies less than 5.0 eV. Above thermal energies, scattering kinematics are based on the stationary nucleus model, while molecular motion and possible chemical binding effects are taken into account for thermal scattering, which can result in an incease in the neutron incident energy. The CENTRM thermal calculation uses Legendre coefficients from the CE kernel library that describes point-to-point energy transfers for incoherent and coherent scattering, as function of temperature, for all moderators that have thermal scattering law data provided in ENDF/B. Thermal kernels for all other materials are generated internally by CENTRM based on the free-gas model.

 Batson Iii committed Jan 19, 2021 172 

Fig. 206 Definition of UMR, PW, and LMR energy ranges.

 Batson Iii committed Nov 17, 2020 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 

Several transport computation methods are available for both MG and PW calculations. These include a space-independent slowing down calculation for infinite homogeneous media, 1D discrete ordinates or P1 methods for slab, spherical, and cylindrical geometries, and a 2D method of chracateristics (MoC) method for lattice unit cells. A simplified two-region collision-probablity method is also available for ther pointwise solution. In general the user may specify different transport methods for the UMR, PW, and LMR, respectively; however, if the 2D MoC method is specified for any range, it will be used for all.

The CENTRM 1D discrete ordinates calculation option has many of the same features as the SCALE MG code XSDRNPM. It represents the directional dependence of the angular flux with an arbitrary symmetric-quadrature order, and uses Legendre expansions up to P5 to represent the scattering source. No restrictions are placed on the material arrangement or the number of spatial intervals in the calculation, and general boundary conditions (vacuum, reflected, periodic, albedo) can be applied on either boundary of the 1D geometry. Lattice cells are represented in the CENTRM discrete ordinates option by a 1D Wigner-Sitz cylindrical or spherical model with a white boundary condition on the outer surface.

Starting with SCALE-6.2, CENTRM also includes a 2D MoC solver for lattice cell geometries consisting of a cylindrical fuel rod (fuel/gap/clad) contained within a rectangular moderator region. The MoC calculation is presently limited to square lattices. The 2D unit cell uses a reflected boundary condition on the outer square surface, which provides a more rigorous treatment than the 1D Wigner-Seitz model; however the MoC option requires a longer execution time than the 1D discrete ordinates method. The MoC option has been found to improve results compared to the 1D Wigner-Seitz cell model for many cases, but in other cases the improvement is marginal.

A variable PW energy mesh is generated internally to accurately represent the fine-structure flux spectrum for the system of interest. This gives CENTRM the capability to rigorously account for resonance interference effects in systems with multiple resonance absorbers. Because CENTRM calculates the space-dependent PW flux spectrum, the spatial variation of the self-shielded cross sections within an absorber body can be obtained. A radial temperature distribution can also be specified, so that space-dependent Doppler broadening can be treated in the transport solution. Within the epithermal PW range, the slowing-down source due to elastic and discrete-level inelastic reactions is computed with the analytical scatter kernel based upon the neutron kinematic relations for s-wave scattering. Continuum inelastic scatter is approximated by an analytical evaporation spectrum, assumed isotropic in the laboratory system. For many thermal reactor and criticality safety problems, self-shielding of inelastic cross sections has a minor impact, and by default these options are turned off for faster execution. As previously discussed, the thermal scatter kernel is based on the ENDF/B scattering law data for bound moderators, and uses the free-gas model for other materials.

Theory and Analytical Models

This section describes the coupled MG and PW techniques used to solve the neutron transport equation.

Energy/lethargy ranges for MG and PW calculations

The combined MG/PW CENTRM calculation is performed over the energy range spanned by the group structure in the input MG library. The energy boundaries for the “IGM” neutron groups specified on the MG library divide the entire energy range into energy intervals. The lowest energy group contained in the UMR is defined to be “MGHI”; while the highest energy group in the LMR is designated “MGLO.” The boundary between the PW and UMR energy intervals is set by the energy value “DEMAX,” while “DEMIN” is the boundary between the PW and LMR. The default values of 0.001 eV and 20 keV for DEMIN and DEMAX, respectively, can be modified by user input, but the input values are altered by the code to correspond to the closest group boundaries. Hence, DEMAX is always equal to the lower energy boundary of group MGHI and DEMIN the upper energy boundary of MGLO. The PW calculation is performed in terms of lethargy (u), rather than energy (E). The origin (u=0) of the lethargy coordinate corresponds to the energy E=DEMAX, which is the top of the PW range. See  Batson Iii committed Jan 19, 2021 246 Fig. 207.

 Batson Iii committed Nov 17, 2020 247 248 249 250 251 252 253 254 255 256 

The highest energy group of the thermal range is defined by the parameter “IFTG,” obtained from the MG library. If DEMIN is less than the upper energy boundary of IFTG, the PW range extends into thermal. In this case, scattering in the PW region of the thermal range is based on the PW scattering kernel data; and the LMR calculation uses 2D transfer matrices for incoherent and coherent scattering on the MG library. Coupling between the MG and PW thermal calculations is treated, and outer iterations are required to address effects of upscattering.

 Batson Iii committed Jan 19, 2021 257 

Fig. 207 Definition of High and Transition regions in upper multigroup region.

 Batson Iii committed Nov 17, 2020 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 

With the exception of hydrogen moderation, elastic down-scattering coupling the UMR and PW ranges, occurs only within a limited sub-range of the UMR called the “transition region”. The highest energy group in the transition region is designated “MGTOP.” The precise definition of the transition region is given in Scattering sources for the PW range.

Energy boundaries of the group structure on the input MG library correspond to the IGM+1 values, { G1, G2, … Gg, Gg+1, …, GIGM+1}. It is convenient to designate the number of groups in the UMR, PW, and LMR ranges equal to NGU, NGP, and NGL, respectively, so that IGM = NGU + NGP + NGL; or in terms of the parameters MGHI and MGLO introduced previously:

NGU = MGHI; NGP = MGLO − MGHI − 1; NGL = IGM − MGLO + 1.

The flux per unit lethargy is calculated for a discrete energy (or lethargy) mesh spanning the MG structure. Groups in the UMR and LMR each contain a single energy mesh point, while groups in the PW range generally contain several points. The number of mesh points in the UMR, PW, and LMR is equal respectively to NGU, NP, and NGL; and the total number of points in the entire energy mesh is designated as “NT,” which is equal to NGU + NP + NGL. Thus the lethargy (u) mesh consists of the set of points: {u:sub:1,….uNGU, uNGU+1,….uNGU+NP, uNGU+NP+1,…uNT}. Based on the lethargy origin at E=DEMAX, the lethargy “un” associated with any energy point “En” is equal to,

un = ln(DEMAX/En).

Lethargy points are arranged in order of increasing value. The lethargy origin is at point NGU+1, the lower energy boundary of group MGHI; i.e., uNGU+1=0. Note that the entire UMR (E>DEMAX) corresponds to negative lethargy values. Lethargy values for the first NGU and the last NGL points in the mesh are defined to be the midpoint lethargies of groups in the UMR and LMR ranges, respectively. For example, for the NGU groups within the UMR,

u1 = 0.5[ln(DEMAX/G1) + ln(DEMAX/G2)];

uNGU = 0.5[ln(DEMAX/GMGHI) + ln(DEMAX/GMGHI + 1)];

and similarly for the NGL groups in the LMR,

uNGU + NP + 1 = 0.5[ln(DEMAX/GMGLO) + ln(DEMAX/GMGLO + 1)]

uNT = 0.5[ln(DEMAX/GIGM) + ln(DEMAX/GIGM + 1)]

The remaining NP points in the mesh (i.e., values uNGU+1 to uNGU+NP) are contained within the NGP groups that span the PW range. By definition the first point in the PW range is the lower energy boundary of group MGHI. The other mesh points are computed internally by CENTRM, based on the behavior of the macroscopic PW total cross sections and other criteria.

The neutron flux, as a function of space and direction, is calculated for each energy/lethargy point in the mesh by solving the Boltzmann transport equation. The transport equation at each lethargy point generally includes a source term representing the production rate due to elastic and inelastic scatter from other lethargies, which couples the solutions at different lethargy mesh points. Except in the thermal range, neutrons can only gain lethargy (lose energy) in a scattering reaction; thus the PW flux is computed by solving the transport equation at successive mesh points, sweeping from low to high lethargy values.

The Boltzmann equation for neutron transport

The steady state neutron transport equation shown below represents a particle balance-per unit phase space, at an arbitrary point ρ in phase space,

 Batson Iii committed Jan 19, 2021 330 (326)$\Omega \cdot \nabla \Psi(\rho)+\sum_{t}(\mathrm{r}, \mathrm{u}) \Psi(\rho)=\int_{0}^{\infty} \int_{0}^{4 \pi} \Sigma\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u} ; \mu_{0}\right) \Psi\left(\mathrm{u}^{\prime}, \Omega^{\prime}\right) \mathrm{d} \mathrm{u}^{\prime}+\mathrm{Q}_{\mathrm{ext}}(\rho)$
 Batson Iii committed Nov 17, 2020 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 

where:

ψ(p) = angular flux (per lethargy) at phase space coordinate ρ;

ρ = (r,u,Ω) = phase space point defined by the six independent variables;

r = (x1,x2,x3) = space coordinates;

u = ln(Eref/E) = lethargy at energy E, relative to an origin (u=0) at Eref;

Ω = (μ,ζ) = neutron direction defined by polar cosine μ and azimuthal angle ζ;

Σt(r,u) = macroscopic total cross section;

Σ(u′→u;μ0) = double differential scatter cross section;

μ0 = cosine of scatter angle, measured in laboratory coordinate system;

Qext(ρ) = external source term, including fission source;

 Batson Iii committed Jan 19, 2021 347 

The left and right sides of (326) respectively, are equal to the neutron  Batson Iii committed Nov 17, 2020 348 349 350 351 352 353 354 355 loss and production rates, per unit volume-direction-lethargy. In CENTRM the spatial distribution of the fission source is input as a component of the external source Q; hence, a fixed source rather than an eigenvalue calculation is required for the transport solution.

The angular dependence of the double-differential macroscopic scatter cross section of an arbitrary nuclide “j” is represented by a finite Legendre expansion of arbitrary order L:

 Batson Iii committed Jan 19, 2021 356 (327)$\Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u} ; \mu_{0}\right)=\sum_{=0}^{\mathrm{L}} \frac{2+1}{2} \Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}\right) \quad \mathrm{P}\left(\mu_{0}\right)$
 Batson Iii committed Nov 17, 2020 357 358 359 360 361 362 

where P(μ:sub:0) = Legendre polynomial evaluated at the laboratory scattering cosine μ0; and

$$\Sigma^{(j)}\left(u^{\prime} \rightarrow u\right)$$ = cross section moments of nuclide j, defined by the expression

 Batson Iii committed Jan 19, 2021 363 (328)$\Sigma^{(j)}\left(u^{\prime} \rightarrow u\right)=\int_{-1}^{1} \Sigma^{(j)}\left(u^{\prime} \rightarrow u ; \mu_{0}\right) P\left(\mu_{0}\right) d \mu_{0}$
 Batson Iii committed Nov 17, 2020 364 365 366 

After substitution of the above Legendre expansions for the scattering data of each nuclide, and applying the spherical harmonic addition theorem in the usual manner, the scattering source on the right side of  Batson Iii committed Jan 19, 2021 367 (326) becomes [BG70]:

 Batson Iii committed Nov 17, 2020 368 
 Batson Iii committed Jan 19, 2021 369 (329)$\mathrm{S}(\mathrm{r}, \mathrm{u}, \Omega) \equiv \int_{0}^{\infty} \int_{0}^{4 \pi} \Sigma\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u} ; \mu_{0}\right) \Psi\left(\mathrm{u}^{\prime}, \Omega^{\prime}\right) \mathrm{d} \Omega^{\prime} \mathrm{d} \mathrm{u}^{\prime}=\sum_{\mathrm{k}=1}^{\mathrm{LK}} \frac{2+1}{2} \mathrm{Y}_{\mathrm{k}}(\Omega) \mathrm{S}_{\mathrm{k}}(\mathrm{r}, \mathrm{u})$
 Batson Iii committed Nov 17, 2020 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 

wherein,

$$\mathrm{Y}_{\mathrm{k}}(\Omega)=\mathrm{Y}_{\mathrm{k}}(\mu, \zeta)$$ = the spherical harmonic function evaluated at direction Ω

Sk = spherical harmonic moments of the scatter source, per unit letharagy.

The summation index “ℓk” indicates a double sum over ℓ and k indices; in the most general case it is defined as:

$\sum_{\mathrm{k}=1}^{\mathrm{LK}}=\sum_{=0}^{\mathrm{L}} \sum_{\mathrm{k}=0}$

where “L” is the input value for the maximum order of scatter (input parameter “ISCT”).

Due to symmetry conditions, some of the source moments may be zero. The parameter LK is defined to be the total number of non-zero moments (including scalar flux) for the particular geometry of interest, and is equal to,

LK = L + 1 for 1D slabs and spheres;

LK = L*(L+4)/4+1 for 1D cylinders, and

LK = L*(L+3)/2+1 for 2D MoC cells

More details concerning the 1-D Boltzmann equation can be found in the XSDRNPM chapter of the SCALE manual.

 Batson Iii committed Jan 19, 2021 392 

The Sℓk moments in (329)  correspond to expansion coefficients in  Batson Iii committed Nov 17, 2020 393 394 395 a spherical harmonic expansion of the scatter source. These can be expressed in terms of the cross section and flux moments by

 Batson Iii committed Jan 19, 2021 396 (330)$\mathrm{S}_{\mathrm{k}}(\mathrm{u})=\sum_{\mathrm{j}} \int_{\mathrm{u}^{\prime}} \mathrm{S}_{\mathrm{k}}^{(\mathrm{j})}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}\right) \mathrm{d} \mathrm{u}^{\prime}=\sum_{\mathrm{j}} \int_{\mathrm{u}^{\prime}} \Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}\right) \Psi_{\mathrm{k}}\left(\mathrm{u}^{\prime}\right) \mathrm{d} \mathrm{u}^{\prime}$
 Batson Iii committed Nov 17, 2020 397 398 399 

where ψℓk(u) = spherical harmonic moments of the angular flux;

 Batson Iii committed Jan 19, 2021 400 (331)$= \int_{0}^{4 \pi} \mathrm{Y}_{\mathrm{k}}(\Omega) \Psi(\Omega) \mathrm{d} \Omega$
 Batson Iii committed Nov 17, 2020 401 402 403 

and Sℓk(j)(u′→u) = moments of the differential scatter rate from lethargy u′ to u, for nuclide “j”;

 Batson Iii committed Jan 19, 2021 404 (332)$= \Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}\right) \Psi_{\mathrm{k}}\left(\mathrm{u}^{\prime}\right)$
 Batson Iii committed Nov 17, 2020 405 406 407 408 409 410 

The ψℓk flux moments are the well known coefficients appearing in a spherical harmonic expansion of the angular flux. These usually are the desired output from the transport calculation. In particular, the ℓ=0, k=0 moment corresponds to the scalar flux [indicated here as Φ(r,u)],

 Batson Iii committed Jan 19, 2021 411 412 (333)$\Psi_{0,0}(\mathrm{r}, \mathrm{u})=\Phi(\mathrm{r}, \mathrm{u})=\int_{0}^{4 \pi} \Psi(\mathrm{r}, \mathrm{u}, \Omega) d \Omega$

In general the epithermal component of the scatter source in (329)  Batson Iii committed Nov 17, 2020 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 contains contributions from both elastic and inelastic scatter reactions; however, inelastic scatter is only possible above the threshold energy corresponding to the lowest inelastic level. The inelastic Q values for most materials are typically above 40 keV; therefore, elastic scatter is most important for slowing down calculations in the resolved resonance range of most absorber materials of interest. For example, the inelastic Q values of 238U, iron, and oxygen are approximately 45 keV, 846 keV, and 6 MeV, respectively; while the upper energy of the 238U resolved resonance range is 20 keV in ENDF/B-VII. The inelastic thresholds of some fissile materials like 235U and 239Pu are on the order of 10 keV; however, with the exception of highly enriched fast systems, these inelastic reactions usually contribute a negligible amount to the overall scattering source. CENTRM assumes that continuum inelastic scatter is isotropic in the laboratory system, while discrete level inelastic scatter is isotropic in the center of mass (CM) coordinate system.

Over a broad energy range, elastic scatter from most moderators can usually be assumed isotropic (s-wave) in the neutron-nucleus CM coordinate system. In the case of hydrogen, this is true up to approximately 13 MeV; for carbon up to 2 MeV; and for oxygen up to 100 keV. However, it is well known that isotropic CM scatter does not result in isotropic scattering in the laboratory system. For s‑wave elastic scatter the average scatter-cosine in the laboratory system is given by: $$\bar{\mu}_{0}=0.667 / \mathrm{~A};^{3}$$ where “A” is the mass number (in neutron mass units) of the scattering material. This relation indicates that s-wave, elastic scattering from low A materials tends to be more anisotropic in the laboratory, and that the laboratory scattering distribution approaches isotropic $$\left(\bar{\mu}_{0}=0 ; \theta_{0}=90\right)$$ as A becomes large. For example, the $$\bar{\mu}_{0}$$ of hydrogen is 0.667 (48.2°); while it is about 0.042 (87.6°) for oxygen. Because s‑wave scattering from heavy materials is nearly isotropic in the laboratory system, the differential scattering cross section (and thus the scattering source) can usually be expressed accurately by a low order Legendre expansion. On the other hand light moderators like hydrogen may require more terms—depending on the flux anisotropy—to accurately represent the elastic scatter source in the laboratory system. The default settings in CENTRM are to use P0 (isotropic lab scatter) for mass numbers greater than A=100, and P1 for lighter masses.

An analytical expression can be derived for the cross-section moments in the case of two-body reactions, such as elastic and discrete-level inelastic scattering from “stationary” nuclei. Stationary here implies that the effect of nuclear motion on neutron scattering kinematics is neglected.

Note

The stationary nucleus approximation for treating scattering kinematics does not imply that the effect of nuclear motion on Doppler broadening of resonance cross sections is ignored, since this effect is included in the PW cross-section data.

In CENTRM the stationary nucleus approximation is applied above the thermal cutoff, typically around 3-5 eV, but is not valid for low energy neutrons. CENTRM has the capability to perform a PW transport calculation in the thermal energy range using tabulated thermal scattering law data for bound molecules, combined with the analytical free-gas kernel for other materials. In this case the cross-section  Batson Iii committed Jan 19, 2021 472 moments appearing in (328)  include upscattering effects. The expressions  Batson Iii committed Nov 17, 2020 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 used in CENTRM to compute the PW scatter source moments in the thermal range are given in The Boltzmann equation within the PW range.

The following two sections discuss the evaluation of the scatter source moments for epithermal elastic and inelastic reactions, respectively.

Epithermal Elastic Scatter

Consider a neutron with energy E′, traveling in a direction Ω′, that scatters elastically from an arbitrary material “j,” having a mass A(j) in neutron-mass units. Conservation of kinetic energy and momentum requires that there be a unique relation between the angle that the neutron scatters (relative to the initial direction) and its final energy E after the collision. If the nucleus is assumed to be stationary in the laboratory coordinate system, then the cosine (μ0) of the scatter angle (θ0) measured in the laboratory system, as a function of the initial and final energies, is found to be

 Batson Iii committed Jan 19, 2021 490 (334)$\mu_{0} \equiv \Omega^{\prime} \cdot \Omega=\mathrm{G}^{(\mathrm{j})}\left(\mathrm{E}^{\prime}, \mathrm{E}\right) ,$
 Batson Iii committed Nov 17, 2020 491 492 493 

where the kinematic correlation function G relating E′, E, and μ0 for elastic scatter is equal to

 Batson Iii committed Jan 19, 2021 494 (335)$\begin{split} \begin{array}{l}  Batson Iii committed Nov 17, 2020 495 496 497 498 499 500 \mathrm{G}^{(\mathrm{j})}\left(\mathrm{E}^{\prime}, \mathrm{E}\right)=\mathrm{a}_{1}^{(\mathrm{j})} \times\left[\mathrm{E} / \mathrm{E}^{\prime}\right]^{\frac{1}{2}}-\mathrm{a}_{2}^{(\mathrm{j})} \times\left[\mathrm{E}^{\prime} / \mathrm{E}\right]^{\frac{1}{2}} \\ \text { and } \mathrm{a}_{1}^{(\mathrm{j})}=\left(\mathrm{A}^{(\mathrm{j})}+1\right) / 2 \quad ; \quad \mathrm{a}_{2}^{(\mathrm{j})}=\left(\mathrm{A}^{(\mathrm{j})}-1\right) / 2 \end{array} .\end{split}$

The final energy E of an elastically scattered neutron is restricted to the range,

 Batson Iii committed Jan 19, 2021 501 (336)$\alpha^{(j)} E^{\prime} \leq E \leq E^{\prime}$
 Batson Iii committed Nov 17, 2020 502 503 

where α(j) =  maximum fractional energy lost by elastic scatter

 Batson Iii committed Jan 19, 2021 504 (337)$= \left[\mathrm{a}_{2}^{(\mathrm{j})} / \mathrm{a}_{1}^{(\mathrm{j})}\right]^{2}$
 Batson Iii committed Nov 17, 2020 505 506 

The corresponding range of scattered neutrons in terms of lethargy is equal to

 Batson Iii committed Jan 19, 2021 507 (338)$\mathrm{u}^{\prime} \leq \mathrm{u} \leq \mathrm{u}^{\prime}+\varepsilon^{(\mathrm{j})}$
 Batson Iii committed Nov 17, 2020 508 509 

where

 Batson Iii committed Jan 19, 2021 510 (339)\begin{split} \begin{aligned}  Batson Iii committed Nov 17, 2020 511 512 513 514 515 516 517 \mathrm{u}, \mathrm{u}^{\prime} &=\mathrm{u}(\mathrm{E}), \mathrm{u}^{\prime}\left(\mathrm{E}^{\prime}\right)=\text { lethargies corresponding to } \mathrm{E} \text { and } \mathrm{E}^{\prime}, \text { respectively; and } \\ \varepsilon^{(\mathrm{j})} &=\text { maximum increase in lethargy, per elastic scatter }=\ln \left[1 / \alpha^{(j)}\right] \end{aligned}\end{split}

The double-differential scatter kernel of nuclide j (per unit lethargy and solid angle) for s-wave elastic scatter of neutrons from stationary nuclei, is equal to

 Batson Iii committed Jan 19, 2021 518 (340)\begin{split} \begin{aligned}  Batson Iii committed Nov 17, 2020 519 520 521 522 523 \Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u} ; \mu_{0}\right) &=\frac{\mathrm{E} \Sigma^{(\mathrm{i})}\left(\mathrm{u}^{\prime}\right)}{\mathrm{E}^{\prime}\left(1-\alpha^{(\mathrm{j})}\right)} \delta\left[\mu_{0}-\mathrm{G}^{(\mathrm{j})}\left(\mathrm{E}^{\prime}, \mathrm{E}\right)\right], \text { for } \mathrm{u}^{\prime} \leq \mathrm{u} \leq \mathrm{u}^{\prime}+\varepsilon^{(\mathrm{j})} \\ &=0 \quad \mathrm{u}<\mathrm{u}^{\prime} \text { or } \mathrm{u}>\mathrm{u}^{\prime}+\varepsilon^{(\mathrm{j})} \end{aligned}\end{split}

The presence of the Dirac delta function completely correlates the angle of scatter and the values of the initial and final energies.  Batson Iii committed Jan 19, 2021 524 525 Substituting the double differential cross-section expression from (340) into (328)  gives the single-differential Legendre moments of the  Batson Iii committed Nov 17, 2020 526 527 cross section, per final lethargy:

 Batson Iii committed Jan 19, 2021 528 (341)\begin{split} \begin{aligned}  Batson Iii committed Nov 17, 2020 529 530 531 532 533 \Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}\right) &=\frac{\mathrm{E} \mathrm{P}\left[\mathrm{G}^{(\mathrm{j})}\right] \Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime}\right)}{\mathrm{E}^{\prime}\left(1-\alpha^{(\mathrm{j})}\right)}, \text { for } \mathrm{u}^{\prime} \leq \mathrm{u} \leq \mathrm{u}^{\prime}+\varepsilon^{(\mathrm{j})} \\ &=0 \quad \mathrm{u}^{\prime} \text { or } \mathrm{u}>\mathrm{u}^{\prime}+\varepsilon^{(\mathrm{j})} \end{aligned}\end{split}

where P = Legendre polynomial evaluated at argument G(j) equal to the scatter cosine.

 Batson Iii committed Jan 19, 2021 534 

When the above expressions are used in (330) , the following is obtained  Batson Iii committed Nov 17, 2020 535 536 537 for the ℓk moment of the epithermal elastic scattering source at lethargy u:

 Batson Iii committed Jan 19, 2021 538 (342)$\mathrm{S}_{\mathrm{k}}^{(\mathrm{j})}(\mathrm{u})=\sum_{\mathrm{j}} \int_{\left.\mathrm{u}-\varepsilon^{(\mathrm{i}}\right)}^{\mathrm{u}} \mathrm{S}_{\mathrm{k}}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}\right) \mathrm{du}^{\prime}=\sum_{\mathrm{j}} \int_{\mathrm{u}-\varepsilon^{(j)}}^{\mathrm{u}} \frac{\mathrm{E} \Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime}\right) \mathrm{P}\left[\mathrm{G}^{(\mathrm{j})}\right]}{\mathrm{E}^{\prime}\left(1-\alpha^{(\mathrm{j})}\right)} \Psi_{\mathrm{k}}\left(\mathrm{u}^{\prime}\right) \mathrm{du}^{\prime} .$
 Batson Iii committed Nov 17, 2020 539 540 541 542 543 544 545 546 547 548 549 550 551 552 

Epithermal Inelastic Scatter

If the input value of DEMAX is set above the inelastic threshold of some materials in the problem, then inelastic scattering can occur in the PW range. The pointwise transport calculation may optionally include discrete-level and continuum inelastic reactions in computing the PW scatter source moments. The multigroup calculations always consider inelastic reactions.

Discrete-level inelastic reactions are two-body interactions, so that kinematic relations can be derived relating the initial and final energies and the angle of scatter. It can be shown that the kinematic correlation function for discrete-level inelastic scatter can be written in a form identical to that for elastic scatter by redefining the  Batson Iii committed Jan 19, 2021 553 parameter a1 in (335)  to be the energy dependent function [Wil00],

 Batson Iii committed Nov 17, 2020 554 
 Batson Iii committed Jan 19, 2021 555 (343)$\mathrm{a}_{1}^{(\mathrm{~m}, \mathrm{j})} = \frac{\left(\mathrm{A}^{(\mathrm{j})}+1\right)}{2}+\frac{\left(-\mathrm{Q}^{(\mathrm{m}, \mathrm{j})}\right) \mathrm{A}^{(\mathrm{j})}}{2 \mathrm{E}}$
 Batson Iii committed Nov 17, 2020 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 

The parameter Q(m, j) is the Q-value for the mth level of nuclide “j”. The Q value is negative for inelastic scattering, while it is zero for elastic scatter. The threshold energy in the laboratory coordinate system is proportional to the Q-value of the inelastic level, and is given by:

$\mathrm{E}_{\mathrm{T}}^{(\mathrm{m}, \mathrm{j})}=\frac{\mathrm{A}^{(\mathrm{j})}+1}{\mathrm{~A}^{(\mathrm{j})}} \times\left(-\mathrm{Q}^{(\mathrm{m}, \mathrm{j})}\right)$

The range of energies that can contribute to the scatter source at E, due to inelastic scatter from the mth level of nuclide j is defined to be [E:sub:L , E:sub:H ], where EH >E:sub:L >E:sub:T . This energy range has a corresponding lethargy range of [u:sub:LO , u:sub:HI ] which is equal to,

$\begin{split}\begin{array}{l} \mathrm{u}_{\mathrm{LO}}=\mathrm{u}-\ln \left(\frac{1}{\alpha_{1}^{(\mathrm{j})}\left(\mathrm{E}_{\mathrm{H}}\right)}\right)=\mathrm{u}-\varepsilon_{1}^{(\mathrm{j})} \\ \mathrm{u}_{\mathrm{HI}}=\mathrm{u}-\ln \left(\frac{1}{\alpha_{2}^{(\mathrm{j})}\left(\mathrm{E}_{\mathrm{L}}\right)}\right)=\mathrm{u}-\varepsilon_{2}^{(\mathrm{j})} \end{array}\end{split}$

The energy-dependent alpha parameters in the above expressions are defined as,

$\begin{split}\begin{array}{l} \alpha_{1}(\mathrm{E})=\left(\frac{\mathrm{A} \Delta^{(\mathrm{m}, \mathrm{j})}(\mathrm{E})-1}{\mathrm{~A}+1}\right)^{2} \\ \alpha_{2}(\mathrm{E})=\left(\frac{\mathrm{A} \Delta^{(\mathrm{m}, \mathrm{j})}(\mathrm{E})+1}{\mathrm{~A}+1}\right)^{2} \end{array}\end{split}$

where

 Batson Iii committed Jan 19, 2021 582 583 (344)$\Delta^{(\mathrm{m}, \mathrm{j})}(\mathrm{E})=\sqrt{1-\frac{\mathrm{E}_{\mathrm{T}}^{(\mathrm{m}, \mathrm{j})}}{\mathrm{E}}}$

Modifying the epithermal elastic scatter source in (342)  to include  Batson Iii committed Nov 17, 2020 584 585 discrete-level inelastic scatter gives the following expression

 Batson Iii committed Jan 19, 2021 586 (345)$\mathrm{S}_{\mathrm{k}}(\mathrm{u})=\sum_{m, j} \int_{u_{L O}}^{u_{H I}} \frac{\mathrm{E}}{\mathrm{E}^{\prime}} \frac{\Sigma^{(\mathrm{m}, \mathrm{j})}\left(\mathrm{E}^{\prime}\right) \mathrm{P}\left[\mathrm{G}^{(\mathrm{m}, \mathrm{j})}\right]}{\left(1-\alpha^{j}\right) \Delta^{(\mathrm{m}, \mathrm{j})}\left(\mathrm{E}^{\prime}\right)} \Psi_{\mathrm{k}}\left(\mathrm{u}^{\prime}\right) \mathrm{d} u^{\prime}$
 Batson Iii committed Nov 17, 2020 587 588 

Detailed expressions for the lethargy limits are given in [Wil00]. Since Δ(m,j) is equal to unity for elastic scatter, the above  Batson Iii committed Jan 19, 2021 589 equation reduces to (339) if there is no discrete-level inelastic  Batson Iii committed Nov 17, 2020 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 contribution.

At high energies, the inelastic levels of the nucleus become a continuum. In this case CENTRM represents the energy distribution of the scattered neutrons by an evaporation spectrum with an isotropic angular distribution in the lab system; thus, only the P0 moment appears in the continuum inelastic scattering source. Including continuum inelastic reactions in the PW calculation usually has a small impact on the spectrum used for resonance self-shielding, and may adversely impact the computer memory requirements and execution time. Therefore, by default, CENTRM does not include continuum inelastic reactions in the pointwise solution; however, it is always included in the UMR solution.

Thermal Scatter

Since thermal neutrons have energies comparable to the mean kinetic energy of molecules in thermal equilibrium, the scattering kernels must account for molecular motion. The scatter moments include both downscatter as well as upscatter contributions; hence, the integration  Batson Iii committed Jan 19, 2021 609 limits appearing in (342)  must be extended from the lowest to the highest  Batson Iii committed Nov 17, 2020 610 611 612 613 614 615 616 617 618 619 energy of the thermal range. Furthermore the cross-section moments correspond to the Legendre expansion coefficients of the thermal scatter kernel, which has a substantially different form than the epithermal kernel discussed in the previous two sections. In general the ℓth Legendre moment of the thermal scattering kernel at temperature T, describing scattering from E to E′, is given by

$\sigma \quad\left(\mathrm{E}^{\prime} \rightarrow \mathrm{E} ; \mathrm{T}\right)=\frac{\sigma_{b}}{T} \sqrt{\frac{\mathrm{E}}{\mathrm{E}^{\prime}}} e^{-\frac{\beta\left(\mathrm{E}^{\prime} \rightarrow \mathrm{E}\right)}{2}} \int \mathrm{P}\left(\mu_{0}\right) \mathrm{S}[\alpha, \beta ; T] \quad d \mu_{0}$

where β(E′→E) and α(E′,E,μ0) are dimensionless variables (functions of temperature) defining the energy and momentum exchange,  Batson Iii committed Jan 19, 2021 620 respectively, of the collision [BG70]; σb is the rigidly  Batson Iii committed Nov 17, 2020 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 bound scatter cross section, which is proportional to the free atom cross section; and S(α, β; T) describes the temperature-dependent thermal scattering law.

If atomic bonding effects are neglected, the atoms of a material behave like a gas in thermal equilibrium at the temperature of the medium. In this case S(α, β) can be expressed by an analytical function. CENTRM uses the free gas model for all nuclides except those materials that have thermal scattering laws available in the ENDF/B nuclear data files. The ENDF/B scattering law data account for the effects of molecular bonding and possible polyatomic crystalline structure. While free-gas kernels are computed internally in CENTRM, the kernel moments describing bound thermal scatterers are stored in a data file that can be accessed by CENTRM.

Bound thermal kernels

Thermal scattering from bound atoms is classified either as an “inelastic reaction,” in which the neutron energy may change, or an “elastic reaction,” in which the neutron changes direction, but does not change energy. In ENDF/B the former reactions are treated as incoherent inelastic scattering with a doubly differential kernel describing the secondary neutron energy and angle distribution. The latter reactions are usual treated as coherent elastic scatter characterized by diffractive interference of the scattered deBroglie waves, although a few materials are modeled by the incoherent elastic approximation. Legendre moments for thermal elastic kernels describe the secondary angular distribution with no energy exchange, at a given neutron energy. Bound scatter kernels have been processed by the AMPX code system for most of the ~25 compounds with thermal scatter laws in ENDF/B, and are stored in individual kinematics files distributed with the SCALE code system. These include materials such as: H in water, H and C in polyethylene, H and Zr in ZrH, C in graphite, deuterium in heavy water, Be metal, Be in BeO, etc. The CRAWDAD module processes scattering kernel data for individual nuclides into a combined library used in CENTRM, and also interpolates the kernels to the appropariate temperatures.

The bound scatter kernels are tabulated at different energy points from the flux solution mesh; therefore it is necessary to map the data onto the desired energy mesh in the CENTRM calculation. Because thermal elastic scattering results in no energy loss, the elastic moments only appear in the within-point term of the scattering source in the CENTRM thermal calculation. Thus the coherent elastic data is easily interpolated since it only involves a single energy index and temperature. However, the incoherent inelastic moments are 2-D arrays in terms of the initial and final energies, so that a 2-D interpolation must be done for each temperature. CENTRM uses a simple type of “unit-base transform” method to interpolate incoherent inelastic kernels onto the flux solution mesh. The method attempts to preserve the absolute peak of the secondary energy distribution, at given initial energy. For water-bound hydrogen and several other moderators, this is quite adequate, since the kernel generally has only a single maximum. However, if more than one local extrema is present, such as for graphite, the other local peaks are not explicitly preserved in the interpolation method. For this reason it is necessary to include a fairly dense set of initial energies in the tabulated kernels of graphite and similar materials, to avoid gross changes in the kernel shape at adjacent initial energy panels.

Free gas thermal kernels

CENTRM computes free-gas kernels using the approach proposed by Robinson [Rob81] as a modification to the original FLANGE [HF71] methodology. Legendre moments of the free-gas scatter kernel per unit lethargy are expressed as,

 Batson Iii committed Jan 19, 2021 685 (346)$\Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}\right)=\mathrm{A}^{(\mathrm{j})} \Sigma_{\mathrm{free}}^{(\mathrm{j})} \frac{\mathrm{E}}{\mathrm{E}^{\prime}} \mathrm{e}^{-\beta / 2} \sum_{\mathrm{n}=0} \mathrm{W}_{\mathrm{n}} \mathrm{H}_{\mathrm{n}}\left(\mathrm{E}^{\prime}, \mathrm{E}\right)$
 Batson Iii committed Nov 17, 2020 686 687 688 689 690 691 692 693 694 695 

where Wℓn are constant coefficients associated with the Legendre polynomial of order ℓ; Σfree is the constant free-atom cross section for the material; and Hn are the α-moments of the free-gas scatter law, given as

$\mathrm{H}_{\mathrm{n}}\left(\mathrm{E}^{\prime}, \mathrm{E}\right)=\frac{1}{\sqrt{\pi}} \int_{\alpha_{\mathrm{L}}}^{\alpha_{\mathrm{H}}} \alpha^{\mathrm{n}} \times\left(\frac{\mathrm{e}^{-\frac{\alpha^{2}+\beta^{2}}{4 \alpha}}}{2 \sqrt{\alpha}}\right) \mathrm{d} \alpha$

The limits on the above integral correspond to:

$\alpha_{\mathrm{L}}\left(\mathrm{E}^{\prime}, \mathrm{E}\right)=\alpha\left(\mathrm{E}^{\prime}, \mathrm{E}, \mu_{0}=-1\right) \quad \text { and } \quad \alpha_{\mathrm{H}}\left(\mathrm{E}^{\prime}, \mathrm{E}\right)=\alpha\left(\mathrm{E}^{\prime}, \mathrm{E}, \mu_{0}=1\right) .$

The alpha moments for n > 0 can be evaluated very efficiently using a  Batson Iii committed Jan 19, 2021 696 recursive relation [Wil00]:

 Batson Iii committed Nov 17, 2020 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 
$\mathrm{H}_{\mathrm{n}}\left(\mathrm{E}^{\prime}, \mathrm{E}\right)=2(2 \mathrm{n}-1) \mathrm{H}_{\mathrm{n}-1}+\beta^{2} \mathrm{H}_{\mathrm{n}-2}-\left[\mathrm{F}_{\mathrm{n}}\left(\sqrt{\alpha_{\mathrm{H}}}, \beta\right)-\mathrm{F}_{\mathrm{n}}\left(\sqrt{\alpha_{\mathrm{L}}}, \beta\right)\right]$

where Fn is the function,

$\mathrm{F}_{\mathrm{n}}(\mathrm{t}, \beta)=\frac{\mathrm{t}^{2 \mathrm{n}-1} \mathrm{e}^{-\frac{1}{4}\left(\frac{\beta^{2}}{\mathrm{t}^{2}}+\mathrm{t}^{2}\right)}}{\sqrt{\pi} / 2}$

Analytical expressions for the initial two moments, H0 and H:sub: −1, are given in [Rob81].

The standard free-gas kernel is based on the assumption of a constant free atom cross section. When averaged over the molecular velocity distribution, this gives a 1/v variation in the effective free-gas cross section at low energies. To approximately account for nuclear structure effects on the energy dependence of the thermal cross section (e.g., low energy resonances), the free-gas moments are multiplied by the ratio σs(E)/σFG(E), where σs is the Doppler broadened scatter cross section processed from ENDF/B data; and σFG is the effective free-gas cross section,

$\sigma_{\mathrm{FG}}\left(\mathrm{E}^{\prime}\right)=\frac{\sigma_{\mathrm{free}}}{\mathrm{y}^{2}}\left[\left(\mathrm{y}^{2}+1 / 2\right) \operatorname{erf}(\mathrm{y})+\frac{\mathrm{y} \mathrm{e}^{-\mathrm{y}^{2}}}{\sqrt{\pi}}\right]$

where $$\mathrm{y}^{2}=\mathrm{~A} \frac{\mathrm{E}}{\mathrm{kT}}$$.

Sub-moment expansion of the epithermal scattering source

One difficulty in computing the epithermal scatter source moments is  Batson Iii committed Jan 19, 2021 722 that the Legendre polynomial in the integrand of (342)  and (345) is a function  Batson Iii committed Nov 17, 2020 723 724 725 726 727 728 729 730 of both the initial and final lethargy (or energy) of the scattered neutrons, due to the correlation function G(j)(E,E′). At each lethargy u this requires that the u′-integral be recomputed over all lower lethargies, for every nuclide and moment. A more efficient algorithm would be possible if the differential scattering moments appearing in the integrand could be factored into a product of a function of u multiplied by a function of u′ such as

 Batson Iii committed Jan 19, 2021 731 (347)$\mathrm{S}_{\mathrm{k}}^{(\mathrm{j})}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}\right)=\mathrm{F}^{(\mathrm{j})}\left(\mathrm{u}^{\prime}\right) \quad \mathrm{H}^{(\mathrm{j})}(\mathrm{u})$
 Batson Iii committed Nov 17, 2020 732 733 734 735 736 737 738 739 

where F and H are the two factors (to be specified later).

If this is done, the u-function can be factored from the scatter source integrals, leaving only integrals over the u′-function as shown below:

$\mathrm{S}^{(\mathrm{j})}(\mathrm{u})=\int_{\mathrm{u}^{\prime}} \mathrm{S}_{\mathrm{k}}^{(\mathrm{j})}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}\right) \mathrm{du}^{\prime}=\mathrm{H}^{(\mathrm{j})}(\mathrm{u}) \int_{\mathrm{u}^{\prime}} \mathrm{F}^{(\mathrm{j})}\left(\mathrm{u}^{\prime}\right) \mathrm{du}^{\prime}$

Because the factored integrand does not depend on the variable u, a running summation over all u′ points can be accumulated and saved as the calculation sweeps from low to high lethargy. For example, note that the  Batson Iii committed Jan 19, 2021 740 ℓ = 0 moment in (342) is already separable into a product of u times u′  Batson Iii committed Nov 17, 2020 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 because P0 is equal to one at all values of G. Thus the isotropic component of the elastic differential scatter rate (per unit lethargy) from u′ to u is proportional to E/E′, where

$\mathrm{E}=\mathrm{E}(\mathrm{u})=\mathrm{E}_{\mathrm{ref}} \mathrm{e}^{-\mathrm{u}}, \quad \text { and } \quad \mathrm{E}^{\prime}=\mathrm{E}^{\prime}\left(\mathrm{u}^{\prime}\right)=\mathrm{E}_{\mathrm{ref}} \mathrm{e}^{-\mathrm{u}^{\prime}}$

Therefore, the two separable factors in the lowest moment, S(j)0.0(u′→u), are identified as,

$\begin{split}\begin{array}{l} \mathrm{H}(\mathrm{u})=\mathrm{E} /\left(1-\alpha^{(\mathrm{j})}\right), \quad \text { and } \\ \mathrm{F}\left(\mathrm{u}^{\prime}\right)=\Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime}\right) \Psi_{00}\left(\mathrm{u}^{\prime}\right) / \mathrm{E}^{\prime} \end{array}\end{split}$

However, the higher order Legendre moments contain the term P(G) in the integrand; and the expression for G(E′,E) is a difference of two terms that depend on both E and E′. A new method called a “sub‑moment expansion” has been developed for CENTRM that allows the Legendre polynomials appearing in the differential scatter moments to be factored into the desired separable form. Each spherical harmonic moment of the scattering source appears expanded in a series of factored “sub‑moments.”

The Legendre polynomial of order ℓ is a polynomial containing terms up to the ℓth power. Applying the binomial expansion theorem and some algebraic manipulation, the standard expression for P evaluated at “G” can be expressed as

 Batson Iii committed Jan 19, 2021 766 (348)$\mathrm{P}_{\ell}(\mathrm{G})=\frac{\mathrm{E}^{\prime}}{\mathrm{E}} \times \mathrm{a}_{1}^{\ell} \sum_{\mathrm{K}=-\ell}^{\ell} \tilde{\mathrm{g}}_{\ell, \mathrm{K}}(\mathrm{E}) \quad \mathrm{h}_{\mathrm{K}}(\mathrm{E}) \mathrm{h}_{\mathrm{K}}^{-1}\left(\mathrm{E}^{\prime}\right)$
 Batson Iii committed Nov 17, 2020 767 768 769 770 

where hk(E)=E1+K/2; and the expansion coefficients $$\tilde{\mathrm{g}}_{\ell, \mathrm{k}}$$ are equal to, $$\tilde{\mathrm{g}}_{\ell \mathrm{K}}=\frac{\mathrm{g}_{\ell \mathrm{K}}}{N_{\ell} \times \alpha_{1}^{\ell}}$$ where the gℓ,K (no tilde)  Batson Iii committed Jan 19, 2021 771 parameters were defined in [WA95] to be:

 Batson Iii committed Nov 17, 2020 772 
 Batson Iii committed Jan 19, 2021 773 (349)\begin{split}\begin{aligned}  Batson Iii committed Nov 17, 2020 774 775 776 777 778 779 780 &\mathrm{g}_{, \mathrm{K}}=\frac{\left(1+(-1)^{+\mathrm{K}}\right)}{2} \sum_{K^{\prime}=0}^{\frac{-K}{2}}(-1)^{K^{\prime}} b_{2 K^{\prime}+K},\left(\begin{array}{r} 2 K^{\prime}+K \\ K^{\prime} \end{array}\right) \quad a_{1}^{K+K^{\prime}} a_{2}^{K^{\prime}} ; \quad \text { for } \quad \mathrm{K} \geq 0\\ &\text { and }\\ &=\left(-a_{2} / a_{1}\right)^{|K|} \quad g_{,|K|} \quad ; \quad \text { for } \quad K<0 \end{aligned}\end{split}
 Batson Iii committed Jan 19, 2021 781 

In (348) –(349), the constants bm,ℓ and N are the standard  Batson Iii committed Nov 17, 2020 782 Legendre constants and normalization factors, respectively, which are  Batson Iii committed Jan 19, 2021 783 tabulated in Table 269 for orders through P7; and $$\left(\begin{array}{c}  Batson Iii committed Nov 17, 2020 784 785 786 787 788 \mathrm{m} \\ \mathrm{i} \end{array}\right)=$$ the binomial expansion coefficients(20) $$= \frac{\mathrm{m} !}{(\mathrm{m}-\mathrm{i}) ! \quad \mathrm{i} !}$$

 Batson Iii committed Jan 19, 2021 789