Mathematical formulae for neutron self-shielding properties of media in an isotropic neutron field

Ateia W. Mahmoud [email protected] [email protected] Physics Department, Faculty of Science, Ain Shams University, Cairo, Egypt. Reactor Physics Department, Nuclear Research Center, Egyptian Atomic Energy Authority, Cairo 13759, Egypt. Elsayed K. Elmaghraby [email protected] [email protected] Experimental Nuclear Physics Department, Nuclear Research Center, Egyptian Atomic Energy Authority, Cairo 13759, Egypt. E. Salama [email protected] Basic Science Department, Faculty of Engineering, The British University in Egypt (BUE), Cairo, Egypt. A. Elghazaly [email protected] Reactor Physics Department, Nuclear Research Center, Egyptian Atomic Energy Authority, Cairo 13759, Egypt. S. A. El-fiki [email protected] Physics Department, Faculty of Science, Ain Shams University, Cairo, Egypt.

(June 30, 2021)

Abstract

The complexity of the neutron transport phenomenon throws its shadows on every physical system wherever neutron is produced or used. In the current study, an ab initio derivation of the neutron self-shielding factor to solve the problem of the decrease of the neutron flux as it penetrates into a material placed in an isotropic neutron field. We have employed the theory of steady-state neutron transport, starting from Stuart’s formula. Simple formulae were derived based on the integral cross-section parameters that could be adopted by the user according to various variables, such as the neutron flux distribution and geometry of the simulation at hand. The concluded formulae of the self-shielding factors comprise an inverted sigmoid function normalized with a weight representing the ratio between the macroscopic total and scattering cross-sections of the medium. The general convex volume geometries are reduced to a set of chord lengths, while the neutron interactions probabilities within the volume are parameterized to the epithermal and thermal neutron energies. The arguments of the inverted-sigmoid function were derived from a simplified version of neutron transport formulation. Accordingly, the obtained general formulae were successful in giving the values of the experimental neutron self-shielding factor for different elements and different geometries.

Neutron self-shielding; Neutron transport and absorption; ab initio approach.

pacs:

28.20.Gd , 25.40.Dn , 25.40.Ep , 25.40.Fq , 28.20.-v , 28.20.Cz , 28.41.-i , 28.41.Pa , 29.25.Dz

I Introduction

Over time, neutron activation analysis has been evolving into a very effective nuclear analytical technique. Such techniques are often utilized for non-destructive elemental concentration measurement in unknown materials and nuclear material interrogation [1, 2]. The constraints include neutron fluence, the fraction of fluence reaches the interior of the sample, sample mass and sample geometry [3, 4, 5, 6, 7, 8, 9, 10, 11, 12]. Furthermore, the neutron’s energy spectrum is varied, but ideally suited to research using the ideal Maxwellian distribution at room temperature, while the other distributions must be altered to match the reference nuclear reaction data [7, 13, 3]. Apart from that, neutrons are deeply employed in two significant geometries, including but not limited to: (1) beam geometry, where the neutron currents are assumed to travel in one direction, and (2) field geometry, where neutrons impact the sample from all directions presuming the material is isotropic. There is an important functional difference between these two geometries, i.e. the effect on the neutron flux itself. For instance, when exposing a sample to a neutron beam, the interior of the sample will be exposed to a lesser neutron fluence than the exterior part, in all circumstances regardless of the geometry of the neutron current. This phenomenally known as self-shielding, and it is a critical element of the neutron transport phenomena. In the case of field geometry, the net neutron current essentially disappears, while the fluence (or flux) becomes the observable quantity. There is an interplay between neutron absorption in the sample and the overall neutron flux [14]. Predominantly, the correlation between neutrons self-shield factors and the set of parameters involved in the calculation of its value had been studied by several scientists [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25], who gave dimensionless variables to identify and encompass the physical and geometric varieties of the samples geometries in order to attain a universal formula for self-shielding. The Montè-Carlo approach effectively calculates self-shielding, but it takes time and an experienced user to achieve acceptable accuracy and efficiency [26], see Appendix A. Empirical expressions, such as those given by researchers in Refs. [17, 18, 20, 21, 19] based on Ref. [15], had became routine in calculating self-shielding, these empirical expressions have been derived for a few specific geometries and limited number of elements.

Herein, we present a complete investigation of the neutron self-shielding phenomenon in different media. Additionally, we had provided a full description of the physics behind the theme, taking into consideration the neutron transport inside the sample, and the absorption and scattering phenomena as a function of neutron energy. We aim to transform the problem from a spectroscopic set of parameters, usually unreachable for the common user, to an integrated set of well-known parameters and factors. The use of detailed spectroscopic parameters, such as ENDF data, cross-section, detailed dimensions and shape, the widths of neutron resonances for either scattering or absorption, etc., requires time and experienced users to make use of them with acceptable accuracy and efficiency, the cost most scientists cannot afford to just calculate a single parameter in their routine work. Our intention is focused on avoiding such cost and enhancing present existing formulae and using of an integral set of the well-known parameters, such as thermal cross-section, resonance integral, average chord length. All remaining factors are calculated from these three parameters. Though, existence of a mathematical formulation of self-shielding in material of different geometries and composition shall deliver additional tool to improve precision of reaction parameters and activation analysis calculations.

II Materials and Methods

Experimentally measured and theoretically calculated data were collected from different sources for the self-shielding factor in In, Au, Co, Cu, and Fe samples. The geometries for these elements were foils, wires, and infinite slabs. The Experimental data of G_th of Mahmoud et al. [27], Taylor & Linacre [28], Carre et al. [29], Hasnain et al. [30], Sola [31], Walker et al. [32], Klema [33], and Crane $\&$ Doerner [34] were digitized from Martinho et al. [20]. The data of Eastwood & Werner [35] for Co Wire was collected from their original values.

For the epithermal neutron energy region, The experimental results from literature of Gonalves et al. [17], Lopes [36, 37], McGarry [38], Brose [39], Yamamoto et al. [40], Jefferies et al. [41], Eastwood $\&$ Werner [35], and Kumpf [42] were used.

Having two energy ranges at thermal neutron region and epithermal neutron energies, the cross-section was taken as the element-averaged thermal neutron cross-section and the resonance integrals for 1/E averaged neutron distribution as given in Table 1 based on previous evaluations [10, 43, 9]. These data shall be used in the specific calculations presented in the present work.

Table 1: Element-averaged cross- section and resonance integral data.

\sigma_{g}

and

\sigma_{s}

, are the capture and scattering cross-section at thermal energies,

\sigma_{a}

\sigma_{t}

are the absorption and total cross-sections at thermal energies. Similarly,

I_{g}

I_{s}

I_{a}

and

I_{t}

have the same respective meaning of its subscript. Data were taken from [10, 9].

Isotope	$\sigma_{g}$	$\sigma_{s}$	$\sigma_{a}$	$\sigma_{t}$	$I_{g}$	$I_{s}$	$I_{a}$	$I_{t}$
	(b)	(b)	(b)	(b)	(b)	(b)	(b)	(b)
Na	0.528	3.3929	0.528	3.9209	0.3021	130.81	0.3021	131.11
Mn	13.275	2.1163	13.275	15.391	13.168	621.33	13.168	634.5
Fe	2.5615	11.35	2.5615	13.912	1.2706	127.09	1.2706	128.36
Co	37.173	6.0319	37.173	43.204	74.78	791.53	74.78	866.31
Cu	3.7531	7.8424	3.7531	11.595	4.0309	129.89	4.0309	133.93
In	194.07	2.5686	194.07	196.64	3088.5	214.12	3088.5	3302.6
Au	98.672	7.9298	98.672	106.6	1567.9	405.52	1567.9	1973.4

The uncertainty of digitized data was difficult to determine due to the use of different linear and logarithmic scales in old graphs and dependence among them. We had used the following formula $\sigma=$ $\sigma_{X}$ + $\sigma_{Y}$ , with $\sigma_{X}$ and $\sigma_{Y}$ are the dependent uncertainties in the digitized X and Y coordinates in the graph. The typical value of digitization uncertainty was less than 2% which was added to the reported uncertainty, if available.

III Results and discussion

According to Stuart [44], the old quantity for an absorbing body is its neutron blackness,

\beta=\frac{j_{in}-j_{out}}{j_{in}},

(1)

based on the neutron current density entering the body ( $j_{in}$ ) or going out of it ( $j_{out}$ ). Stuart [44] had derived the formula for $\beta$ based on variational principle and assuming uniform isotropic neutron field with scattering that do not change energy spectrum on the neutrons (change of energy is treated as absorption). Blaauw [16, 45] began with Stuart’s formula [44];

\beta=\frac{\frac{\Sigma_{a}}{\Sigma_{t}}P_{0}}{1-\frac{\Sigma_{s}}{\Sigma_{t}}\left(1-\frac{1}{\Sigma_{t}\bar{\ell}}\,P_{0}\right)}

(2)

Here, $P_{0}$ is the probability of the first interaction derived from the transport kernal [45, 44] in steady-state.

P_{0}=\frac{\sigma_{t}}{S}\int_{V}\Pi_{0}(\vec{r})d\vec{r},

(3)

where $\Pi_{0}(\vec{r})$ is the unscattered flux within the material;

\Pi_{0}(\vec{r})=4\int_{S}(\vec{n}\cdot\vec{\lambda})\mathbb{G}(\vec{r},\vec{r}^{\prime})dS,

(4)

Here, $\vec{n}$ and $\vec{\lambda}$ are unit vectors in the direction of the normal to the surface and the neutron wavevector, respectively. The point symmetry neutron collision kernel has the form of Green’s function [46, 47, 48, 25]:

\mathbb{G}(\vec{r},\vec{r}^{\prime})=\frac{\exp\left(-\Sigma_{t}|\vec{r}-\vec{r}^{\prime}|\right)}{4\pi|\vec{r}-\vec{r}^{\prime}|^{2}},

(5)

that gives the probability a neutrons shifts between phase-space coordinates $\vec{r}$ and $\vec{r}^{\prime}$ in one collision in point geometry. Time reversal applied in such cases by interchange or $\vec{r}$ and $\vec{r}^{\prime}$ . In general, other geometries had asymptotic form as point like geometry [47, 25]. The value of $\Pi_{0}(\vec{r})$ equals to 4 in case of nonexistence of the material in the medium because in Eq. 4 become unity. And that is the condition that must be satisfied by the transport kernels in all geometries. Multiple collision probability may be obtained through recursive relation [44]

P_{n}=1-\frac{\int_{V}\Pi_{n-1}(\vec{r})\Pi_{n}(\vec{r})d\vec{r}}{4\int_{V}\Pi_{n-1}(\vec{r}))d\vec{r}},

(6)

\Pi_{n}=\frac{1}{\bar{\ell}}\int_{S}\Pi_{n-1}(\vec{r^{\prime}})\mathbb{G}(\vec{r},\vec{r}^{\prime})dS,

(7)

The value of $\bar{\ell}$ is the average of Chord Length Distribution (CLD) may be weighted with the cosine of the angle between chord and the normal to the surface,

\bar{\ell}=\frac{\int_{0}^{\pi}\int_{0}^{\pi/2}V(\underline{r},\underline{\theta},\underline{\phi})\cos(\underline{\theta})d\underline{\theta}d\underline{\phi}}{\int_{0}^{\pi}\int_{0}^{\pi/2}S_{\bot}(\underline{r},\underline{\theta},\underline{\phi})\cos(\underline{\theta})d\underline{\theta}d\underline{\phi}}.

(8)

Here, $S_{\bot}$ is the surface area perpendicular to the direction of the neutron. The value of mean CLD for a convex body is related to the first Cauchy formula [49] of the integration for cylindrical shape comprises 4V/S [50, 51], where $S$ is the surface area inclosing a volume $V$ . This value may be used as the upper limit of one of the integrations. Note that the coordinate variables were underlined in order to avoid confusion among symbols.

There was equivalent definition of the sample blackness, the self-shielding factor denoted $G_{(\mathrm{energy~{}domain})}$ ; which is defined as the ratio between the volume-averaged fluence rate within the material’s volume that may absorb or scatter neutrons and the fluence rate within the same volume considering absence of the interaction with neutrons. According to Blaauw [16],

G_{(\mathrm{energy~{}domain})}=\frac{\int_{V}\psi(\vec{r})d\vec{r}}{V}+(\mathrm{higher~{}order~{}terms})

(9)

the higher order terms were introduced by Blaauw [16] for extended neutron velocity distributions – denoted here $\mathcal{R}$ . Remembering that the $G_{(\mathrm{energy~{}domain})}$ is neutron energy specific parameter, any perturbation of the neutron energy distribution shall affect the experimental results as discussed in earlier work [52].

III.1 Mathematical model

In accordance to previous constraint of self-shielding formulae, we shall use Eq. 2, include the high order in Eq. 9, and use the relation of blackness and self-shielding [45] (i.e. $\Sigma_{a}\bar{\ell}G$ = $\beta$ as first approximation). The combined formulae comprises;

\begin{split}G_{(\mathrm{energy~{}domain})}&=\frac{\left(\frac{\Sigma_{t}}{\Sigma_{s}}\right)}{1+\Sigma_{t}\frac{\bar{\ell}}{P_{0}}\left(\frac{\Sigma_{a}}{\Sigma_{s}}\right)}+\mathcal{R}.\end{split}

(10)

The value of $\mathcal{R}$ was found to has negligible contributions except when relying on the entire range of Maxwellian thermal neutron distribution, as proven in Appendix B. However, and for the practical of constraining the thermal neutron energy range by the cadmium cutoff energy around 0.5 eV, this term can be neglected. Hence,

G_{(\mathrm{energy~{}domain})}=\left(\frac{\Sigma_{t}}{\Sigma_{s}}\right)\frac{1}{1+\Sigma_{t}\frac{\bar{\ell}}{P_{0}}\left(\frac{\Sigma_{a}}{\Sigma_{s}}\right)}

(11)

The value of the first term in Eq. 11 is correct when it less than 1, i.e. under the condition:

\left(\frac{\Sigma_{t}}{\Sigma_{s}}\right)\leq 1+\Sigma_{t}\frac{\bar{\ell}}{P_{0}}\left(\frac{\Sigma_{a}}{\Sigma_{s}}\right)

(12)

\Sigma_{t}\leq\Sigma_{s}+\frac{\Sigma_{t}\bar{\ell}}{P_{0}}\Sigma_{a}

(13)

i.e. the parameters in the first term in Eq. 11 must satisfy the condition.

\frac{\Sigma_{t}\bar{\ell}}{P_{0}}\geq 1

(14)

Precise choice of the value of $P_{0}$ is given in Section III.3.

Eq. 11 can be rewritten as follows:

G_{(\mathrm{energy~{}domain})}=\left(\frac{\Sigma_{t}}{\Sigma_{s}}\right)\times\frac{1}{1+\mathcal{Z}}

(15)

where

\mathcal{Z}=\underbrace{\Omega(\bar{\ell},\Sigma_{a},\Sigma_{s},\Sigma_{t})}_{\mathrm{Geometry}}\underbrace{\chi(\Sigma_{t})}_{\mathrm{Composition}}\underbrace{\eta(\Sigma_{a},\Sigma_{s})}_{\mathrm{Probability}}.

(16)

The factor $\left(\frac{\Sigma_{t}}{\Sigma_{s}}\right)$ represents the weight of the total interaction cross-section versus the scattering contribution. The dimensionless parameter, $\mathcal{Z}$ , is expressed as product of three functions, geometry function ( $\Omega=\cfrac{\bar{\ell}}{P_{0}}$ ) as a function of the dimensions of the sample in the unit of [cm], macroscopic cross-section function ( $\chi=\Sigma_{t}$ ) in the units of [cm^-1] which depends on the isotopic content of the sample, and a dimensionless neutron energy correcting factor ( $\eta=\cfrac{\Sigma_{a}}{\Sigma_{s}}$ ) which is a function of the neutron absorption and the scattering cross-sections. The neutron-chord length is not the only parameter in the transport equation that depends on geometry; the first and higher order interaction probabilities, i.e. $P_{0}$ , are also dependent on both geometry and medium contents. These are the fundamental morphological descriptors of the media that describes the mean intercept length and relative to the mean free-paths of neutron within the medium. Here, $\Omega$ contains the factor of a shift in the Euclidean distance value and reflecting the total distance of interest within the medium, while $\chi$ and $\eta$ determine the slope of the steeping part of the curve.

Macroscopic cross-section function is expressed as

\chi(E_{n})=\Sigma(E_{1},E_{2})=\frac{\rho N_{A}\theta_{i}}{M}\tilde{\sigma},

(17)

\tilde{\sigma}=\frac{\int_{E_{1}}^{E_{2}}\sigma_{c,\,i}(E_{n})\varphi(E_{n})dE_{n}}{\int_{E_{1}}^{E_{2}}\varphi(E_{n})dE_{n}}

(18)

where $N_{A}$ is the Avogadro’s number [mol^-1], $\rho$ is the density of the material [g cm^-3], and $M$ is its atomic mass [g mol^-1]. $\theta_{i}$ is the isotopic abundance of the absorbing isotope in which it should be multiplied by the fraction of the element in the material if a compound material is used. Here, $\tilde{\sigma}$ is the integral cross-section in the energy domain between $E_{1}$ and $E_{2}$ . For practical purposes, the thermal neutron energy range is bounded by the cadmium cutoff energy around 0.5 eV while the epithermal range extends from 0.5 to few MeVs [3, 10, 7]. Here, $\sigma_{c_{i}}(E_{n})$ is the $i^{th}$ -isotope’s cross-section for the reaction channel $c$ at the neutron energy $E_{n}$ in this energy domain [cm²]. In case of compounds, this formula becomes a summation over $i^{th}$ -isotope.

The $\eta$ term is the scattering to absorption ratio;

\eta(\Sigma_{s},\Sigma_{a})=\frac{\Sigma_{a}}{\Sigma_{s}}

(19)

In general, values of thermal cross-section and resonance integral are well known from tables [10, 53, 54] or integration of spectroscopic cross-section [9].

III.2 The geometry factor

The geometry factor depends on the neutron-chord length and the probability of interaction as;

\Omega(\text{shape parameters})=\frac{\bar{\ell}}{P_{0}}.

(20)

There were great efforts to parameterize this factor through years. Recent efforts had been made by Trkov et al. [51] especially in the extended range of the neutron resonances.

In integral geometry, obtaining the orientation-dependent chord lengths of a convex body, in general, is a complicated mathematical argument; most researchers treat the problem with the body that has the minimal volume in a class of convex bodies having the same dimensions [55]. As long as we want to escape the rigorous derivations of the average neutron-chord length (cf. Refs. [56, 57, 50, 58, 59, 60] for details), we shall use simple formulation of the average neutron-chord length based on the fact that the trajectories of the incident isotropic neutrons traverse different lengths within the body due to scattering. The derivation proceeded under this assumption in which the convex bodies are in an isotropic neutron field are the geometry which the most materials have.

Considering an Euclidean space ( $\mathbb{E}^{3}$ ) where the irradiated body is located with a center-of-mass at the origin of the coordinate system. The three coordinate vectors are in orthogonal directions – denoted $\hat{1}$ , $\hat{2}$ , and $\hat{3}$ . The orientations of these coordinates were chosen as follows: At least one of these vectors ( $\hat{1}$ ) shall intersects the body surface in direction of the shortest length between the center-of-mass (at the $\mathbb{E}^{3}$ origin) and a point on the surface of the body, see Fig. 1.

Refer to caption — Figure 1: (Color on-line) Most known shapes of neutron activated materials. $L_{1}$ , $L_{2}$ , and $L_{3}$ are the measurable lengths, see text.

The next coordinate vector ( $\hat{2}$ ) shall lay in a plan perpendicular to the first one to the shortest point on the body surface. The third coordinate vector ( $\hat{3}$ ) shall be perpendicular to the plan containing the first and second vectors and had length equals to the distance to the surface. The lengths of the distances between the center-of-mass and the actual surfaces along the coordinate vectors are denoted $L_{1}$ , $L_{2}$ , and $L_{3}$ , in their respective order. Due to scattering, the coordinates are transformed to another virtual coordinate system that is suitable to the situation and the shape on hand. The distances travelled by the neutron along the new virtual coordinate system in the body are the neutron-chord lengths, denoted $\ell_{1}$ , $\ell_{2}$ , and $\ell_{3}$ , which need to be determined using transport equations. The average neutron-chord length is taken, in the present work, as the harmonic mean of these three distances; i.e.

\frac{1}{\bar{\ell}}=\frac{1}{3}\left(\frac{1}{\ell_{1}}+\frac{1}{\ell_{2}}+\frac{1}{\ell_{3}}\right).

(21)

Here, $\bar{\ell}$ is the average neutron-chord length.

Due to symmetry operations in the diffusion equation, we shall take the condition that if there were a symmetry making one or more of these virtual lengths undeterminable, it should take an infinity value. The behavior of the absorption and scattering transport are strongly dependent on each other and depend on the neutron energy [52]. However, a single group that cover the entire energy range of neutron (both thermal and epi-thermal) can be used to obtain meaningful value of the average neutron-chord length. The time dependent diffusion equation comprises;

\frac{1}{{\mathrm{v}_{av}}}\frac{\partial\varphi(\vec{r},t)}{\partial t}=\nabla\cdot J(\vec{r},t)-{{\Sigma}_{a}}\varphi(\vec{r},t)+Q(\vec{r},t),

(22)

where $\varphi(\vec{r},t)$ is the flux, $J(\vec{r},t)=$ $D(\vec{r},t)$ $\nabla\varphi(\vec{r},t)$ is the neutron current, $Q(\vec{r},t)$ is the neutron production rate within the medium in units of [n cm^-3s^-1]. The sample, however, is embedded, presumably, within a uniform neutron field in which the flux outside it, denoted $\varphi_{o}$ , is isotropic, uniform and does not depend on the diffusion within the sample. Considering our situation of sample absorbing neutrons, the solution of the problem comes as difference-problem in the steady-state where $Q(\vec{r},t)$ equated to $\varphi_{o}$ , and considering only the difference within the sample;

\varphi(\vec{r},t)=\varphi_{o}-\phi(\vec{r},t),

(23)

Under the condition of steady-state, the time derivative vanishes; while Eq. 22 is reduced to:

	$\displaystyle\nabla\cdot D(\vec{r})\nabla\phi(\vec{r})\text{-}{{\Sigma}_{a}}\phi(\vec{r},t)$	$\displaystyle=0$	inside the sample,		(24)
	$\displaystyle\phi(\vec{r},t)$	$\displaystyle=0$	outside the sample.		(25)

In the homogenous isotropic medium, $D(\vec{r})$ and $\Sigma_{a}(\vec{r})$ are constants, so Eq. 24 can be rewritten as:.

\nabla^{2}\phi(\vec{r})-B^{2}\phi(\vec{r})=0,

(26)

The $B^{2}$ factor depends on body materials usually called the geometric buckling factor in units of [cm^-2s^-1], see Appendix C. In steady-state, and in the present work, we rewrite the factor $B^{2}$ as reciprocal squared dimensions multiplied by (1 s^-1);

B^{2}=\sum_{i=1}^{3}B_{i}=\sum_{i=1}^{3}\frac{{\pi}^{2}}{\ell_{i}^{2}}(1\,s^{-1}),

(27)

where the values of $\ell_{i}$ are the virtual dimensions of the sample in the transformed coordinates.

The diffusion length is the mean square distance that a neutron travels in the one direction from the source to its absorption point. The steady-state condition requires that the neutron current through the body’s surfaces at its measurable dimensions, denoted $L_{1}$ , $L_{2}$ up to $L_{3}$ , be constant, i.e. the change in the leakage of current ( $\left.\nabla\cdot J(\vec{r})\right|_{\vec{r}_{B}}$ ) at these boundaries, $\vec{r}_{B}$ , be zero. For convex surface, a neutron leaving the region through the surface cannot intersect the surface again. Consequently, general solution at the surface shall satisfy the conditions;

	$\displaystyle\left.\nabla\cdot D(\vec{r})\nabla\phi(\vec{r})\right\|_{\vec{r}_{B}}$	$\displaystyle\equiv$	$\displaystyle\left.\nabla\cdot J(\vec{r})\right\|_{\vec{r}_{B}}=0,$
	$\displaystyle J(\vec{r}>\vec{r}_{B})$	$\displaystyle=$	$\displaystyle J(\vec{r}\leq\vec{r}_{B}).$		(28)

Which refer to the continuity of flux at the surface and constancy of current, respectively. Under this condition, the boundary in the difference-problem of Eqs. 24 and 25 is considered a vacuum boundary [61]. At the specific boundary vector $\vec{r}_{B}$ , the flux become $\varphi(\vec{r}_{B})=\varphi_{o}$ , or, according to Eq. 23,

\phi(\vec{r}_{B})=0

(29)

The solution of Eq. 26, as derived in Appendix C.1 for cartesian coordinates, resembles;

\phi(\underline{x},\underline{y},\underline{z})=\mathrm{Const.}e^{-iB_{1}\,\underline{x}-iB_{2}\,\underline{y}-iB_{3}\,\underline{z}}

(30)

Note that the coordinate variables are underlined from now one in order to avoid confusion among symbols. In order to satisfy the condition Eq. 29, the values of $\ell_{i}$ should be related to the measurable coordinate dimensions as follows: $\ell_{1}=2L_{1}$ , $\ell_{2}=2L_{2}$ and , $\ell_{3}=2L_{3}$ . As presented in Fig. 1, $L_{1}=W/2$ , similarly all other dimensions. i.e.

\bar{\ell}=3\left(\frac{1}{W}+\frac{1}{D}+\frac{1}{H}\right)^{-1}

(31)

Note that for infinite foils and sheets, there is only one measurable length exits, the thickness with $L_{1}$ $=t/2$ , $L_{2}=\infty$ , and $L_{3}=\infty$ making $\ell_{1}$ $=t$ , $\ell_{2}=\infty$ , $\ell_{3}=\infty$ , and

\bar{\ell}=3t

(32)

For cylindrical geometries the solution of Eq. 26, see Appendix C.2, becomes;

$\displaystyle\phi(\underline{\rho},\underline{\phi},\underline{z})=$	$\displaystyle\sum_{m=0}^{\infty}\sum_{n=0}^{\infty}C_{mn}J_{m}(\sqrt{n^{2}+B_{1}^{2}}\underline{\rho})\times$
	$\displaystyle\cos{m\underline{\phi}}e^{-\sqrt{n^{2}-B_{2}^{2}}\underline{z}}\text{, if }\,n^{2}>B_{2}^{2},$
$\displaystyle=$	$\displaystyle\sum_{m=0}^{\infty}\sum_{n=0}^{\infty}C_{mn}J_{m}(\sqrt{n^{2}+B_{1}^{2}}\underline{\rho})\times$	(33)
	$\displaystyle\cos{m\underline{\phi}}\cos{\sqrt{B_{2}^{2}-n^{2}}\underline{z}}\text{, if }\,n^{2}\leq B_{2}^{2}$	(33)

where $J_{m}$ is Bessel functions of the first kind, $C_{mn}$ is a constant and m and n are integers. Generally, m=0 and n=0 have the largest contribution. Hence

\phi(\underline{\rho},\underline{\phi},\underline{z})\simeq C_{00}J_{0}(B_{1}\,\underline{\rho})\cos{B_{2}\,\underline{\phi}}

(34)

The first root of $J_{0}$ in Eq. 34 is when the argument $B_{1}\underline{\rho}=$ 2.4048 while the root of the cosine function is when its argument $B_{2}\,\underline{\phi}$ = $\frac{\pi}{2}$ . In order to satisfy the condition Eq. 29;

	$\displaystyle\ell_{1}$	$\displaystyle=\frac{\pi}{2.4048}L_{1}$	$\displaystyle=\frac{\pi}{2.4048}R$		(35)
	$\displaystyle\ell_{2}$	$\displaystyle=2L_{2}$	$\displaystyle=H,$		(36)

and none-existence of $\underline{z}$ dependance Eq. 34 gives $\ell_{3}=\infty$ . i.e.

\bar{\ell}=\frac{3\pi RH}{2.4048H+\pi R}

(37)

Note that for finite disk of the radius $R$ and the thickness $t$ .

\bar{\ell}=\frac{3\pi Rt}{2.4048t+\pi R}

(38)

For infinite wire and cylinders, there is only one measurable coordinate length, the radius $R$

\bar{\ell}=\frac{3\pi R}{2.4048}=3.9191R

(39)

The solution of Eq. 26, as derived in the Appendix C.3, for spherical shape resembles;

\phi(\underline{r},\underline{\theta},\underline{\phi})=\sum_{l=0}^{\infty}\sum_{m=0}^{\infty}C_{mn}j_{l}(B_{1}\underline{r})P_{l}^{m}(\cos\underline{\theta})\cos{m\underline{\phi}},

(40)

where, $P_{l}^{m}(\cos\underline{\theta})$ is the associated Legendre polynomial, $j_{l}(B\,\underline{r})$ is the spherical Bessel function, and $C_{mn}$ is the integration constant. Again $C_{00}$ has largest contribution. Hence;

\phi(\underline{r},\underline{\theta},\underline{\phi})\simeq C_{00}\;j_{0}(B_{1}\,\underline{r})\;P_{0}^{0}(\cos\underline{\theta})

(41)

The condition in Eq. 29 is satisfied if the argument of the spherical Bessel function ( $B_{1}\,\underline{r}$ ) equal the root of the spherical Bessel function at $\pi$ ; i.e. the condition is satisfied if $\ell_{1}=L_{1}=R$ . Here, $P_{0}^{0}(\cos\underline{\theta})$ =1 and due to symmetry of the body, there exists none- $\underline{\theta}$ dependance which reveal $\ell_{2}=\infty$ , while none- $\underline{\phi}$ dependance requires $\ell_{3}=\infty$ ; i.e.

\bar{\ell}=3R

(42)

For any other irregular shape the average neutron-chord length can be calculated in the same manner using box geometry as approximation. Table 2 showed a comparison between our simple procedure and others.

Table 2: The geometry factor

\Omega(\bar{\ell},\Sigma_{t})

deduced from simple convex geometries presented in Fig. 1 in comparison to reported factors from different researches. Reported factors are from Refs. Ref. [17, 18, 20, 21].

Body		The geometry factor $\Omega(\bar{\ell},\Sigma_{t})$
Shape	Dimensions	Present work	Reported
Sheet-infinite	$t$	${3t}/{P_{0}}$	$1.5t$
Disc	$t$ , $R$	$\frac{3\pi tR}{(P_{0}(2.4048t+\pi R))}$	–
Box/Slap	$H$ , $W$ , $D$	$\frac{3HDW}{(P_{0}(HW+HD+DW))}$	–
Sphere	$R$	${3R}/{P_{0}}$	$R$
Cylinder	$H$ , $R$	$\frac{3\pi RH}{(P_{0}(2.4048H+\pi R))}$	$1.65\frac{HR}{H+R}$
Cylinder-infinite	$R$	${3.9191R}/{P_{0}}$	–
General shape	$L_{1}$ , $L_{2}$ , $L_{3}$	$\frac{3}{P_{0}}\left(\frac{1}{\ell_{1}}+\frac{1}{\ell_{2}}+\frac{1}{\ell_{3}}\right)^{-1}$	–

The average neutron-chord length is larger than the dimensions of the body due to the irregular path of neutrons in the body’s material. For a sphere with radius R is $\bar{\ell}=$ 3R, not the value of 2R, while for infinite foil it is also three times its thickness, not the value of 1.5t, due to the average of the cosine in the neutron scattering path length inside the volume. However, when Sjostrand et al. [62] calculated the average neutron-chord length for a sphere assuming an isotropic flux distribution, the result was equal to the radius R due to use of different weighting factors.

III.3 Probability of the neutron interaction

The next step is to obtain a mathematical formula for the probability of single interaction within the volume ( $P_{0}$ ), i.e. the probability that a neutron will suffer at least one more interaction. In the case of thermal energies, the domain of Maxwellian distribution below cadmium cut-off energy may be considered as that of the averaged energy at 0.025 eV for which scattering and multiple scattering shall not disturbs the overall neutron energy distribution [63]. The neutron absorptions are the result of the various neutron resonances which are predominant in the epi-thermal region including those of capture and possible fission components while scattering components cause escape of neutrons from this region.

The neutron-escape probability, the factor $p$ in reactor physics, measures the fraction of neutrons that have escaped absorption and still exist after having been slowed-down from their epi-thermal energies to – say thermal energies – due to these “resonance traps” and reduces the absorption losses [64, 65]. Several authors had tried calculate the resonance escape probabilities from first principles [66, 67, 68, 65] while other calculate directly from thermal utilization factor of reactors (cf. Ref. [69]).

In the present work and under the condition in Eq. 14, the probability of interaction is obtained from the attenuation relation ( $\psi=$ $\exp$ $(-\Sigma_{t}$ $\times$ $\mathrm{mean~{}distance})$ ). The scattered neutron continues to exist within the body. The mean travelled distance is the averaged neutron-chord length which allows us to write directly and according to Rothenstein [65], and for approximation, the following

P_{0}\sim(1-\exp(-\Sigma_{t}\bar{\ell})),

(43)

which satisfies the condition in Eq. 14. In the extended range of epi-thermal neutrons, the flux varies as 1/ $E_{n}$ from cadmium cut-off energy at 0.5 eV to the end of the neutron spectrum – say 1 MeV. Multiple scattering disturbs the energy distribution by reducing the number of neutrons in the epi-thermal region. In reactor physics, this phenomenon is described by resonance escape probability, which is the probability that a neutron will slow down from fission energy to thermal energies without being captured by a nuclear resonance. This phenomenon depends on the diffusion properties of the medium. In the comparison given in Fig. 3, the values of $G_{epi}$ vary as in Eq. 43 but with $\bar{\ell}^{2/3}$ replaces $\bar{\ell}$ . There was no clear reason about this dependence. However, there is a simple experimental remark in neutron physics: whenever the neutron energy distribution repudiate the proper thermalization distribution (Maxwellian+1/E dependance in epi-thermal region) by any mean such as absorption, the neutrons rapidly redistribute its velocity population within the diffusion distance to follow the proper distribution – up to thermalization. [70, 71, 2]. To compensate such dependence, in the present work we had introduced a parameterized factor to enhance formula in Eq. 43 in the epi-thermal range as follows;

P_{0}\sim\frac{1}{p_{\mathrm{escape}}}\left(1-\exp\left(-\Sigma_{t}\bar{\ell}{p_{\mathrm{escape}}}\right)\right),

(44)

where

{p_{\mathrm{escape}}}=2\frac{\sqrt[2]{|\Sigma_{a}-\Sigma_{s}|}}{\sqrt[2]{\Sigma_{a}}\sqrt[3]{\Sigma_{s}\bar{\ell}}}.

(45)

which, also, satisfies the condition in Eq. 14.

The subscript $(\mathrm{energy~{}domain})$ is to replaced by “ $th$ ” in case of thermal neutron energies below cadmium cutoff energy ( $\sim$ 0.5 eV) or by “ $epi$ ” for epi-thermal neutrons having energy domain greater than the cadmium cutoff energy. Here, we have used two notions of flux $\varphi_{o}$ and $\varphi$ which stand for unperturbed neutron flux for which the material is diluted or absent [14] and the measured self-shielded neutron flux in the vicinity of the material, respectively.

III.4 Verification with experiment

The obtained mathematical values with present ab initio model were compared with experimental values in Fig. 2. The exact parameters of the experimental data such as foil thickness and wire radius, cylinder hight were obtained from the original sources (whether these were literature or our previous experiments). In the thermal energy range, the derived formula in Eqs. 15 and 43 gave a good representation of the experimental data for In, Au and Co within the experimental uncertainty, whether it was wires or foils .

In the epi-thermal region, the integral cross-section of Eq. 18 is replaced by the resonance integral. Table 1 contains the element-averaged resonance integrals for 1/E averaged neutron distribution together with the thermal cross- section data based on Maxwellian distribution of neutron energies for capture reactions [10] and scattering reactions [9]. These integral data were used to calculate epi-thermal self-shielding factor, G_epi, and represented by lines in Fig. 3. Experimental results from literature of Gonalves et al. [17], Lopes [36, 37], McGarry [38], Brose [39], Yamamoto et al. [40], Jefferies et al. [41], Eastwood $\&$ Werner [35], and Kumpf [42] were used for comparison. Elements were Au, Co, Mn and Na in a form of wire, foils or infinite slabs. Note that our model calculations were based on the derived formula in Eq. 15 and the interaction probability of Eq. 44. Parameters of such as foil thickness and wire radius, cylinder hight were obtained from the original sources of the experimental data. With the adaptation in the Eq. 45, our model gave a good representation of the experimental data for Au, Co, Mn and Na within the experimental uncertainty. Otherwise, the model curve shall be more steeper and messes the experimental data. As shown in Fig. 3, there is a slightly different between experimental, calculated values and our model in gold and Indium wires and foils.

It is clear that the neutrons self-shielding factor depends not only on the properties and geometry of the material but also on the neutron energy range as shown in Fig. 2 and 3. A comparison between the present approach of ab initio calculations, considering the extreme cases of infinite wire and infinite foil, and the empirical equation given in Refs. [17, 18, 20, 21, 19] are given in Appendix D of the present work.

IV Conclusion

In the vicinity of neutron-absorbing elements within the sample, neutron flux shall be modified continuously with the depth. The activation formulae that take the flux as a constant value shall be corrected by a self-shielding factor. The self-shielding corrected neutron flux factor are often obtained from numerous approaches, both empirical based on fitting or analytic analysis as presented within the present work. Understanding the physics behind self-shielding enabled the extension of a simple thermal neutron picture into the epithermal energies with the possibility for application to high-energy neutrons. Equation 15 together with its descriptive parameters in Eqs. 17, 19, 20, 43, and 44 were satisfactory in the determination of the self-shielding factors when the average chord lengths are calculated from our derived formulae in Table 2. The analytical formulae enable its implantation the longer-term application in the analysis of neutron activation and neutron-induced effects of materials for different materials in different geometries, especially neutron shields, using integral parameter representation, instead of spectroscopic one.

Appendix A Advantage compared to Monté Carlo methods

The use of Monté Carlo simulation software (MC) for calculating the self-shielding factors is feasible but not yet efficacious. The principle underlying MC is to avoid the direct analytical solution of the problem. The goal of MC is to simulate and average a sufficiently large number of particle histories to obtain estimates of the flux which include rigorous approximations. According to Larson, MC of difficult problems are often very costly to set up and run. To make the MC code run with acceptable efficiency, the code users must specify a large number of biasing parameters, which are specialized to each different problem. Determining these parameters can be difficult and time-consuming. Also, even when the biasing parameters are well-chosen, MC converges slowly and non-monotonically with increasing run time. Thus, while MC solutions are free of truncation errors, they are certainly not free of statistical errors, and it is challenging to obtain MC solutions with sufficiently small statistical errors, and with acceptable cost. Finally, the non-analog techniques that have been developed for making MC simulations acceptably efficient and were useful for source-detector problems – in which a detector response in a small portion of phase space is desired –are not useful for obtaining efficient global solutions, over all of phase space. Generally, MC solutions work best when very limited information about the flux (e.g. a single detector response) is desired in a given simulation. MC is feasible for calculating self-shielding of a single sample, it requires time and an experienced user to make the acceptable accuracy and efficiency, the cost most scientists cannot afford to just calculate a single parameter in their routine work. For example, considering set of different samples need to be analyzed using neutron activation for the purpose of elemental analysis, MC SIM needs experience and a lot of time to reduce fluctuation, adopt the geometry, consideration of the neutron transport inside and outside the sample, and the absorption and scattering phenomenon as a function of neutron energy.

Although our intention while deriving and validating the present mathematical formulae was focused on avoiding such costs and enhancing present existing empirical equations and transforming all the problems from the spectroscopic set of parameters, such as thermal cross-section, thickness, width, height, radius, shape, a width of neutron first resonance, the width of the first gamma resonance, etc. into an integrated set of well-known parameters, thermal cross-section, resonance integral, average chord length. All remaining factors are calculated from these three parameters. Of course, the thermal cross-section and sample mass and composition are common.

Appendix B Contribution of velocity distribution

According to the results of Blaauw [16], calculation of the reaction rate is need to be with neutron density averaged macroscopic cross-section (function of velocity) instead of flux averaged macroscopic cross-section (energy dependent). Blaauw found that the self-shielding factors calculated for mono energetic neutrons yields the same results as if they are used with the flux averaged macroscopic cross-section provided that the neutron density averaged macroscopic cross-section given by

<\Sigma>=\frac{2}{\sqrt{\pi}}\sqrt{\frac{T_{o}}{T}}\Sigma_{o},

(46)

is used instead of the flux averaged capture cross-section given by

<\Sigma>=\frac{\sqrt{\pi}}{2}\sqrt{\frac{T_{o}}{T}}\Sigma_{o}

(47)

Blaauw [16] results showed that the volume-averaged attenuation self-shielding factor in extended neutron distributions, has an extra term that depends on the statistical moments of deviation in reciprocal velocity average. The contribution of this extra factor had been estimated by Goncalves et al. [21] to be around 6 $\pm$ 1%.

The higher order terms in Eq. 9 (Denoted $\mathcal{R}$ ) can be obtained from Blaauw [16] and adopted to our notions as;

\mathcal{R}\cong\sum_{i=1}^{\infty}\frac{(-1)^{i}}{i!}(\Sigma_{a})^{i}v_{o}^{i}\left(\left\langle\frac{1}{v^{i}}\right\rangle-\left\langle\frac{1}{v}\right\rangle^{i}\right)\frac{\int_{V}(\vec{r})^{i}d\vec{r}}{V}

(48)

For the first term, $i=1$ , the value of $\left(\left\langle\frac{1}{v}\right\rangle-\left\langle\frac{1}{v}\right\rangle\right)$ vanishes. While for approximate spherical symmetry the integral yields the average squared length over volume of sphere, The first term comprises;

	$\displaystyle\int_{V}(\vec{r})^{2}d\vec{r}$	$\displaystyle=$	$\displaystyle\int_{0}^{\bar{\ell}}\int_{0}^{2\pi}\int_{0}^{\pi}\underline{r}^{2}(\underline{r}\sin\underline{\phi})d\underline{\theta}d\underline{\phi}\underline{r}d\underline{\theta}dr$		(49)
		$\displaystyle=$	$\displaystyle\frac{4}{5}\pi\bar{\ell}^{5}=\frac{4}{5}\frac{3}{4}\frac{4}{3}\pi\bar{\ell}^{3}\bar{\ell}^{2}=\frac{3}{5}V\bar{\ell}^{2}$		(50)

For the first approximation, only term of i=2 has an effective contribution. Hence,

\mathcal{R}\cong\frac{+1}{2}(\Sigma_{a})^{2}v_{o}^{2}\left(\left\langle\frac{1}{v^{2}}\right\rangle-\left\langle\frac{1}{v}\right\rangle^{2}\right)\frac{3}{5}\bar{\ell}^{2}

(51)

For Maxwellian velocity distribution

	$\displaystyle\left\langle\frac{1}{v}\right\rangle$	$\displaystyle=$	$\displaystyle\frac{2}{\sqrt{\pi}}\frac{1}{v_{o}}$		(52)
	$\displaystyle\left\langle\frac{1}{v^{2}}\right\rangle$	$\displaystyle=$	$\displaystyle 2\frac{1}{v_{o}^{2}}$		(53)

$\displaystyle\mathcal{R}$	$\displaystyle\cong(\Sigma_{a}\bar{\ell})^{2}\frac{3}{5}\left(1-\frac{2}{\pi}\right)\approx 0.218(\Sigma_{a}\bar{\ell})^{2},\,\,\,Maxwellian$	(54)
	$\displaystyle\cong 0,\,\,\,monoenergetic$
	$\displaystyle\cong 0,\,\,\,20v_{o}\leq v\leq 4\times 10^{7}v_{o}$

The first formulae is valid only for the entire range of Maxwellian neutron distribution. However, for practical use, only the epi-thermal range between cadmium cutoff at 0.5 eV upto about 1 MeV is used. In such cases, the difference $\left(\left\langle\frac{1}{v^{2}}\right\rangle-\left\langle\frac{1}{v}\right\rangle^{2}\right)$ practically vanishes, $\mathcal{R}\cong 0$ .

There is an additional reason why this value is ignored within the sample in the present work. The Blaauw [16] derivation is based on the idea that the neutron flux distribution have constant shape as it passes through the depth $\bar{\ell}$ . However, there is a simple experimental remark in neutron physics: whenever the neutron energy distribution repudiate the proper thermalization distribution (Maxwellian+1/E dependance in epi-thermal region) by any mean such as absorption, the neutrons rapidly redistribute its velocity population within the diffusion distance to follow the proper distribution – up to thermalization. The idea is, even if there is absorption of neutrons having a velocity $v$ in Eq. 48 at some distance, there were a sort of recovery of that distribution. And hence, the difference in Eq. 51 has much less value than expected by Blaauw.

Appendix C Determination of chord lengths based on neutron transport formulae

The time dependent diffusion equation comprises;

\frac{1}{{\mathrm{v}_{av}}}\frac{\partial\varphi(\vec{r},t)}{\partial t}=\nabla\cdot J(\vec{r},t)-{{\Sigma}_{a}}\varphi(\vec{r},t)+Q(\vec{r},t),

(55)

where $J(\vec{r},t)=$ $D(\vec{r},t)$ $\nabla\varphi(\vec{r},t)$ is the neutron current, $Q(\vec{r},t)$ is the neutron production rate within the medium in units of n cm^-3s^-1. Under the condition of steady-state ( $\frac{\partial\varphi(\vec{r},t)}{\partial t}=0$ ) and Considering:

•

the sample is embedded within a uniform neutron field in which the flux outside it, $\varphi_{o}$ , to be isotropic, uniform and does not depend on the diffusion within the sample,
•

our situation of sample absorbing neutrons not generating it,

the solution of the problem comes as difference-problem in the steady-state where $Q(\vec{r},t)$ equated to $\varphi_{o}$ , and considering only the difference replacing $\varphi(\vec{r})$ by $\varphi_{o}-phi(\vec{r})$ : Then

\nabla\cdot D(\vec{r})\nabla\phi(\vec{r})-\Sigma_{a}(\vec{r})\phi(\vec{r})=0

(56)

In the homogenous isotropic medium, $D(\vec{r})$ and $\Sigma_{a}(\vec{r})$ are constants.

\nabla^{2}\phi(\vec{r})-B^{2}\phi(\vec{r})=0,

(57)

The $B$ factor is the geometric buckling factor in reactor physics. Taking into consideration that $\Sigma_{a}=$ $1/\lambda_{a}$ , and $D=\lambda_{tr}/3$ where $\lambda_{a}$ is the absorption mean-free path and $\lambda_{tr}$ is the transport scattering diffusion length given by the more advanced transport theory in terms of transport and absorption cross-sections equation as [72, 73];

\lambda_{tr}=\frac{1}{\Sigma_{a}+\Sigma_{s}(1-\bar{\mu})},

(58)

where $\bar{\mu}=\frac{2}{3A}$ is average value of the cosine of the angle in the lab system. So,

B^{2}=\frac{\Sigma_{a}[cm^{-1}]}{D[cm\,s^{-1}]},

(59)

C.1 Rectangular geometries

In cartesian coordinates, Eq. 57 is reduced to three independent equations by separation of variables assuming $\phi(\underline{x},\underline{y},\underline{z})=$ $\phi_{x}(\underline{x})$ $\phi_{y}(\underline{y})$ $\phi_{z}(\underline{z})$ . I.e.

$\displaystyle\left(\frac{d^{2}}{d\underline{x}^{2}}-B_{1}^{2}\right)\phi_{x}(\underline{x})$	$\displaystyle=$	$\displaystyle 0$	(60)
$\displaystyle\left(\frac{d^{2}}{d\underline{y}^{2}}-B_{2}^{2}\right)\phi_{y}(\underline{y})$	$\displaystyle=$	$\displaystyle 0$	(61)
$\displaystyle\left(\frac{d^{2}}{d\underline{z}^{2}}-B_{3}^{2}\right)\phi_{z}(\underline{z})$	$\displaystyle=$	$\displaystyle 0,$	(62)

The general solution is:

\phi(\underline{x},\underline{y},\underline{z})=\mathrm{Const.}e^{-iB_{1}\,\underline{x}}e^{-iB_{2}\,\underline{y}}e^{-iB_{3}\,\underline{z}}

(63)

C.2 Cylindrical geometries

For a definite convex shapes of having cylindrical geometries, Eq. 57 becomes the Helmholtz differential equation;

\displaystyle\left(\frac{1}{\underline{\rho}}\frac{\partial}{\partial\underline{\rho}}\left(\underline{\rho}\frac{\partial}{\partial\underline{\rho}}\right)\text{+}\frac{1}{\underline{\rho}^{2}}\frac{\partial^{2}}{\partial\underline{\phi}^{2}}\text{+}\frac{\partial^{2}}{\partial\underline{z}^{2}}\text{+}B^{2}\right)\phi(\underline{\rho},\underline{\phi},\underline{z})=0.

(64)

In which the metric tensor scale factors are 1, $\rho$ , and 1 for the coordinates $\rho$ , $\phi$ , $z$ , respectively. Separation of variables is done by writing $\phi(\underline{\rho},\underline{\phi},\underline{z})=\phi_{\rho}(\underline{\rho})\phi_{\phi}(\underline{\phi})\phi_{z}(\underline{z})$ . Eq. 64 becomes;

	$\displaystyle\left({\frac{\underline{\rho}^{2}}{\phi_{\rho}}\frac{d^{2}\phi_{\rho}}{d\underline{\rho}^{2}}+\frac{\underline{\rho}}{\phi_{\rho}}\frac{d\phi_{\rho}}{d\underline{\rho}}}\right)+\frac{1}{\phi_{\phi}}\frac{d^{2}\phi_{\phi}}{d\underline{\phi}^{2}}+\frac{\underline{\rho}^{2}}{\phi_{z}}\frac{d^{2}\phi_{z}}{d\underline{z}^{2}}-$
	$\displaystyle\underline{\rho}^{2}B^{2}=0.$		(65)

Solution requires negative separation constants – say $m^{2}$ in order to maintain the periodicity in $\varphi$ ; hence,

\frac{1}{\phi_{\phi}}\frac{d^{2}\phi_{\phi}}{d\underline{\phi}^{2}}=-m^{2},

(66)

which has a general solution

\phi(\underline{\phi})=\mathrm{C}_{1}\cos{m\,\underline{\phi}}+\mathrm{C}_{2}\sin{m\,\underline{\phi}}.

(67)

where $C_{1}$ , and $C_{2}$ are constants. Hence,

\left({\frac{1}{\phi_{\rho}}\frac{d^{2}\phi_{\rho}}{d\underline{\rho}^{2}}+\frac{1}{\underline{\rho}\phi_{\rho}}\frac{d\phi_{\rho}}{d\underline{\rho}}}\right)-\frac{m^{2}}{\underline{\rho}^{2}}+B^{2}+\frac{1}{\phi_{z}}\frac{d^{2}\phi_{z}}{d\underline{z}^{2}}=0,

(68)

i.e. for finite cylinder, there are two coordinate lengths, radius and hight

	$\displaystyle\left({\frac{1}{\phi_{\rho}}\frac{d^{2}\phi_{\rho}}{d\underline{\rho}^{2}}+\frac{1}{\underline{\rho}\phi_{\rho}}\frac{d\phi_{\rho}}{d\underline{\rho}}}\right)-\frac{m^{2}}{\underline{\rho}^{2}}+$
	$\displaystyle B_{1}^{2}+B_{2}^{2}+\frac{1}{\phi_{z}}\frac{d^{2}\phi_{z}}{d\underline{z}^{2}}=0$		(69)

The $\phi_{z}$ must not be sinusoidal at $\pm\infty$ which lead to positive separation constant – say $n^{2}$ ;

\frac{1}{\phi_{z}}\frac{d^{2}\phi_{z}}{d\underline{z}^{2}}=n^{2}-B_{2}^{2},

(70)

\frac{d^{2}\phi_{\rho}}{d\underline{\rho}^{2}}+\frac{1}{\underline{\rho}}\frac{d\phi_{\rho}}{d\underline{\rho}}+\left({n^{2}+B_{1}^{2}-\frac{m^{2}}{\underline{\rho}^{2}}}\right)\phi_{\rho}=0.

(71)

The solutions are

	$\displaystyle\phi_{Z}(\underline{z})=\mathrm{C}_{3}e^{-\sqrt{n^{2}-B_{2}^{2}}\underline{z}}+\mathrm{C}_{4}e^{\sqrt{n^{2}-B_{2}^{2}}\underline{z}}$
	$\displaystyle\,if\,n^{2}>B_{2}^{2}$		(72)
	$\displaystyle=\mathrm{C}_{5}\cos{\sqrt{B_{2}^{2}-n^{2}}\underline{z}}+\mathrm{C}_{6}\sin{\sqrt{B_{2}^{2}-n^{2}}\underline{z}}$
	$\displaystyle\,if\,n^{2}\leq B_{2}^{2}$		(73)

\displaystyle\phi_{\rho}(\underline{\rho})=\mathrm{C}_{7}J_{m}(\sqrt{n^{2}\text{+}B_{1}^{2}}\underline{\rho})\text{+}\mathrm{C}_{8}Y_{m}(\sqrt{n^{2}\text{+}B_{1}^{2}}\underline{\rho}),

(74)

where $J_{m}$ and $Y_{m}$ are Bessel functions of the first kind and second Kind, respectively. These results requires that $n$ and $m$ be integers. $Y_{n}(0)$ = $-\infty$ which lead to un-physical solution, hence, $\mathrm{C}_{2}$ =0 and $\mathrm{C}_{8}$ =0. Similarly, $\mathrm{C}_{4}$ =0 and $\mathrm{C}_{6}$ . The solution is reduced to;

	$\displaystyle{\phi}(\underline{\rho},\underline{\phi},\underline{z})=\sum_{m=0}^{\infty}\sum_{n=0}^{\infty}C_{mn}J_{m}(\sqrt{n^{2}+B_{1}^{2}}\underline{\rho})\cos{m\,\underline{\phi}}$
	$\displaystyle e^{-\sqrt{n^{2}-B_{2}^{2}}\underline{z}}\,if\,n^{2}>B_{2}^{2}$		(75)
	$\displaystyle=\sum_{m=0}^{\infty}\sum_{n=0}^{\infty}C_{mn}J_{m}(\sqrt{n^{2}+B_{1}^{2}}\underline{\rho})\cos{m\,\underline{\phi}}$
	$\displaystyle\cos{\sqrt{B_{2}^{2}-n^{2}}\underline{z}}\,if\,n^{2}\leq B_{2}^{2}$		(76)

where $C_{mn}$ is a constant that depends on the values of m and n.

C.3 Spherical geometries

In spherical coordinates, Eq. 57 resembles;

	$\displaystyle\left(\frac{1}{\underline{r}^{2}}\frac{\partial}{\partial\underline{r}}\left(\underline{r}^{2}\frac{\partial}{\partial\underline{r}}\right)+\frac{1}{\underline{r}^{2}\sin\underline{\theta}}\frac{\partial}{\partial\underline{\theta}}\left(\sin\underline{\theta}\frac{\partial}{\partial\underline{\theta}}\right)+\right.$
	$\displaystyle\left.\frac{1}{\underline{r}^{2}\!\sin^{2}\underline{\theta}}\frac{\partial^{2}}{\partial\underline{\phi}^{2}}\right)\phi+B_{1}^{2}\phi=0.$		(77)

Because of the spherical symmetry there is only one value of $B$ = $B_{1}$ . Values of $B_{2}$ and $B_{3}$ vanishes; i.e. $\ell_{2}=\infty$ and $\ell_{3}=\infty$ .

	$\displaystyle\left(\frac{1}{\phi_{r}}\frac{d}{dr}\left(\underline{r}^{2}\frac{d\phi_{r}}{d\underline{r}}\right)+\frac{1}{\phi_{\theta}}\frac{1}{\sin\underline{\theta}}\frac{d}{d\underline{\theta}}\left(\sin\underline{\theta}\frac{d\phi_{\theta}}{d\underline{\theta}}\right)+\right.$
	$\displaystyle\left.\frac{1}{\phi_{\phi}}\frac{1}{\sin^{2}\underline{\theta}}\frac{d^{2}\phi_{\phi}}{d\underline{\phi}^{2}}\right)+\underline{r}^{2}B_{1}^{2}=0.$		(78)

Separation of variable requires substitution of $\phi(\underline{r},\underline{\theta},\underline{\phi})$ by $\phi_{r}(\underline{r})\phi_{\theta}(\underline{\theta})\phi_{\phi}(\underline{\phi})$ . Separating the $r$ term with separation constant $l(l+1)$

\frac{1}{\phi_{r}}\frac{d}{d\underline{r}}\left(\underline{r}^{2}\frac{d\phi_{r}}{d\underline{r}}\right)-l(l+1)+\underline{r}^{2}B_{1}^{2}=0

(79)

	$\displaystyle\frac{1}{\phi_{\theta}}\frac{1}{\sin\underline{\theta}}\frac{d}{d\underline{\theta}}\left(\sin\underline{\theta}\frac{d\phi_{\theta}}{d\underline{\theta}}\right)+\frac{1}{\phi_{\varphi}}\frac{1}{\sin^{2}\underline{\theta}}\frac{d^{2}\phi_{\varphi}}{d\underline{\varphi}^{2}}$
	$\displaystyle+l(l+1)=0$		(80)

Eq. C.3 is separated by separation constant $m^{2}$ , then

\displaystyle\frac{1}{\phi_{\theta}}\frac{1}{\sin\underline{\theta}}\frac{d}{d\underline{\theta}}\left(\sin\underline{\theta}\frac{d\phi_{\theta}}{d\underline{\theta}}\right)+l(l+1)\sin\underline{\theta}-m^{2}=0

(81)

\frac{1}{\phi_{\phi}}\frac{d^{2}\phi_{\phi}}{d\underline{\phi}^{2}}+m^{2}=0

(82)

solution of Eqs. 79, 81, and 82 yield the following solution;

\phi(\underline{r},\underline{\theta},\underline{\phi})=\sum_{l=0}^{\infty}\sum_{m=0}^{\infty}C_{mn}\;j_{l}(B_{1}\,\underline{r})\;P_{l}^{m}(\cos\underline{\theta})\cos{m\,\underline{\phi}}

(83)

by ignoring the anti symmetric terms. Here, $P_{l}^{m}(\cos\underline{\theta})$ gives the associated Legendre polynomial while $j_{l}(B\,\underline{r})$ is the spherical Bessel function.

Appendix D Comparison with empirical formula

Figures 4 and 5 represent the thermal self-shielding factor calculated using empirical equations given by Refs. [17, 18, 20, 21, 19].

Their curves were calculated with the specific values of cross-section given in their manuscripts, which does not equal to the recommended cross-sections in literatures. Our approach was calculated using the extreme approximation of infinite foil (having only one variable with is the thickness) and with the general cross-section values in Refs. [53, 54, 10]. Based on this comparison based on extreme cases of infinite wire and infinite foil, the results of our approach matched the empirical equation in most cases, where it was already succeeded. Note that: there is no such infinite foil or infinite wire in experimental situations.

Conflict of interest

The authors declare that they have no known source for conflict of interest with any person.

References

Tohamy et al. [2019] M. Tohamy, E. K. Elmaghraby, and M. Comsan, Nucl. Instrum. Meth. Phys. Res. A 942, 162387 (2019).
Tohamy et al. [2021] M. Tohamy, E. K. Elmaghraby, and M. N. H. Comsan, Phys. Scr. 96, 045304 (2021).
Tohamy et al. [2020] M. Tohamy, E. K. Elmaghraby, and M. Comsan, Appl. Radiat. Isotopes 165, 109340 (2020).
Nakamura et al. [2003] S. Nakamura, H. Wada, O. Shcherbakov, K. Furutaka, H. Harada, and T. Katoh, J. Nucl. Sci. Technol. 40, 119 (2003).
Ali and Elmaghraby [2020] A. Ali and E. K. Elmaghraby, Nucl. Instrum. Meth. Phys. Res. B 471, 63 (2020).
Farina-Arboccò et al. [2012] F. Farina-Arboccò, P. Vermaercke, K. Smits, L. Sneyers, and K. Strijckmans, J. Radioanal. Nucl. Chem. 296, 931 (2012).
Elmaghraby et al. [2019a] E. K. Elmaghraby, E. Salem, Z. Yousef, and N. El-Anwar, Phys. Scr. 94, 015301 (2019a).
Chilian et al. [2010] C. Chilian, R. Chambon, and G. Kennedy, Nucl. Instrum. Meth. Phys. Res. A 622, 429 (2010).
Elmaghraby [2019a] E. K. Elmaghraby, Phys. Scr. 94, 065301 (2019a).
Elmaghraby [2017a] E. K. Elmaghraby, Euro. Phys. J. Plus 132, 249 (2017a).
Jacimovic et al. [2010] R. Jacimovic, A. Trkov, G. Zerovnik, L. Snoj, and P. Schillebeeckx, Nucl. Instrum. Meth. Phys. Res. A 622, 399 (2010), proceedings of the 5th International k0 Users Workshop.
Elmaghraby [2017b] E. K. Elmaghraby, Nucl. Instrum. Meth. Phys. Res. B 398, 42 (2017b).
Elmaghraby et al. [2019b] E. K. Elmaghraby, G. Y. Mohamed, and M. Al-abyad, Nucl. Phys. A984, 112 (2019b).
Mahmoud et al. [2022a] A. W. Mahmoud, E. K. Elmaghraby, A. H. M. Soliman, E. Salama, A. Elghazaly, and S. A. El-fiki, in Procedings of the 1^st International Conference on Pure and Applied Physics (ICPAP2021), edited by M. Abdel-Harith and et al. (2022).
Fleming [1982] R. F. Fleming, The Int. J. Appl. Radiat. Isotopes 33, 1263 (1982).
Blaauw [1995] M. Blaauw, Nucl. Instrum. Meth. Phys. Res. A 356, 403 (1995).
Gonalves et al. [2001] I. Gonalves, E. Martinho, and J. Salgado, Appl. Radiat. Isotopes 55, 447 (2001).
Martinho et al. [2003] E. Martinho, I. Gonnçalves, and J. Salgado, Appl. Radiat. Isotopes 58, 371 (2003).
Salgado et al. [2004] J. Salgado, I. F. Goncalves, and E. Martinho, Nucl. Sci. Eng. 148, 426 (2004).
Martinho et al. [2004] E. Martinho, J. Salgado, and I. Gonalves, J. Radioanal. Nucl. Chem. 261, 637 (2004).
Goncalves et al. [2004] I. Goncalves, E. Martinho, and J. Salgado, Nucl. Instrum. Meth. Phys. Res. B 213, 186 (2004).
Sudarshan et al. [2005] K. Sudarshan, R. Tripathi, A. Nair, R. Acharya, A. Reddy, and A. Goswami, Analytica Chimica Acta 549, 205 (2005).
Nasrabadi et al. [2007] M. Nasrabadi, M. Jalali, and A. Mohammadi, Nucl. Instrum. Meth. Phys. Res. B 263, 473 (2007).
Bolyatko et al. [1983] V. V. Bolyatko, M. Y. Vyrskiı, A. A. Ilyushkin, G. Manturov, V. Mashkovich, M. Nikolaev, V. Sakharov, and A. Suvorov, Error estimation in reactor self shielding calculations, edited by V. Mashkovich, AIP translated series (AIP, 1983).
Moll et al. [2020] E. Moll, C. Aplin, and D. Henderson, Ann. Nucl. Energy 136, 106990 (2020).
Larsen [2006] E. W. Larsen, in Computational Methods in Transport, edited by F. Graziani (Springer Berlin Heidelberg, Berlin, Heidelberg, 2006) pp. 513–534.
Mahmoud et al. [2022b] A. W. Mahmoud, E. K. Elmaghraby, E. Salama, A. Elghazaly, and S. A. El-fiki, Br. J. Phys. (2022b).
Taylor and Linacre [1964] N. Taylor and J. Linacre, The use of Cobalt as an Accurate Thermal Neutron Flux Monitor, Tech. Rep. (United Kingdom Atomic Energy Authority, 1964).
Carre et al. [1965] J. Carre, F. Roullier, and R. Vidal, CEA–Département des Études de Piles, Service des Expériences Critiques, Rapport MIN 76 (1965).
Hasnain et al. [1961] S. Hasnain, T. Mustafa, and T. Blosser, Thermal Neutron Density Perturbations by Foils in Water, Report ORNL-3193 (Oak Ridge National Laboratory, 1961).
SOLA [1960] A. SOLA, Nucleonics 18, 78 (1960).
Walker et al. [1963] J. V. Walker, J. D. Randall, and R. C. Stinson Jr, Nuclear Science and Engineering 15, 309 (1963).
Klema and Ritchie [1952] E. Klema and R. Ritchie, Physical Review 87, 167 (1952).
Crane and Doerner [1963] J. Crane and R. Doerner, Nucl. Sci. Eng. 16, 259 (1963).
Eastwood and Werner [1962] T. Eastwood and R. Werner, Nucl. Sci. Eng. 13, 385 (1962).
Lopes et al. [1990] M. Lopes, A. Molina, et al., Kerntechnik 55, 49 (1990).
Lopes [1991] M. Lopes, Sensitivity of self-powered neutron detectors to thermal and epithermal neutrons with multiple collision treatment, Ph.D. thesis, PhD Thesis, University of Coimbra, 1991 (in Portuguese) (1991).
McGarry [1964] E. McGarry, Transactions of the American Nuclear Society (US) 7 (1964).
Brose [1964] M. Brose, Nukleonik 6, 134 (1964).
Yamamoto and Yamamoto [1965] H. Yamamoto and K. Yamamoto, Journal of Nuclear Science and Technology 2, 421 (1965).
Jefferies et al. [1983] S. Jefferies, T. Mac Mahon, J. Williams, A. Ahmad, and T. Ryves, in Nuclear Data for Science and Technology (Springer, 1983) pp. 681–684.
Kumpf [1986] H. Kumpf, Nucl. Instrum. Meth. Phys. Res. A 251, 193 (1986).
Elmaghraby [2019b] E. K. Elmaghraby, Phys. Scr. 94, 065301 (2019b).
Stuart [1957] G. W. Stuart, Nucl. Sci. Eng. 2, 617 (1957).
Blaauw [1996] M. Blaauw, Nucl. Sci. Eng. 124, 431 (1996).
Davison and Sykes [1957] B. Davison and J. B. Sykes, Neutron transport theory, International series of monographs on physics (Clarendon Press, Oxford, 1957).
Henderson and Maynard [1989] D. L. Henderson and C. W. Maynard, Nucl. Sci. Eng. 102, 172 (1989).
Hörmander [2018] L. Hörmander, The diffusion approximation in neutron transport theory, in Unpublished Manuscripts : from 1951 to 2007 (Springer International Publishing, 2018) pp. 47–51.
Bair et al. [2020] J. Bair, P. Blaszczyk, P. Heinig, V. Kanovei, M. G. Katz, and T. McGaffey, arXiv e-prints , arXiv:2003.00438 (2020), arXiv:2003.00438 [math.HO] .
de Kruijf and Kloosterman [2003] W. de Kruijf and J. Kloosterman, Ann. Nucl. Energy 30, 549 (2003).
Trkov et al. [2009] A. Trkov, G. Zerovnik, L. Snoj, and M. Ravnik, Nucl. Instrum. Meth. Phys. Res. A 610, 553 (2009).
Mahmoud et al. [2022c] A. W. Mahmoud, E. K. Elmaghraby, A. H. M. Solieman, E. Salama, A. Elghazaly, and S. A. El-fiki, Stimulated perturbation on the neutron flux distribution in the mutually-dependent source-to-absorber geometry (2022c), arXiv:2204.13246 [hep-ph] .
Sukhoruchkin et al. [1998] S. I. Sukhoruchkin, Z. N. Soroko, and H. Schopper, Tables of Neutron Resonance Parameters, 1st ed., Landolt-Börnstein - Group I Elementary Particles, Nuclei and Atoms 16B : Elementary Particles, Nuclei and Atoms, Vol. 16C (Springer-Verlag Berlin Heidelberg, 1998).
Sukhoruchkin et al. [2009] S. I. Sukhoruchkin, Z. Soroko, F. Gunsing, V. Pronyaev, and H. Schopper, Neutron Resonance Parameters, 1st ed., Landolt-B¨rnstein - Group I Elementary Particles, Nuclei and Atoms 24 : Elementary Particles, Nuclei and Atoms (Springer-Verlag Berlin Heidelberg, 2009).
Horvath [2020] A. G. Horvath, Arnold Mathematical Journal 6, 1 (2020).
Mazzolo et al. [2003] A. Mazzolo, B. Roesslinger, and C. M. Diop, Ann. Nucl. Energy 30, 1391 (2003).
Roberts and Torquato [1999] A. P. Roberts and S. Torquato, Phys. Rev. E 59, 4953 (1999).
Zoia et al. [2019] A. Zoia, C. Larmier, and D. Mancusi, EPL (Europhysics Letters) 127, 20006 (2019).
El Khaldi and Saleeby [2017] K. El Khaldi and E. G. Saleeby, Monte Carlo Methods and Applications 23, 13 (2017).
Zhang [1999] G. Zhang, Trans. Amer. Math. Soc. 351, 985 (1999).
Almenas and Lee [1992] K. Almenas and R. Lee, Nuclear Engineering: An Introduction (Springer Berlin Heidelberg, Berlin, Heidelberg, 1992).
Sjostrand [2002] N. Sjostrand, Ann. Nucl. Energy 29, 1607 (2002).
Beckurts and Wirtz [1964] K. H. Beckurts and K. Wirtz, Neutron Physics Translated by L. Dresner (Springer Verlag OHG Berlin, 1964).
Liverhant [1960] S. Liverhant, Elementary Introduction to Nuclear Reactor Physics (Wiley, 1960).
Rothenstein [1960] W. Rothenstein, Nucl. Sci. Eng. 7, 162 (1960).
Leoncini et al. [2018] X. Leoncini, A. Vasiliev, and A. Artemyev, Phys. D 364, 22 (2018).
Leslie et al. [1965] D. C. Leslie, J. G. Hill, and A. Jonsson, Nucl. Sci. Eng. 22, 78 (1965).
Christy et al. [1992] R. F. Christy, A. M. Weinberg, and E. P. Wigner, Resonance escape probability in lattices, in Nuclear Energy, edited by A. M. Weinberg (Springer Berlin Heidelberg, Berlin, Heidelberg, 1992) pp. 475–486.
Laramore et al. [2018] D. Laramore, M. P. Pfeifer, J. Lindstrom, and H. Bindra, Ann. Nucl. Energy 111, 255 (2018).
Todorov and Bloch [2017] P. Todorov and D. Bloch, The Journal of Chemical Physics 147, 194202 (2017).
Elmaghraby et al. [2020] E. K. Elmaghraby, S. Abdelaal, A. Abdelhady, S. Fares, S. Salama, and N. Mansour, Nucl. Instrum. Meth. Phys. Res. A 949, 162889 (2020).
Espinosa-Paredes et al. [2008] G. Espinosa-Paredes, J. B. Morales-Sandoval, R. Vázquez-Rodríguez, and E.-G. Espinosa-Martínez, Ann. Nucl. Energy 35, 1963 (2008).
Tzika and Stamatelatos [2004] F. Tzika and I. Stamatelatos, Nucl. Instrum. Meth. Phys. Res. B 213, 177 (2004).

Graphical Abstract