Neutron total cross section calculation within the framework of quasi-harmonic approximation

In this section, we briefly introduce the scattering theory. The assumptions made in the process of formula development for cross section calculation are emphasised. More detailed introduction to the theory can be found in [3, 23].

Given by Van Hove [1], the double differential scattering cross section can be expressed as

$$\begin{eqnarray}\displaystyle \frac{{\partial }^{2}\sigma }{\partial {E}^{{\prime} }\partial {\rm{\Omega }}} & = & {b}^{2}\sqrt{\displaystyle \frac{{E}^{{\prime} }}{E}}\displaystyle \int \displaystyle \int G(r,t)\exp ({\rm{i}}\kappa r)\exp ({\rm{i}}\omega t){\rm{d}}r{\rm{d}}t\\ & = & {b}^{2}\sqrt{\displaystyle \frac{{E}^{{\prime} }}{E}}\displaystyle \int I({\boldsymbol{\kappa }},t)\exp ({\rm{i}}\omega t){\rm{d}}t\\ & = & {b}^{2}\sqrt{\displaystyle \frac{{E}^{{\prime} }}{E}}S({\boldsymbol{\kappa }},\omega ),\end{eqnarray} \tag{ 1 }$$

where b is the scattering length, ω is the phonon frequency, ${\boldsymbol{\kappa }}$ is the momentum transfer, G is the so-called time-dependent pair-correlation function, which is the probability of finding an atom at position r at time t. I is the intermediate scattering function, which is the Fourier transform of $G(r,t)$ in space. $S({\boldsymbol{\kappa }},\omega )$ is the Fourier transform of $I({\boldsymbol{\kappa }},t)$ in time. I can be expressed as

$$\begin{eqnarray}&&I({\boldsymbol{\kappa }},t)=\langle \displaystyle \sum _{d=0}^{N}\exp [-{\rm{i}}{\boldsymbol{\kappa }}\cdot {u}_{0}(r,0)]\exp [{\rm{i}}{\boldsymbol{\kappa }}\cdot {u}_{r}({r}_{d},t)]{\rangle }_{T},\end{eqnarray} \tag{ 2 }$$

where d is the atom index, u is the instantaneous displacement.

Strictly speaking, the exponential factors in equation (2) do not commute, and the evaluation of equation (2) in anharmonic crystals was given by Maradudin and Fein [24]. In the limiting case where the crystal is harmonic, only the leading terms of the anharmonic expression need to be evaluated. In that case, equation (2) can be simplified as [24, 36],

$$\begin{eqnarray}&&I({\boldsymbol{\kappa }},t)=\displaystyle \sum _{r}\exp (-{\rm{i}}{\boldsymbol{\kappa }}\cdot r)\exp \langle [{\boldsymbol{\kappa }}\cdot u(r,0)][{\boldsymbol{\kappa }}\cdot u(r,t)]{\rangle }_{T}.\end{eqnarray} \tag{ 3 }$$

In the special case that r = 0, equation (3) becomes the intermediate function for the incoherent scattering, I_s, which describes the correlation of an atom with itself at different time. The Fourier transform of I_s in time is strictly valid to represent the incoherent scattering of the atom in a general harmonic lattice, implying that the total incoherent scattering of a lattice can be decomposed into the independent contributions from individual atoms.

The Debye–Waller factor, $\exp (-2W)$, is a measure of the vibrational effects to the elastic scattering cross sections. In harmonic crystals, where the displacements are small compared to atomic distances, it can be expressed as

$$\begin{eqnarray}&&{W}_{d}=\displaystyle \frac{{\rm{\hslash }}}{4{M}_{d}N}\displaystyle \sum _{j}\displaystyle \frac{| {\boldsymbol{\kappa }}\cdot {\hat{{\boldsymbol{e}}}}_{{dj}}{| }^{2}}{{\omega }_{j}}\langle 2{n}_{s}+1\rangle ,\end{eqnarray} \tag{ 4 }$$

where $\hat{{\boldsymbol{e}}}$ is the phonon eigenvector, n_s is the Bose–Einstein function, M_d is the mass of the d_th nucleus, j is the normal mode index, N is the number of unit cells.

In a cubic Bravais crystal (i.e. elemental cubic crystal), the square sum of the eigenvectors at any k-point in the Brillouin zone is unity. In that case, atomic displacement is isotropic, so equation (4) can be simplified as

$$\begin{eqnarray}&&W=\displaystyle \frac{{{\rm{\hslash }}}^{2}{\kappa }^{2}}{4M}{\int }_{0}^{{\omega }_{m}}\displaystyle \frac{1}{\omega }\coth ({\rm{\hslash }}\omega /2{k}_{b}T)\rho ({\rm{\hslash }}\omega ){\rm{d}}{\rm{\hslash }}\omega ,\end{eqnarray} \tag{ 5 }$$

where ρ is the VDOS, which is a probability density function, and can be expressed as

$$\begin{eqnarray}&&\rho ({\rm{\hslash }}\omega )=\displaystyle \frac{1}{N}\displaystyle \sum _{{\boldsymbol{q}}j}\delta [{\rm{\hslash }}{\omega }_{{\boldsymbol{q}}j}-{\rm{\hslash }}\omega ],\end{eqnarray} \tag{ 6 }$$

where ${\rm{\hslash }}\omega$ is the vibrational energy. For crystals in the ideal powder form, equation (5) is expected to be close to exact in a general lattice, because the isotropically randomised crystallite orientations largely remove the angular dependence effects.

Sjölander [3] showed that the self-intermediate scattering function in a harmonic cubic Bravais lattice can be expressed in a simple form as

$$\begin{eqnarray}&&{I}_{s}({\boldsymbol{\kappa }},t)=\exp \left[-\displaystyle \frac{{{\rm{\hslash }}}^{2}{\kappa }^{2}}{2M}\gamma (0)\right]\exp \left[\displaystyle \frac{{{\rm{\hslash }}}^{2}{\kappa }^{2}}{2M}\gamma (t)\right],\end{eqnarray} \tag{ 7 }$$

where

$$\begin{eqnarray}&&\gamma (t)={\int }_{0}^{\infty }\displaystyle \frac{\rho ({\rm{\hslash }}\omega )}{\omega }\left[\coth \left(\displaystyle \frac{{\rm{\hslash }}\omega }{2{k}_{b}T}\right)\cos (\omega t)+{\rm{i}}\sin (\omega t)\right]{\rm{d}}\omega .\end{eqnarray} \tag{ 8 }$$

Comparing equation (5) and equation (8), it is straight forward to obtain the Debye–Waller factor in an isotropic material as $2W={{\rm{\hslash }}}^{2}{\kappa }^{2}\gamma (0)/2M$, where $\gamma (0)$ is linear with $2W$, so with the mean-square displacement. It is possible to integrate the Fourier transform of equation (7) directly to obtain the total inelastic cross section [37]. The challenge is that the integrand of the equation is highly oscillating and the error is less straight-forward to estimate. Moreover, the information of the contribution from a specific number of phonon scattering is missed in this way, it is then disabled the possibility to study the contribution from multiple phonons scattering. We therefore are more in favour of evaluating this equation in the series of Taylor expansion [3, 38].In that method, the second factor in equation (7) is broken down into Taylor series as,

$$\begin{eqnarray}&&{I}_{s}(\kappa ,t)={{\rm{e}}}^{-2W}\displaystyle \sum _{n=0}^{n}\displaystyle \frac{1}{n!}{\left[\displaystyle \frac{{{\rm{\hslash }}}^{2}{\kappa }^{2}\gamma (t)}{2M}\right]}^{n}.\end{eqnarray} \tag{ 9 }$$

This technique is known as phonon expansion, where each term of n represents the contribution from the scattering process involving n phonons. It is also often labelled as the multiple phonon correction in the case of removing the multiple phonon scattering background from the raw data of neutron scattering measurements (e.g. in [39]).

The Fourier transform in time for each power of $\gamma (t)/\gamma (0)$ is [3]

$$\begin{eqnarray}{H}_{n}(\epsilon )=\left\{\begin{array}{ll}\delta (\epsilon ) & ,\quad n=0\\ {[2\epsilon \gamma (0)\sinh (\epsilon /2{k}_{b}T)\exp (\epsilon /2{k}_{b}T)]}^{-1}\rho (\epsilon ) & ,\quad n=1\\ \displaystyle \int {H}_{1}({\epsilon }^{{\prime} }){H}_{N-1}(\epsilon -{\epsilon }^{{\prime} }){\rm{d}}{\epsilon }^{{\prime} } & ,\quad n\gt 1\end{array}\right..\end{eqnarray} \tag{ 10 }$$

The partial incoherent scattering kernel ${S}_{n}({\boldsymbol{\kappa }},\epsilon )$ for each phonon term is then

$$\begin{eqnarray}&&{S}_{n}({\boldsymbol{\kappa }},\epsilon )={{\rm{e}}}^{-2W}\displaystyle \frac{{(2W)}^{n}}{n!}{H}_{n}(\epsilon )\end{eqnarray} \tag{ 11 }$$

The intermediate function of the coherent inelastic scattering is the case where both r and t are finite in equation (3). To ease the calculation, the incoherent approximation is employed in this work. In this approximation, the coherent inelastic scattering in a general lattice is treated as if it was the incoherent inelastic scatterings in a cubic Bravais lattice. In other words, the self-intermediate scattering function, equation (7), is chosen to represent the coherent scattering as well. The validity of this approximation is discussed separately in section 4.1. To apply this approximation in the total cross section calculation, b in equation (1) is substituted by ${b}_{i}^{2}+{b}_{c}^{2}$, where b_i and b_c are the bound incoherent and coherent scattering length, respectively.

The incoherent approximation is known to be unsuitable for the elastic coherent scattering (i.e. the Bragg diffraction), because this type of scattering depends highly on the crystal structure. In high quality single crystals, the Bragg diffraction is neglectable in the total cross section, since the diffraction condition is approaching a Dirac delta function, to which very few neutrons can fulfil. On the other hand, the materials used in the nuclear industry and many other places are often in their powder or polycrystal forms, where, ideally, the crystallites are isotropically distributed, so the Dirac delta function is significantly broadened. As a result, each symmetry-equivalent reflection families introduces a stepping rise, which is known as the Bragg edge, at the neutron wavelength that is twice of the corresponding atomic plane spacings. Unfortunately, because of the lack of handling crystal symmetries in our numerical implementation, the cross section of elastic coherent scattering is not included in our calculated total cross sections. However, this component is very distinctive and can be easily identified when comparing our calculated cross sections with those from transmission measurements.

To simplify the evaluation of total cross section, we define the momentum and energy transfers in the dimensionless form, as shown in equation (12).

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{E+{E}^{{\prime} }-2\mu \sqrt{{E}^{{\prime} }E}}{{k}_{b}T}\qquad \quad \beta =\displaystyle \frac{{E}^{{\prime} }-E}{{k}_{b}T}.\end{eqnarray} \tag{ 12 }$$

The total inelastic scattering cross section is the integral of the double differential cross section (i.e. equation (1)) over the boundary that defined by equation (13). It can be seen that the integration boundary is only defined by $E/{k}_{b}T$, where the factor of temperature is introduced by the variable change from $(Q,\omega )$ to $(\alpha ,\beta )$ for which the Jacobian is k_bT. Within this boundary, the scattered angle cosine μ is in the range $[-1,1]$, and the scattered energy ${E}^{{\prime} }$ is non-negative.

$$\begin{eqnarray}&&{(\alpha -\beta )}^{2}=\displaystyle \frac{4E\alpha }{{k}_{b}T}.\end{eqnarray} \tag{ 13 }$$

The energy transfer distribution function P is the integral over α

$$\begin{eqnarray}&&P(\beta )={k}_{b}T\exp (-\beta /2){\int }_{\alpha -}^{\alpha +}\displaystyle \sum _{n=1}^{\infty }{S}_{n}(\alpha ,\beta ){\rm{d}}\alpha ,\end{eqnarray} \tag{ 14 }$$

where S_n is non-symmetric and obeys the detailed balance, and has the analytic property of $S(Q,\omega )=\exp (-{\rm{\hslash }}\omega /{k}_{b}T)S(Q,-\omega )$. In this work, this integral is considered to be converged when the nth term contributes less than 0.01% to the sum of the lower terms For the calculations performed on the materials that studied in this work, this integral is well converged below 50 orders, because the neutron energies of interest are relatively small and the atoms are heavy. Equation (14) also indicates that the temperature dependency of the cross section is linear in first approximation. The total inelastic scattering cross section ${\sigma }_{i}$ at incident energy E is

$$\begin{eqnarray}&&{\sigma }_{i}(E)=\displaystyle \frac{\pi ({b}_{i}^{2}+{b}_{c}^{2})}{E}{\int }_{-E/{k}_{b}T}^{\infty }P(\beta ){\rm{d}}\beta .\end{eqnarray} \tag{ 15 }$$

The scattering lengths, b_i and b_c, used in this work are taken from the table compiled in [40].

In the epithermal energy region, ${\sigma }_{i}(E)$ approaches the limiting free incoherent scattering cross section $4\pi {a}_{i}^{2}$, where a_i is the free incoherent scattering length [41].

For the incoherent elastic scattering, indicated by equation (10) and (11), the scattering function has a simple form as

$$\begin{eqnarray}&&{S}_{0}({\boldsymbol{\kappa }},\epsilon )={{\rm{e}}}^{-2W}\delta (\epsilon ).\end{eqnarray} \tag{ 16 }$$

The Dirac delta function here emphasises zero energy transfer in this process. Integrating equation (16) over 4π and scaling by the incoherent cross section, we obtain the total cross section.

$$\begin{eqnarray}&&{\sigma }_{i,e}=4\pi {b}_{i}^{2}{\left(\displaystyle \frac{2B}{{\lambda }^{2}}\right)}^{-1}\left[1-\exp \left(-\displaystyle \frac{2B}{{\lambda }^{2}}\right)\right],\end{eqnarray} \tag{ 17 }$$

where $B=8{\pi }^{2}\langle {u}^{2}\rangle$, it is the so-called B factor in crystallography.

The dependence of phonon energy on the unit cell volume is traditionally described by the Grüneisen parameter defined in equation (18).

$$\begin{eqnarray}&&{\gamma }_{{\boldsymbol{\kappa }},j}=-\displaystyle \frac{V}{\omega ({\boldsymbol{\kappa }},j)}\displaystyle \frac{\partial \omega ({\boldsymbol{\kappa }},j)}{\partial V}.\end{eqnarray} \tag{ 18 }$$

At a fixed pressure, the Grüneisen parameter of a phonon is often regarded as a constant, so that the deviation of phonon energy is linear with the deviation of the volume. Therefore, we can obtain the phonon energy as

$$\begin{eqnarray}&&{\omega }_{{\boldsymbol{q}}j}({V}_{T})={\omega }_{{\boldsymbol{q}}j}({V}_{o})\left(1-\displaystyle \frac{{V}_{T}-{V}_{o}}{{V}_{o}}{\gamma }_{{\boldsymbol{q}}j}\right),\end{eqnarray} \tag{ 19 }$$

where V_o and V_T are the reference volume and the volume at a finite temperature T, respectively.

Generally, two methods can be used to obtain the VDOS defined in equation (6). For applications that are interested in high energy resolution, the tetrahedron method [42] is often employed. For a crystal at high temperatures (or even at room temperature), the sub-meV structures in the VDOS may not be observed, because as a type of quasi-particle, the phonon energy is broadened by the presence of different isotopes, phonon–phonon and electron–phonon scatterings. As a practical and simpler solution, the delta function in equation (6) is often replaced by a Gaussian. However, in this way, artificial contribution can be introduced to the part of VDOS where energy is approaching zero. It can be observed in equation (10) that the H₁ function is extremely divergent in that region. Therefore, the width of the Gaussian kernel need to be tuned carefully to be small enough to avoid large artifact at low energies; while large enough to ensure no artificial structures can be produced at higher energies. In this work, the delta function is implemented in the same form as a classical damped oscillator [43]. The analytic property of such broadening function is naturally in consistent with the requirement in this work: the VDOS is vanishing when energy is approaching zero. Therefore, we can write the temperature dependent VDOS as

$$\begin{eqnarray}&&\rho ({\rm{\hslash }}\omega )=\displaystyle \frac{1}{{\rm{\hslash }}N}\displaystyle \sum _{{\boldsymbol{q}}j}\displaystyle \frac{{\omega }^{2}{{\rm{\Gamma }}}_{{\boldsymbol{q}}j}({\rm{\hslash }}\omega ,T)}{{[{\omega }^{2}-{\omega }_{{\boldsymbol{q}}j}^{2}({V}_{T})]}^{2}+{\omega }^{2}{{\rm{\Gamma }}}_{{\boldsymbol{q}}j}^{2}({\rm{\hslash }}\omega ,T)},\end{eqnarray} \tag{ 20 }$$

where ${{\rm{\Gamma }}}_{{\boldsymbol{q}}j}$ is the full-width-half-maximum of the broadening, and it is set to 0.5 meV for all the phonons. We justify this simple treatment of the broadening in section 3.4.

Detailed DFT calculation procedures of phonon mode Grüneisen parameters and crystal thermal expansion in the framework of QHA are presented in the appendix. Using the calculated results of Al, we compare the softening of four phonons with the experimental data from [44] in figure 1. The calculated phonon energies differs from measured energies by less than 5% in all cases, except when temperature is very close to melting point (i.e. 12% for the k-point $q/{q}_{\max }=0.3$ at 933.5 K).

Zoom In Zoom Out Reset image size

Figure 2 shows the calculated dispersion curves of Al at 80 and 300 K using equation (19). Our QHA calculation predicts that the volume of Al unit cell expanded about 0.0769% at 80 K and 1.35% at 300 K. As the mode Grüneisen parameters of Al are roughly 2 (see figure A2(a) in the appendix), the solid curve of 80 K is expected to be hardly distinguishable from the dispersion curve predicted by the harmonic approximation; the phonons at 300 K are expected to be softened by roughly 3%, which is observable in the regions near the peaks of the curves. The calculated curves closely follow the structure of the data, e.g. the relatively small slopes of the blue curves for one of the transverse phonon branches in the regions that between K and X and that close to L. The absolute differences between the measured and calculated phonon energies are probably close to the error introduced in the process of extracting the data points. As the dispersion curve at 300 K is clearly distinguishable from that at low temperature, we expect QHA to be effective already at this temperature for the inelastic scattering in Al.

Zoom In Zoom Out Reset image size

In this section, we benchmark our calculated total cross sections against the data compiled in EXFOR [35].

The predicted and measured cross sections of Al are compared in figure 3. The sudden rise of the measured cross sections at about 5 meV is due to the coherent elastic scattering in the powder sample used in the measurement. Such stepping behaviour can not be observed in the 298 K data from Grabcev, which is measured on a piece of high quality single crystal sample. It is notable that the coherent elastic and inelastic scatterings are two competing processes (see e.g. [16]). When the coherent elastic scattering becomes valid above its energy threshold, a moderate deviation on the VDOS may change the contribution ratio of the processes, but have a minor impact to the total cross sections. In Al, the capture cross section is the dominant process at low temperatures and low incident energies. It contributes over 97% to the total cross sections when the incident energy is 1 meV at 100 K. The predicted cross sections at 100 K systematically underestimate the measurement by about 0.06 b. Such discrepancies may be due to large capture cross section impurities in the sample. It can be observed that the measured points are all situated slightly above the curve predicted by the QHA model at higher temperatures.

Zoom In Zoom Out Reset image size

The comparisons on Mg are shown in figure 4. All the energies in the measurement are below the coherent elastic scattering threshold at about 2.66 meV, which is defined by the $\{1,0,0\}$ atomic planes, so the Bragg edge can not be observed in the measured data. Similar to the case of Al, good agreement is observed at all temperatures, most of the measured points intersect with the predicted curves of the QHA model. The static model curves are systematically lower than the QHA curves in all cases, and the underestimation increases with temperature, indicating the phonon softening contributes positively to the cross sections. It can be understood as that with the softening of the VDOS, $\gamma (t)$ expressed in equation (8) is increased, and subsequently leads to the increase of S. As the integration boundary on S for the total cross section is unchanged, the final cross section is increased. At 781 K, the total cross sections are underestimated by about 7% around 1 meV.

Zoom In Zoom Out Reset image size

As beryllium lattice is less sensitive to temperature, we observe that the differences between the static and QHA models are only distinguishable above 700 K. Figure 5 shows the comparisons between the predictions of the QHA model and the measured data at 100 and 300 K along with the weak capture cross section. At 300 K, the predicted curve is compared with data taken from two independent measurements. About 10% overestimation in the prediction can be observed below 1 meV. The cross sections at high temperatures are compared in figure 6. Similar overestimations are observed in Palevsky's data at 440 K. Predictions from both of the QHA and static models slightly underestimated the measured data from Kanda at higher temperatures, but most of the predictions are within the experimental uncertainties.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We emphasise that the small discrepancies between the calculated and measured cross sections observed in this section may be already close to, if not small than, the uncertainty introduced by the incoherent approximation. Nevertheless, it is clear that the QHA model can effectively capture the additional growth of the cross section introduced by phonon softening, at least in the crystals studied in this work. We discuss the validity of this approximation in section 4.1.

The increase of calculated cross sections due to phonon softening can readily be seen in the results presented in section 3.2. Because of the insignificant incoherent elastic contribution, those increases are almost entirely contributed by the inelastic scattering process, which is the principle process of neutron slowing down. In this section, we quantitatively study the contribution originated from phonon softening.

Figure 7 shows the quasi-harmonic contribution to the integral inelastic scattering cross section as a function of temperature and incident energy. Such contribution is defined as the ratio between the increased cross section, when applying the QHA in the place of the static model, and the cross section from QHA. The upper limit of the temperature in the figures is about 90%–95% of the temperature of the corresponding material melting point. Overall, the contribution to the cross section of Be is relatively weak. It can be observed from the colour bar, the maximum value in the figure is less than 7.5%. On the other hand, Al is very sensitive to temperature, the maximum discrepancy is about 19% at 900 K, and 7% at room temperature.

Zoom In Zoom Out Reset image size

As expected, the global maximum for each material is observed at the highest temperature calculated. However, the energy location of the global maximum is material dependent, which indicates the dependency to the energy range of the phonon distribution. Observed from the figures, the locations are at about 0.06 eV, 0.016 eV and 0.02 eV in Be, Mg, and Al, respectively. At low neutron energies, in the $1/v$ region, the quasi-harmonic contribution grows with temperature, but is insensitive to the incident energy. Above a few meV, the contribution also shows a strong dependency on the incident energy, likely due to the shift of the global minimum of the total cross section (see e.g. cross section curves at 800 K in figure 3). In the region close to eV, the predictions from both of the harmonic and quasi-harmonic models should approache the limiting free scattering cross section, regardless of temperature. Consistently, it can be observed in figure 7 that the quasi-harmonic contribution is diminishing at incident energies close to 1 eV.

In this section, we assess the simple treatment of using a fixed global phonon linewidth in equation (20). Showed by Tang et al [46] in Al, accounting for the finite lifetime due to the three-phonon scattering process, the broadened quasi-harmonic VDOS agrees well with the VDOS that is derived from inelastic neutron scattering measurements at 750 K. For this particular metal, this finding can be interpreted as, firstly, the difference between the bare-phonon energy and the renormalized phonon energy is small in comparison with the energy shift due to volume expansion; secondly, the dominant contributor for phonon lifetime is the process of three-phonon scattering. Therefore, for the purpose of assessing our simple treatment for the linewidth, we study the effects of including the dominant contributor of the phonon linewidth in equation (20) in Al.

The calculated phonon linewidths at 300 and 800 K are shown in figure 8. At 300 K, phonons are highly concentrated around 35 meV, where the mean line width is about 1.45 meV. Such observation is in excellent agreement with the calculated results given in [46]. At both temperatures, the energetic longitudinal phonons have much broader linewidth than in the other two phonon branches. Phonon linewidths at 800 K are about twice of the those at 300 K, and the overall structures are broadly similar.

Zoom In Zoom Out Reset image size

We found that when including the three-phonon scattering process in equation (20), the atomic mean-square displacement is reduced by a few percent at high temperatures. The calculated B factors from different models are compared with experimental data [47] in figure 9. As expected, the results from the harmonic model show linear relationship with temperature. The other two models show additional higher order T behaviour. However, none of them increase as fast as the experimental results with temperature. The results based on the broadened VDOS are systematically lower than that from the QHA model. As for the total cross sections, the results based on the broaden VDOS systematically underestimate the results from the QHA model. However, the maximum discrepancy we observed at 800 K is about 1% in the energy range between 10⁻⁵ and 1 eV.

Zoom In Zoom Out Reset image size

Based on the findings in this section, we conclude that our total cross section model is insensitive to phonon linewidth.

One of the major sources of the uncertainty in our model is the incoherent approximation, in which the inelastic scatterings are treated as purely incoherent. So the angular distribution does not contain any structure peaks. For neutron transport problems in an infinite scattering media, it is arguable if the interference effects are important at all. After a number of scatterings, the sharp peaks in the angular distributions arising from the single phonon scattering smear out. Therefore, we expect that the angular distribution modelled in the incoherent approximation is realistic when the material size is at least a few times of the inelastic mean free path.

However, the absolute error of the incoherent approximation in total cross section calculation is lattice dependent and not well understood. The main uncertainty is most likely introduced by the treatment to the single phonon scattering, which is the dominant inelastic process at low energies. For example, figure 10 shows the contribution of the single phonon scattering to the total inelastic cross section in Al. The single phonon scattering contribution is significant in the low energy region, where the flat plateaus are formed due to the $1/v$ behaviour of the inelastic processes. Its contribution gradually decreases with the increased temperature, from about 98% at 100 K to about 73% at 800 K. Also the decay of the plateaus starts from a lower energy at increased temperature, from at a few tens of meV at 100 K to at a few meV at 800 K. This figure clearly shows the promotion of multiple phonon scatterings as a result of increased phonon population at higher temperatures. More importantly, it can be observed that the contribution from the single phonon scattering is significant in the energy region of interest in this work.

Zoom In Zoom Out Reset image size

Suggested in [48, p 314], when applying the incoherent approximation to the energy differential single phonon scattering cross sections, the error rarely exceeds 20% if the integral domain of the solid angle is sufficiently large. We believe the net effect to the total cross section is effectively reduced in the process of integrating the differential cross section. Indicated by the good performance of the incoherent approximation shown in section 3.2, we expect the error introduced by this approximation is only a few percent for the materials studied in this work.

However, significant errors introduced by incoherent approximations have been reported previously in Be [49] and graphite [50] near room temperature, where the anharmonic effects are insignificant. Both [49, 50] calculate the coherent single phonon scattering explicitly, and employ incoherent approximation for the higher order scatterings. We shall refer their model as the quasi-coherent approximation.

For the graphite case, it has been reported that the calculated total cross sections using the incoherent approximation are only about 50% of those from measurement below the threshold energy of coherent elastic scattering [50]. Indicated by Bowman et al in a recent study [51], the measured data in [50] may be based on reactor grade graphite samples, of which the density could be lower than 80% of the density of ideal pyrolytic graphite due to the high porosity concentration. It has been shown in [51] that the measured total cross sections from high quality graphite closely agree with those predicted by incoherent approximation. Very recently, Hawari and Gillete [52] showed that using the VDOS calculated by molecular dynamics techniques, even the cross section of reactor grade graphite can be accurately calculated in the incoherent approximation.

Borgonovi and Sprevak [49] reported that the incoherent approximation can overestimate the cross section of Be by up to 20%. We compared the cross sections predicted in this work with the experimental data and the predictions by Borgonovi and Sprevak in figure 11. The quasi-coherent approximation curve from [49] overlay with experimental points. However, the incoherent approximation cross section curve (in red) from [49] shows less good agreement with the measured cross section than the curve (in green) calculated in this work. Indicated by figure 7 in [53], the phonon dispersion curve predicted by the lattice dynamics model used in [49] underestimates the measured data typically by 10%–20%, which can explain the overestimation observed in the red curve. So the good performance given by Borgonovi and Sprevak's quasi-coherent approximation is likely the result of the error in their numerical implementation cancelles out the error in the lattice dynamics model.

Zoom In Zoom Out Reset image size

Both of the works [49, 50] likely suffer from the error introduced by the numerical implementations. Nevertheless, the model they proposed is certainly more physically sound than the incoherent approximation, as the single phonon scattering can be treated much more realistically. Facilitated by modern ab initio methods and improved numerical techniques (e.g. the tetrahedron method), overcoming the technical hurdles for such advanced model and beyond should be feasible.

From the results presented in section 3.2 we can see that the increase of total scattering cross sections in the materials studied in this work are closely correlated to temperature as well as to the softening of phonon energy. Under these circumstances, studying neutron transmission cross sections of a crystal at different sample environments may provide a useful insight to the deviation of VDOSs of the material. As an example, we recover the VDOS of Al at 800 K from experimental transmission cross sections of Al powder at 100 and 800 K shown in figure 3.

We formulate the recovering model into a simple optimisation problem, assuming the cross sections at 100 K can be correctly calculated by using the harmonic VDOS. The VDOS deviation at 800 K can be described by the modification of the harmonic VDOS, of which the energy axis is scaled by a scalar value X₁ , which is a parameter to characterise phonon softening in the place of mode Grüneisen parameters and volume expansion (see equation (19)), as a measure of the thermal dependency of the cross section. Finally, the discrepancies between calculated and measured cross sections are insensitive to neutron energy and temperature, so they can be described by a second scalar value X₂. The quantity that we are trying to minimise is the sum of the mean squared errors between the measured cross sections and the corresponding calculated results at 100 and 800 K. Experimental cross section above $3\times {10}^{-3}\,\mathrm{eV}$ are removed in the calculation to eliminate coherent elastic scattering contributions. The higher order anharmonicity and experimental background may be oversimplified in our parameterisation, but this optimisation problem is primarily formulated to capture the softening of VDOS. We believe it can serve our purpose of demonstration.

We use the particle swarm optimisation [55] to find the best solution. The error is minimised when the softening parameter X₁ equals 0.862 and the constant cross section error X₂ equals 0.0508 b. The recovered VDOS is compared with the VDOS from our direct DFT calculations and that from a measurement at 775 K [54] in figure 12. The phonon softening effects can be observed in all the measured and calculated VDOS. The recovered VDOS shows good agreements with the VDOS obtained from the direct DFT QHA calculation, because our recovering model is quite quasi-harmonic. For instance, the finite lifetime of the phonon is not considered, and the anharmonic contribution to the Debye–Waller factor is also disregarded. Our recovered and the calculated QHA VDOS both overestimate the softening of the peak of the longitudinal phonons, of which the centre is at about 36meV in the harmonic VDOS. Such overestimation is expected. An excellent discussion regarding the QHA behaviour at high temperature can be found in [56].

Zoom In Zoom Out Reset image size

VDOS is very important for modelling thermal neutron transport and material thermodynamics. The direct measurement of VDOS by neutron scattering techniques can provide an initial VDOS before the VDOS deviation; while the VDOS deviation obtained from analysing transmission cross sections can be used to correct the background in the incoherent approximation. So these two methods are complementary. In principle, it is possible to measure the double differential cross sections and the total cross sections can be measured in one measurement. After all, enabling transmission measurement in an existing TOF neutron instrument amounts to introducing some relatively trivial add-on detectors in front of the instrument beamstop.