1. Background on the Linear Boltzmann Equation

The following textbooks are a good source of information on transport theory Duderstadt and Martin [DM79], computational methods for transport Lewis and Miller [LM84], and nuclear reactor theory Bell and Glasstone [BG79], Duderstadt and Hamilton [DH76]. For a review article on transport approximations, consult Sanchez and McCormick [SM82], Azmy and Sartori [AS10].

1.1. Definitions

In radiation transport, particles are described by their position in space and by their momentum vector. This leads to a six-dimensional phase-space. This increased dimensionality is what makes radiation transport problems computationally expensive.

The six-dimensional phase-space is often described in terms of the particle’s position $\overrightarrow{r}$, the particle’s direction of flight using the unit direction vector $\overrightarrow{\mathrm{\Omega }}$, and the particle’s energy $E$. Hence, to describe the time-dependent distribution of particles in phase-space, we introduce their phase-space density as:

$$n(\overrightarrow{r},\overrightarrow{\mathrm{\Omega }},E,t).$$

It is customary to introduce the following quantities:

• The angular flux $\psi$:

  $$\psi (\overrightarrow{r},\overrightarrow{\mathrm{\Omega }},E,t)=v(E)n(\overrightarrow{r},\overrightarrow{\mathrm{\Omega }},E,t),$$

  where $v(E)$ is the magnitude of the particle’s velocity.

• The angular current $\overrightarrow{j}$:

  $$\overrightarrow{j}(\overrightarrow{r},\overrightarrow{\mathrm{\Omega }},E,t)=\overrightarrow{\mathrm{\Omega }}\psi (\overrightarrow{r},\overrightarrow{\mathrm{\Omega }},E,t).$$

The unit direction vector in $xyz$ in Cartesian coordinate system is

$$\begin{array}{r}\overrightarrow{\mathrm{\Omega }}=\left[\begin{array}{c}\sin \theta \cos \phi \\ \sin \theta \sin \phi \\ \cos \theta \end{array}\right]\end{array}$$

where $\theta$ is the polar angle and $\phi$ the azimuthal angle. We often use $\mu =\cos \theta$.

1.2. The Linear Boltzmann Equation

The linear Boltzmann transport equation for neutral particles is a conservation statement over the phase-space, stating that the rate of change of the particle’s density is equal to gain terms minus loss terms. It is typically written as follows:

$$\begin{array}{r}\begin{array}{c}\frac{1}{v(E)}\frac{\partial \psi (\overrightarrow{r},\overrightarrow{\mathrm{\Omega }},E,t)}{\partial t}=\\ \int d{E}^{\prime }{\int }_{4\pi }d{\mathrm{\Omega }}^{\prime }{\sigma }_s(\overrightarrow{r},{\overrightarrow{\mathrm{\Omega }}}^{\prime }\to \overrightarrow{\mathrm{\Omega }},E^{\prime }\to E,t)\psi (\overrightarrow{r},{\overrightarrow{\mathrm{\Omega }}}^{\prime },E^{\prime },t)+Q_{\text{ext}}(\overrightarrow{r},\overrightarrow{\mathrm{\Omega }},E,t)\\ -\overrightarrow{\mathrm{\Omega }}\cdot \overrightarrow{\mathrm{\nabla }}\psi (\overrightarrow{r},\overrightarrow{\mathrm{\Omega }},E,t)-{\sigma }_t(\overrightarrow{r},E,t)\psi (\overrightarrow{r},\overrightarrow{\mathrm{\Omega }},E,t)\end{array}\end{array}$$

for $\overrightarrow{r}\in \mathcal{D}$ (the spatial domain), $E\in \mathcal{E}=[0,E_{\text{max}}]$ (the energy range), and $\overrightarrow{\mathrm{\Omega }}\in {\mathcal{S}}^2$ (the unit sphere).

• $Q_{\text{ext}}$ represents the external source rate density,

• ${\sigma }_t(\overrightarrow{r},E,t)$ is the total interaction cross section,

• ${\sigma }_s(\overrightarrow{r},{\overrightarrow{\mathrm{\Omega }}}^{\prime }\to \overrightarrow{\mathrm{\Omega }},E^{\prime }\to E,t)$ is the double differential scattering cross section; note that additional energy-angle redistribution events can be modeled, such as $(\text{n},\text{xn})$ neutron interactions; finally, note that for isotropic materials (which is an assumption we make), the angular distribution of the outgoing particles only depends on the cosine of the incoming and outgoing directions, ${\mu }_0={\overrightarrow{\mathrm{\Omega }}}^{\prime }\cdot \overrightarrow{\mathrm{\Omega }}$, so we will replace ${\sigma }_s(\overrightarrow{r},{\overrightarrow{\mathrm{\Omega }}}^{\prime }\to \overrightarrow{\mathrm{\Omega }},E^{\prime }\to E,t)$ with ${\sigma }_s(\overrightarrow{r},{\overrightarrow{\mathrm{\Omega }}}^{\prime }\cdot \overrightarrow{\mathrm{\Omega }},E^{\prime }\to E,t)$.

In succinct operator notation, the above equation can be re-written as

$$\frac{1}{v}\frac{\partial \mathrm{\Psi }}{\partial t}=H\mathrm{\Psi }+Q_{\text{ext}}-L\mathrm{\Psi },$$

where $L$ is the streaming + interaction operator and $H$ is the collision operator.

If fission interactions and delayed neutron production are to be included, the above equation is amended as follows

$$\frac{1}{v}\frac{\partial \mathrm{\Psi }}{\partial t}=H\mathrm{\Psi }+P_p\mathrm{\Psi }+Q_{\text{ext}}+\sum _i{S}_{d,i}-L\mathrm{\Psi },$$

where the prompt fission operator is

$$P_p\mathrm{\Psi }=\frac{1}{4\pi }\int d{E}^{\prime }{\int }_{4\pi }d{\mathrm{\Omega }}^{\prime }{\nu }_p{\sigma }_f(\overrightarrow{r},E^{\prime })p_f(E^{\prime }\to E,t)\psi (\overrightarrow{r},{\overrightarrow{\mathrm{\Omega }}}^{\prime },E^{\prime },t),$$

(often, the dependence of $p_f(E’to E)$ is assumed to be weak in the incident energy $E’$ and thus $p_f(E’to E)approx p_f(E):=chi_p(E)$, where $chi_p$ is known as the prompt fission spectrum).

and the delayed neutron operator is

$$S_{d,i}=\frac{{\chi }_{d,i}(E)}{4\pi }{C}_i(\overrightarrow{r},t),$$

where $C_i$, the neutron precursor concentration in delayed group $i$, satisfies its own conservation law.

1.3. Boundary Conditions

Boundary conditions are usually supplied either as a known incoming flux (including a zero function for vacuum boundaries) or as an albedo condition. We denote by $\mathrm{\Gamma }$ the boundary of the spatial domain $\mathcal{D}$ (that is, $\mathrm{\Gamma }=\partial \mathcal{D}$).

The following types of boundary conditions are employed:

• The imposed-flux conditions can be expressed as

  $$\psi (\overrightarrow{r},\overrightarrow{\mathrm{\Omega }},E,t)={\psi }_{\text{inc}}(\overrightarrow{r},\overrightarrow{\mathrm{\Omega }},E,t)\mathrm{\forall }\overrightarrow{r}\in {\mathrm{\Gamma }}^{-}$$

  where

  $${\mathrm{\Gamma }}^{-}=\{\overrightarrow{r}\in \mathrm{\Gamma }\text{ such that }\overrightarrow{\mathrm{\Omega }}\cdot \overrightarrow{n}(\overrightarrow{r})<0\}$$

  with $\overrightarrow{n}(\overrightarrow{r})$ the outward pointing unit normal vector at position $\overrightarrow{r}$.

• The albedo condition can be expressed as

  $$\psi (\overrightarrow{r},\overrightarrow{\mathrm{\Omega }},E,t)=\mathcal{A}(\overrightarrow{r},{\overrightarrow{\mathrm{\Omega }}}^{\prime }\to \overrightarrow{\mathrm{\Omega }},E^{\prime }\to E,t)\psi (\overrightarrow{r},{\overrightarrow{\mathrm{\Omega }}}^{\prime },E^{\prime },t)$$

  with

  • $\mathcal{A}(\overrightarrow{r},{\overrightarrow{\mathrm{\Omega }}}^{\prime }\to \overrightarrow{\mathrm{\Omega }},E^{\prime }\to E,t)$ the albedo operator, and

  • ${\overrightarrow{\mathrm{\Omega }}}^{\prime }\cdot \overrightarrow{n}(\overrightarrow{r})>0$ (outgoing direction) and the incoming direction is $\overrightarrow{\mathrm{\Omega }}={\overrightarrow{\mathrm{\Omega }}}^{\prime }-2(\overrightarrow{\mathrm{\Omega }}\cdot \overrightarrow{n}(\overrightarrow{r}))\overrightarrow{n}(\overrightarrow{r})$.

For $k$-eigenvalue problems, the boundary conditions usually devolve to a zero-incoming flux (${\psi }_{\text{inc}}=0$) or a unity albedo (i.e., symmetry line, with $\mathcal{A}=1$).

1.4. Initial Conditions

For time-dependent problems, initial conditions are supplied as

$$\psi (\overrightarrow{r},\overrightarrow{\mathrm{\Omega }},E,t=0)=f_0(\overrightarrow{r},\overrightarrow{\mathrm{\Omega }},E)\mathrm{\forall }\overrightarrow{r}\in \mathcal{D},\mathrm{\forall }E\in \mathcal{E},\mathrm{\forall }\overrightarrow{\mathrm{\Omega }}\in {\mathcal{S}}^2$$

1.5. Expansion of the Angle Redistribution Term

The angle redistribution term

$${\int }_{4\pi }d{\mathrm{\Omega }}^{\prime }{\sigma }_s(\overrightarrow{r},{\overrightarrow{\mathrm{\Omega }}}^{\prime }\cdot \overrightarrow{\mathrm{\Omega }},E^{\prime }\to E,t)\psi (\overrightarrow{r},{\overrightarrow{\mathrm{\Omega }}}^{\prime },E^{\prime },t)$$

can be expanded in angle on the unit sphere ${\mathcal{S}}^2$ using (real-valued) spherical harmonic functions to yield

$$\sum _{\ell =0}^{L_{\text{max}}}\sum _{m=-\ell }^{m=\ell }\frac{2\ell +1}{4\pi }{\sigma }_{s,\ell }(\overrightarrow{r},E^{\prime }\to E,t)Y_{\ell ,m}(\overrightarrow{\mathrm{\Omega }}){\varphi }_{\ell ,m}(\overrightarrow{r},E^{\prime },t)$$

where

• $L_{\text{max}}$ is the highest order of scattering anisotropy retrained in the cross section expansion (usually, a user-supplied value)

• the Legendre moments of the scattering cross section are defined as

  $${\sigma }_{s,\ell }(\overrightarrow{r},E^{\prime }\to E,t)=2\pi {\int }_{-1}^{1}d{\mu }_0{\sigma }_s(\overrightarrow{r},{\mu }_0,E^{\prime }\to E,t)P_{\ell }({\mu }_0)$$

  where ${\mu }_0={\overrightarrow{\mathrm{\Omega }}}^{\prime }\cdot \overrightarrow{\mathrm{\Omega }}$ and $P_{\ell }(\mu )$ is the Legendre polynomial of degree $\ell$.

• the moment of the angular flux are defined as

  $${\varphi }_{\ell ,m}(\overrightarrow{r},E,t)={\int }_{4\pi }d\mathrm{\Omega }{Y}_{\ell ,m}(\overrightarrow{\mathrm{\Omega }})\psi (\overrightarrow{r},\overrightarrow{\mathrm{\Omega }},E,t)$$

• the real-valued spherical harmonics are

  $$\begin{array}{r}{Y}_{\ell ,m}(\overrightarrow{\mathrm{\Omega }})=\{\begin{array}{ll}(-1{)}^m\sqrt{2}\sqrt{\frac{(2\ell +1)}{4\pi }\frac{(\ell -|m|)!}{(\ell +|m|)!}}{P}_{\ell }^{|m|}(\cos \theta )\sin |m|\phi & \text{if }m<0\\ \\ \sqrt{\frac{(2\ell +1)}{4\pi }}{P}_{\ell }^m(\cos \theta )& \text{if }m=0\\ \\ (-1{)}^m\sqrt{2}\sqrt{\frac{(2\ell +1)}{4\pi }\frac{(\ell -m)!}{(\ell +m)!}}{P}_{\ell }^m(\cos \theta )\cos m\phi & \text{if }m>0\end{array}\end{array}$$

  where $P_{\ell }^m(\mu )$ is the associated Legendre function of degree $\ell$ and of order $m$.

In operator notation, the steady-state, source-driven problem $L\mathrm{\Psi }=H\mathrm{\Psi }+Q_{\text{ext}}$ can now be re-written as

$$L\mathrm{\Psi }=M\mathrm{\Sigma }\mathrm{\Phi }+Q_{\text{ext}}\text{with }\mathrm{\Phi }=D\mathrm{\Psi }$$

where

• $\mathrm{\Phi }$ are the flux moments,

• $D$ is the discrete-to-moment operator and denotes the angular integration of the angular flux to yield the flux moment ($\mathrm{\Phi }=D\mathrm{\Psi }$). The term discrete stems from the fact that the integration is performed using a quadrature rule, hence at discrete directions for the angular flux,

• $\mathrm{\Sigma }$ denotes the matrix containing the Legendre moments of the scattering cross section, and

• $M$ denotes the moment-to-discrete operator, which takes the source moments ($\mathrm{\Sigma }\mathrm{\Phi }$) and evaluates that term in direction.

If we denote by $N_{\text{mom}}$ the maximum number of flux moments, we have

$$\begin{array}{r}\begin{array}{rlr}{N}_{\text{mom}}& =L_{\text{max}}+1& \text{in 1D,}\\ N_{\text{mom}}& =\frac{(L_{\text{max}}+1)(L_{\text{max}}+2)}{2}& \text{in 2D,}\\ N_{\text{mom}}& =(L_{\text{max}}+1{)}^2& \text{in 3D.}\end{array}\end{array}$$

With the introduction of the moment-to-discrete operator, we can also update the fission production operator

$$\text{from}P\mathrm{\Psi }\text{to}MF\mathrm{\Phi }$$

where, hereafter, $F$ will denote the fission operator acting on flux moments.

1.6. Short summary of Transport Equations Solved in OpenSn

• Time-dependent transport problem:

  $$\frac{1}{v}\frac{\partial \mathrm{\Psi }}{\partial t}=M\mathrm{\Sigma }\mathrm{\Phi }+M{F}_p\mathrm{\Phi }+Q_{\text{ext}}+\sum _i{S}_{d,i}-L\mathrm{\Psi },$$

• Steady-state, subcritical, source-driven transport problem:

  $$L\mathrm{\Psi }=M\mathrm{\Sigma }\mathrm{\Phi }+MF\mathrm{\Phi }+Q_{\text{ext}},$$

  where $F$ is the total (prompt+delayed) fission production operator.

• $k$-eigenvalue transport problem:

  $$L\mathrm{\Psi }=M\mathrm{\Sigma }\mathrm{\Phi }+\frac{1}{k_{\text{eff}}}MF\mathrm{\Phi }$$

1.7. Streaming Term in Cartesian and Curvilinear Coordinate Systems

In $xyz$ coordinates, the streaming term is given by

$$\overrightarrow{\mathrm{\Omega }}\cdot \overrightarrow{\mathrm{\nabla }}\psi ={\mathrm{\Omega }}_x{\partial }_x\psi +{\mathrm{\Omega }}_y{\partial }_y\psi +{\mathrm{\Omega }}_z{\partial }_z\psi$$

In $r\theta z$ cylindrical coordinates, the streaming term is given by

$$\overrightarrow{\mathrm{\Omega }}\cdot \overrightarrow{\mathrm{\nabla }}=\frac{\xi }{r}{\partial }_r(r\psi )-\frac{1}{r}{\partial }_{\phi }(\eta \psi )+\mu {\partial }_z\psi$$

with $\xi =\overrightarrow{\mathrm{\Omega }}\cdot {\overrightarrow{e}}_r=\sqrt{1-{\mu }^2}\cos \phi$, $\eta =\overrightarrow{\mathrm{\Omega }}\cdot {\overrightarrow{e}}_{\theta }=\sqrt{1-{\mu }^2}\sin \phi$, $\mu =\overrightarrow{\mathrm{\Omega }}\cdot {\overrightarrow{e}}_z$

1.8. References

[DM79]

James J Duderstadt and William Russell Martin. Transport theory. John Wiley & Sons, 1979.

[LM84]

Elmer Eugene Lewis and Warren F Miller. Computational methods of neutron transport. John Wiley and Sons, Inc., New York, NY, 1984. ISBN 0-89448-452-4.

[BG79]

George I Bell and Samuel Glasstone. Nuclear reactor theory. RE Krieger Publishing Company, 1979. ISBN 978-0-442-20684-0.

[DH76]

James J Duderstadt and Louis J Hamilton. Nuclear reactor analysis. Wiley, 1976.

[SM82]

R Sanchez and Norman J McCormick. Review of neutron transport approximations. Nuclear Science and Engineering, 1982.

[AS10]

Y. Azmy and E. Sartori. Nuclear Computational Science: A Century in Review. Springer Netherlands, 2010. ISBN 9789048134113. URL: https://books.google.com/books?id=4imnyVPQ-0AC.