Resonance Scattering Treatment with the Windowed Multipole Formalism

Abstract

A new method for directly sampling the resonance upscattering effect is presented. Alternatives have relied on inefficient rejection sampling techniques or large tabular storage of relative velocities. None of these approaches, which require pointwise energy data, are particularly well suited to the windowed multipole cross-section representation. The new method, called multipole analytic resonance scattering, overcomes these limitations by inverse transform sampling from the target relative velocity distribution where the cross section is expressed in the multipole formalism. The closed-form relative speed distribution contains a novel special function we deem the incomplete Faddeeva function, and we present the first results on its efficient numerical evaluation.

Keywords:

• Resonance upscatter
• windowed multipole

I. INTRODUCTION

Early continuous-energy Monte Carlo neutron transport programs sampled scattering from nuclei in thermal motion assuming that the scattering cross section is effectively constant within the scattering kernel.^[1] However, as Refs. [2,3] detail, the resulting scattering kernel implied by the constant cross-section approximation may be far from the actual double-differential cross section near a scattering resonance. As shown in Refs. [4,5], the resulting error tends to cause a worst-case 11% underestimation of the Doppler feedback coefficient in a pressurized water reactor (PWR), with even larger discrepancies in high-temperature gas-cooled reactor problems.

As exhibited by the Power Reactor Analysis using GPU-based Monte Carlo Algorithm (PRAGMA) project,^[6,7] the conventional methods for treating this effect leave something to be desired on graphics processing unit (GPU) architectures, which constitute the majority of computational power on most leading supercomputers. Due to the unique architecture of the GPU, algorithmic modifications to standard Monte Carlo algorithms for neutron tracking can tangibly accelerate computation.^[7] In the same direction, we herein present a GPU-friendly method for handling resonance upscattering when the windowed multipole (WMP)^[8] formalism is employed to represent cross sections. In particular, the heuristic for fast GPU code is to avoid rejection sampling and accesses to distantly spaced places in memory, which the new method achieves. To give context to the new method, we first recall some of the conventional methods for modeling resonance upscattering.

The Doppler Broadening Rejection Correction (DBRC)^[9] was one of the early proposed techniques to treat the effect of strong variations of the interaction cross section within the energetic vicinity of a scattering neutron, whereas $\mathit{S}(\mathit{\alpha },\mathit{\beta })$ tables had been used prior.^[10] The method has since been implemented in numerous continuous-energy Monte Carlo neutron transport programs^[9,11–13] and has been shown to successfully model the resonance upscattering effect. However, DBRC suffers from rejection probabilities as high as 99.995%^[14] for neutron energies in the vicinity of a resonance.

The weight correction method (WCM)^[4] also successfully models the effect of resonances on the double-differential cross section of nuclei in thermal motion. This method adjusts the weight of particles to provide numerically correct results even when the constant cross-section double-differential free gas distribution is employed. WCM carries the same benefit as our newly proposed method of not requiring an additional rejection loop or tables; however, adjustments to the particle weights introduce substantial variance to the overall Monte Carlo simulation, thus degrading estimates on quantities of interest.^[6]

Another technique known as target motion sampling^[15] can be used to model the resonance upscattering effect. However, this method is only applicable to Monte Carlo neutron transport programs employing the delta tracking technique. Unfortunately, the performance of delta tracking appears to be lackluster on GPUs.^[16]

The relative velocity sampling (RVS) method^[14,17] was created to ameliorate the high rejection rates characteristic of the rejection algorithms used to model resonance upscattering. These schemes, in essence, sample probability distributions proportional to $\mathit{f}(\mathit{x})\mathit{g}(\mathit{x})$, where $\mathit{f}(\mathit{x})$ is a distribution and $\mathit{g}(\mathit{x})\in [0,1]$. One then samples from $\mathit{f}(\mathit{x})$ and accepts the sample with probability $\mathit{g}(\mathit{x})$. The RVS method moves the direct sampling from the thermal motion term to the cross-section term, thus worsening the average case rejection rate but massively improving the worst-case rejection rate.

The relative speed tabulation (RST) method^[6] was developed to address the rejection sampling performance impact on GPUs. RST is the first resonance upscattering modeling technique not requiring a rejection loop or bivariate scattering distribution tables. The key observation underlying RST is that the target velocity distribution can be factorized into a marginal relative speed distribution encapsulating the information about the resonances and a simple distribution of the target polar angle conditioned on the target relative speed. Consequently, the univariate relative speed cumulative distribution function (CDF) can be tabulated in select areas of the pointwise cross section, and the conditional polar angle distribution is then directly sampled without requiring any additional data or rejection step.

Despite its simplicity and efficacy on GPUs, the RST method comes with some clear disadvantages. Gigabytes of additional memory are used to store the relative speed CDFs at each energy point and temperature on a pointwise cross-section representation if all nuclides have the resonance upscattering effect treated. The resonance influence on the double-differential cross section is thus restricted for practical reasons to a select few nuclides in the problem to avoid extreme memory usage. Moreover, the method introduces some discretization errors in temperature, although this has been shown to be a reasonable approximation.

If one could avoid pretabulated relative speed distributions, this could substantially reduce the memory needs and avoid costly divergent memory accesses. Reducing serialized memory accesses on GPUs typically leads to significant speedups.

We propose a new method that does just that, and achieves this using the WMP cross-section representation.^[8] Our new method introduces a novel special function we have deemed the incomplete Faddeeva function, which encodes the behavior of the temperature-dependent influence of resonances on the double-differential scattering distribution. It relies on the numerical inversion of the analytic representation of the relative speed distribution under the WMP formalism and polar angle sampling in the same manner as the RST method, but without any precomputed tables.

This contrasts the WMP-based target velocity sampling technique presented in Ref. [18], which is similar in nature to the newly presented method in this work. However, the method in Ref. [18] relies on the separation of the target velocity distribution into a zero-kelvin cross-section component and a Maxwell-Boltzmann component. This method thus requires a numerical inversion step of the integrated scattering cross-section function wrapped in a rejection loop, similar to Ref. [14], but instead using a functional representation of the integrated cross section rather than tabular.

The algorithm presented in Ref. [19] shows how the relative speed can be sampled in the WMP framework similarly to our work. However, Ref. [19] relies on fitting a sum of Gaussian to replace the multipole expansion, thus representing the zero-kelvin scattering cross section in the form

(1) $\mathit{\sigma }(\mathit{E})=\frac{1}{E}\sum _k\sum _j\left[h_{\mathit{s},\mathit{k},\mathit{j}}{e}^{{\left(\sqrt{E}-u_k\right)}^2/w_{\mathit{s},\mathit{k},\mathit{j}}^2}+h_{\mathit{a},\mathit{k},\mathit{j}}{e}^{{\left(\sqrt{E}-u_k\right)}^2/w_{\mathit{a},\mathit{k},\mathit{j}}^2}\right].$(1)

While it remains to be seen that Gaussians can be used to approximate all poles appearing in a WMP library, this proposed approximation introduces a considerable number of degrees of freedom to the problem, with 18 unknowns for each pole. The optimization problem thus encountered is highly nonlinear and nonconvex, leading to fitting difficulties. Additionally, scattering kernels in the formalism in Ref. [19] cannot be straightforwardly differentiated with respect to WMP parameters, which would be needed to perform sensitivity analyses.

Our new method, multipole analytic resonance scattering (MARS), only requires the same WMP data as would be used in a calculation without any treatment of the resonance upscattering effect, thus avoiding the need for additional tables or fitting steps. We demonstrate the new method’s negligible performance overhead compared to other methods for sampling the resonance upscattering effect on CPU architectures, with future work exploring its optimized implementation on GPUs.

II. THEORY

Rothenstein and Dagan^[2] expound on the rigorous Doppler-broadened double-differential cross section. Since the exact expression for the lab frame scattering distribution is quite complicated, authors presenting algorithms to model it typically choose to forgo the lab frame expression and instead reason in terms of joint distributions of the target velocity and direction cosine relative to the direction of projectile motion. After sampling the target speed and direction cosine, standard two-body collision kinematics for elastic scattering are employed, where the center of the mass angular distribution comes from the nuclear data file. The resulting physics matches the complicated expressions of Ref. [2].

Similarly, considering the distribution of collision target velocities, this joint distribution is

(2) $\mathit{f}(\mathit{V},\mathit{\mu })=\mathit{C}{v}_r\mathit{\sigma }(v_r)\mathit{M}(\mathit{T},\mathit{V})f^{\prime }(\mathit{\mu }),$(2)

where $\mathit{M}$ is the Maxwell-Boltzmann distribution of target speeds

(3) $\mathit{M}(\mathit{T},\mathit{V})=\frac{4}{\sqrt{\mathit{\pi }}}{\beta }^3{V}^2{e}^{-{\beta }^2{V}^3},$(3)

the variable with units of inverse velocity $\mathit{\beta }=\sqrt{\frac{A}{2\mathit{kT}}}$ parameterizes the target velocities, and the distribution $f^{\prime }(\mathit{\mu })$ is a uniform distribution between −1 and 1. The other variables are $\mathit{C}$, the distribution’s normalizing constant; $v_r$, the relative speed of the neutron with respect to the target; $\mathit{\sigma }$, the zero-kelvin scattering cross section; $\mathit{A}$, the target mass; $\mathit{k}$, the Boltzmann constant; $\mathit{T}$, the absolute temperature; and $\mathit{V}$, the target speed.

Classically, the approximation that $\mathit{\sigma }(v_r)$ is constant has been employed. However, it was shown in Refs. [3,4] that this approximation is incorrect in the vicinity of resonances, where $\mathit{\sigma }(v_r)$ varies over a few orders of magnitude, preferentially causing scattering with targets of relative velocity more closely matching the scattering resonance peaks. The various aforementioned techniques are all just methods for sampling from the distribution of Eq. (2)(2) $\mathit{f}(\mathit{V},\mathit{\mu })=\mathit{C}{v}_r\mathit{\sigma }(v_r)\mathit{M}(\mathit{T},\mathit{V})f^{\prime }(\mathit{\mu }),$(2) with arbitrary forms of the function $\mathit{\sigma }(v_r)$.

More specific knowledge about the form of $\mathit{\sigma }(v_r)$ can be employed. It has been shown extensively^[20,21] at this point that the cross section is accurately represented as a sum of poles in addition to a low-order Laurent expansion $(\mathit{N}\approx 7)$, namely,

(4) $\mathit{\sigma }(\mathit{E})=\frac{1}{E}\left(\mathrm{\Re }\left[\sum _j\frac{r_j}{p_j-\sqrt{\mathit{E}}}\right]+\sum _{\mathit{n}=0}^N{a}_n{E}^{\mathit{n}/2}\right).$(4)

In fact, for the purposes of sampling the resonance upscattering effect, we claim and later numerically demonstrate that the narrow range of attainable $v_r$ leads the zero-kelvin cross section to be accurately represented as a single pole and a linear term over a sufficiently narrow range of energies:

(5) $\mathit{\sigma }(\mathit{E})=\frac{1}{E}\mathrm{\Re }\left[\frac{r_j}{p_j-\sqrt{\mathit{E}}}\right]+{\sigma }_0+{\sigma }_1\sqrt{\mathit{E}}.$(5)

The ${\sigma }_0$ and ${\sigma }_1$ terms are calculated through a linearization process to avoid some complexities with higher-order polynomial fitting. An algorithm for finding the best values of ${\sigma }_0$ and ${\sigma }_1$ is presented in Sec. II.D.

With this approximation, we use the technique developed in Ref. [6]. Rather than attempting to sample the target speed $\mathit{V}$ and direction cosine $\mathit{\mu }$, one instead samples first the relative velocity $v_r$, and then samples $\mathit{\mu }$ from the distribution of $\mathit{\mu }$ conditioned on $v_r$. The distribution in this form for projectile speed $\mathit{v}$ is thus^[6]

(6) $\mathit{P}(v_r|\mathit{v})=\left(e^{-{\beta }^2{(v-v_r)}^2}-e^{-{\beta }^2{(v+v_r)}^2}\right)v_r^2{\sigma }_0(v_r).$(6)

Next, as previously shown for the Doppler-broadening of the WMP format^[8] and the original full multipole format,^[22] we note the term $e^{-{\beta }^2(\mathit{v}+v_r)}$ to be negligible and use the following:

(7) $\mathit{P}(v_r|\mathit{v})\approx e^{-{\beta }^2{(v-v_r)}^2}{v}_r^2{\sigma }_0(v_r).$(7)

At this point, we write the multipole cross section in terms of the relative velocity rather than in terms of energy. The formula to use is

(8) $\mathit{\sigma }(v_r)=\frac{2}{m_n{v}_r^2}\mathrm{\Re }\left[\frac{r_j}{p_j-\sqrt{2m_n}{v}_r}\right]+{\sigma }_0+{\sigma }_1\sqrt{2m_n}{v}_r.$(8)

If we define the auxiliary variables $\mathit{x}=\mathit{\beta }(v_r-\mathit{v})$ and $\mathit{y}=\mathit{\beta v}$, and insert Eq. (8)(8) $\mathit{\sigma }(v_r)=\frac{2}{m_n{v}_r^2}\mathrm{\Re }\left[\frac{r_j}{p_j-\sqrt{2m_n}{v}_r}\right]+{\sigma }_0+{\sigma }_1\sqrt{2m_n}{v}_r.$(8) into the marginal target collision rate distribution in terms of relative speed [Eq. (7)(7) $\mathit{P}(v_r|\mathit{v})\approx e^{-{\beta }^2{(v-v_r)}^2}{v}_r^2{\sigma }_0(v_r).$(7) ], some algebra reveals that

(9) $\mathit{p}(\mathit{x}|\mathit{y})=e^{-x^2}\left(\mathrm{\Re }\left[\frac{\beta r_j}{\mathit{z}-\mathit{x}}\right]+{\beta }^{-2}{(\mathit{x}+\mathit{y})}^2\left({\sigma }_0+{\sigma }_1\mathit{x}\right)\right),$(9)

where $\mathit{z}=\mathit{\beta }{p}_j-\mathit{y}$, which represents a dimensionless measure of the energy gap between the center of the Maxwell-Boltzmann distribution and the location of the resonance.

At this point, we can consider the CDF for the random variable $\mathit{x}$ (dimensionless target velocity) conditioned on $\mathit{y}$ (dimensionless projectile velocity). Integrating Eq. (9)(9) $\mathit{p}(\mathit{x}|\mathit{y})=e^{-x^2}\left(\mathrm{\Re }\left[\frac{\beta r_j}{\mathit{z}-\mathit{x}}\right]+{\beta }^{-2}{(\mathit{x}+\mathit{y})}^2\left({\sigma }_0+{\sigma }_1\mathit{x}\right)\right),$(9) yields

(10) $\mathit{CP}(\mathit{x}|\mathit{y})={\int }_{-\mathrm{\infty }}^x{e}^{-{x^{\prime }}^2}\left(\mathrm{\Re }\left[\frac{\beta r_j}{\mathit{z}-x^{\prime }}\right]+{\beta }^{-2}{(x^{\prime }+\mathit{y})}^2\left({\sigma }_0+{\sigma }_1{x}^{\prime }\right)\right)\mathrm{d}{x}^{\prime },$(10)

where $\mathit{C}$ is the normalizing constant. After distributing the integral and interchanging the $\mathrm{\Re }$ operator with integration, we obtain

(11) $\begin{array}{rl}& \mathit{CP}(\mathit{x}|\mathit{y})=\mathrm{\Re }\left[\frac{r_j\mathit{\pi }}{\mathit{i}{\beta }^{-1}}\mathit{w}(\mathit{z},\mathit{x})\right]+\frac{{\mathit{\beta }}^{-2}{\sigma }_0}{4}(-2e^{-x^2}(\mathit{x}+2\mathit{y})\\ & +\sqrt{\mathit{\pi }}(1+2y^2)(1+\mathrm{erf}(\mathit{x})))\\ & +\frac{1}{2{\beta }^2}{\sigma }_1{e}^{-x^2}\left(1+{(\mathit{x}+\mathit{y})}^2+\sqrt{\mathit{\pi }}\mathit{y}(1+\mathrm{erf}(\mathit{x}))\right),\end{array}$(11)

where the normalizing constant for the distribution is

(12) $\mathit{C}=\mathrm{\Re }\left[\frac{r_j\mathit{\pi }}{\mathit{i}{\beta }^{-1}}\mathit{w}(\mathit{z})\right]+\frac{\sqrt{\mathit{\pi }}}{2{\beta }^2}\left({\sigma }_0(1+2y^2)+{\sigma }_1\mathit{y}\right),$(12)

which we point out is nothing more than the Doppler-broadened scattering cross section at temperature $\mathit{T}$ under the single pole approximation.

The new special function we deem the incomplete Faddeeva function is defined as

(13) $\mathit{w}(\mathit{z},\mathit{x})=\frac{i}{\pi }{\int }_{-\mathrm{\infty }}^x\frac{e^{-t^2}}{\mathit{z}-\mathit{t}}\mathrm{d}\mathit{t}.$(13)

It can easily be seen that $\mathit{w}(\mathit{z},\mathrm{\infty })=\mathit{w}(\mathit{z})$, as per the definition of the Faddeeva function for $\mathrm{\Im }[\mathit{z}]>0$:

(14) $\mathit{w}(\mathit{z})=\frac{i}{\pi }{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\frac{e^{-t^2}}{\mathit{z}-\mathit{t}}\mathrm{d}\mathit{t}.$(14)

Indeed, for this application, $\mathrm{\Im }[\mathit{z}]>0$ and we maintain this assumption going forward. A specialized root finder has been developed to quickly invert this CDF and consequently sample the relative velocity. This thus constitutes a method to sample target velocities without rejection sampling or extensive tables.

The novelty in our approach lies entirely in the treatment of the zero-kelvin cross section and analytical representation of the relative speed CDF. After sampling from the relative speed distribution, we must sample the target polar angle distribution conditioned on the relative speed, as done in Ref. [6]. For completeness, we conclude with the CDF of the target speed

(15) $\mathit{C}(\mathit{V}|v_r)\propto \left\{\begin{array}{cc}0& \mathit{V}\le |v_r-\mathit{v}|\\ 1-e^{-{\beta }^2{V}^2}& |v_r-\mathit{v}|<\mathit{V}<v_r+\mathit{v}\\ 1& v_r+\mathit{v}\le \mathit{V}\end{array}\right.$(15)

for which Ref. [6] provides a straightforward sampling technique. The remaining work is purely numerical, particularly in requiring an efficient, reasonably accurate algorithm for the incomplete Faddeeva function $\mathit{w}(\mathit{z},\mathit{x})$.

II.A. Incomplete Faddeeva Function

The forthcoming discussion explores the properties of the incomplete Faddeeva function, with a particular focus on properties that can be leveraged to obtain efficient numerical approximations to it. We advise the reader that absorbing this section in depth is not necessary to grasp the claims and algorithms made herein.

The incomplete Faddeeva function, as defined by Eq. (13)(13) $\mathit{w}(\mathit{z},\mathit{x})=\frac{i}{\pi }{\int }_{-\mathrm{\infty }}^x\frac{e^{-t^2}}{\mathit{z}-\mathit{t}}\mathrm{d}\mathit{t}.$(13) is $\mathit{w}:\mathbb{C}\times \mathbb{R}\to \mathbb{C}$. This section attempts to build some intuition as to how this function behaves as $\mathit{z}$ and $\mathit{x}$ individually vary. In resonance upscattering treatment, $\mathit{z}$ parametrizes the location, height, and width of the resonance with respect to the Maxwell-Boltzmann distributed velocities. Values of $\mathrm{\Re }[\mathit{z}]=0$ correspond to scattering resonances exactly situated at the mean component of relative speed along the neutron’s line of flight predicted by the Maxwell-Boltzmann distribution. Values of $\mathrm{\Re }[\mathit{z}]>0$ correspond to resonances at higher energies than the particle’s energy, and therefore induce preferential scattering with relative velocities higher than the incident velocity, and vice-versa for $\mathrm{\Re }[\mathit{z}]<0$. $\mathit{x}$ parametrizes the dimensionless target relative speed. Small values of $\mathrm{\Im }[\mathit{z}]$ imply tall, narrow resonances, with a limit of $\mathrm{\Im }[\mathit{z}]=0$ being a singularity representing a resonance of infinite cross section. Large values of $\mathrm{\Im }[\mathit{z}]$ model wider, weaker resonances.

Figure 1 plots the incomplete Faddeeva function as a function of $\mathit{z}$ for a few values of $\mathit{x}$, and Fig. 2 shows the behavior of its real part in $\mathit{x}$ for a few values of $\mathit{z}$. Both illustrate the rapidly varying behavior of this function for small values of $\mathrm{\Im }[\mathit{z}]$ when $\mathrm{\Re }[\mathit{z}]\approx \mathit{x}$, where a sharp peak follows the value of $\mathit{x}$ near the real line. Figure 1 particularly shows the approach to the familiarly shaped $\mathit{w}(\mathit{z})$ as $\mathit{x}\to \mathrm{\infty }$. The figures also show that $\mathrm{\Re }[\mathit{w}(\mathit{z},\mathit{x})]$ is an increasing function in $\mathit{x}$, as would be expected of a CDF. This behavior is explained by the following asymptotic analysis.

Fig. 1. $\mathit{w}(\mathit{z},\mathit{x})$ for a few values of $\mathit{x}$. The height represents the magnitude and the coloring is done by phase.

Fig. 2. $\mathrm{\Re }[\mathit{w}(\mathit{z},\mathit{x})]$ for a few values of $\mathit{z}$. The legend is the imaginary number added to the real part specified in each figure’s caption. The plotted quantity is normalized by $\mathrm{\Re }[\mathit{w}(\mathit{z})]$ so all lines tend to unity as $\mathit{x}$ grows.

Humliĉek^[23] provides the following asymptotic formula for the Faddeeva function without derivation:

(16) $\mathit{w}(\mathit{z})\approx \frac{i}{\sqrt{\mathit{\pi }}\mathit{z}}.$(16)

Considering the definition of the Faddeeva function once more in Eq. (14)(14) $\mathit{w}(\mathit{z})=\frac{i}{\pi }{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\frac{e^{-t^2}}{\mathit{z}-\mathit{t}}\mathrm{d}\mathit{t}.$(14) , this approximation results from supposing that if $\mathit{z}$ is large, and that only values of $\mathit{t}$ near zero contribute to the integral, we can approximate $\mathit{w}(\mathit{z})$ as

(17) $\mathit{w}(\mathit{z})\approx \frac{i}{\pi }{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\frac{e^{-t^2}}{\mathit{z}}\mathrm{d}\mathit{t},$(17)

which immediately yields the asymptotic estimate of Eq. (16)(16) $\mathit{w}(\mathit{z})\approx \frac{i}{\sqrt{\mathit{\pi }}\mathit{z}}.$(16) . Proceeding with the same approximation in the context of the incomplete Faddeeva function, we thus obtain

(18) $\mathit{w}(\mathit{z},\mathit{x})\approx \frac{i}{2\sqrt{\mathit{\pi }}\mathit{z}}\left(\mathrm{erf}(\mathit{x})+1\right)\mathrm{if}|\mathit{z}|\gg 1,$(18)

suggesting a close connection between the incomplete Faddeeva function’s behavior in $\mathit{x}$ and the error function. In fact, this hunch is confirmed by the following identity, which connects the incomplete Faddeeva to the standard Faddeeva function:

(19) $\mathit{w}(\mathit{z},\mathit{x})=\frac{1}{2}(1+\mathrm{erf}(\mathit{x}))\mathit{w}(\mathit{z})+\frac{i{e}^{-x^2}}{\mathit{\pi }}{\int }_0^{\mathrm{\infty }}\frac{e^{-t^2}{e}^{2\mathit{itz}}\mathrm{d}\mathit{t}}{\mathit{i}(\mathit{x}-\mathit{z})+\mathit{t}}.$(19)

Proof of this relation is provided in Appendix A.

An even more accurate asymptotic estimate can be obtained for large $\mathrm{\Re }[\mathit{z}]$, approximating the pole term as

(20) $\frac{1}{\mathit{z}-\mathit{t}}\approx \frac{1}{2\mathit{z}}\left(1+e^{2\mathit{t}/z}\right).$(20)

This matches the value, slope, and curvature with respect to $\mathit{t}$ of the pole about $\mathit{t}=0$. Substituting this back to Eq. (13)(13) $\mathit{w}(\mathit{z},\mathit{x})=\frac{i}{\pi }{\int }_{-\mathrm{\infty }}^x\frac{e^{-t^2}}{\mathit{z}-\mathit{t}}\mathrm{d}\mathit{t}.$(13) and adjusting the expression such that $\mathrm{\Re }[\mathit{w}(\mathit{z},\mathit{x})]$ is strictly increasing [as suggested by Eq. (18)(18) $\mathit{w}(\mathit{z},\mathit{x})\approx \frac{i}{2\sqrt{\mathit{\pi }}\mathit{z}}\left(\mathrm{erf}(\mathit{x})+1\right)\mathrm{if}|\mathit{z}|\gg 1,$(18) ], and matching the asymptotic value of $\mathit{w}(\mathit{z})$ as suggested by Eq. (19)(19) $\mathit{w}(\mathit{z},\mathit{x})=\frac{1}{2}(1+\mathrm{erf}(\mathit{x}))\mathit{w}(\mathit{z})+\frac{i{e}^{-x^2}}{\mathit{\pi }}{\int }_0^{\mathrm{\infty }}\frac{e^{-t^2}{e}^{2\mathit{itz}}\mathrm{d}\mathit{t}}{\mathit{i}(\mathit{x}-\mathit{z})+\mathit{t}}.$(19) , results in

(21) $\mathit{w}(\mathit{z},\mathit{x})\approx \frac{1}{2}\left(1+\mathrm{erf}\left(\mathit{x}-\mathrm{\Re }{[\mathit{z}]}^{-1}\right)\right)\mathit{w}(\mathit{z}).$(21)

Equation (21)(21) $\mathit{w}(\mathit{z},\mathit{x})\approx \frac{1}{2}\left(1+\mathrm{erf}\left(\mathit{x}-\mathrm{\Re }{[\mathit{z}]}^{-1}\right)\right)\mathit{w}(\mathit{z}).$(21) is sufficiently accurate to be used in practical computations, as shown by Fig. 3.

Fig. 3. Accuracy of the $\mathit{w}(\mathit{z},\mathit{x})$ approximation for $|\mathrm{\Re }[\mathit{z}]|>5$. An unshifted error function is shown for comparison, representing a more naive asymptotic approximation to $\mathit{w}(\mathit{z},\mathit{x})$.

In the context of resonance scattering, Eq. (21)(21) $\mathit{w}(\mathit{z},\mathit{x})\approx \frac{1}{2}\left(1+\mathrm{erf}\left(\mathit{x}-\mathrm{\Re }{[\mathit{z}]}^{-1}\right)\right)\mathit{w}(\mathit{z}).$(21) shows that the relative speed distribution gets shifted forward by a nondimensionalized factor of $1/\mathrm{\Re }[\mathit{z}]$ with its influence scaling by $\mathrm{\Re }[r_j\mathit{w}(\mathit{z})]$ times the resonance’s residue, as suggested by Eq. (11)(11) $\begin{array}{rl}& \mathit{CP}(\mathit{x}|\mathit{y})=\mathrm{\Re }\left[\frac{r_j\mathit{\pi }}{\mathit{i}{\beta }^{-1}}\mathit{w}(\mathit{z},\mathit{x})\right]+\frac{{\mathit{\beta }}^{-2}{\sigma }_0}{4}(-2e^{-x^2}(\mathit{x}+2\mathit{y})\\ & +\sqrt{\mathit{\pi }}(1+2y^2)(1+\mathrm{erf}(\mathit{x})))\\ & +\frac{1}{2{\beta }^2}{\sigma }_1{e}^{-x^2}\left(1+{(\mathit{x}+\mathit{y})}^2+\sqrt{\mathit{\pi }}\mathit{y}(1+\mathrm{erf}(\mathit{x}))\right),\end{array}$(11) .

Another useful property of the incomplete Faddeeva function is a simple connection between its derivative in the complex plane and its value. This is similar in nature to the derivative of the Faddeeva function^[24]:

(22) $\frac{\mathrm{d}w(z)}{\mathrm{d}z}=\frac{2\mathit{i}}{\sqrt{\mathit{\pi }}}-2\mathit{zw}(\mathit{z}).$(22)

The relation we have obtained generalizes this as

(23) $\frac{\mathrm{d}w(z,x)}{\mathrm{d}z}=\frac{i}{\sqrt{\mathit{\pi }}}\left(1+\mathrm{erf}(\mathit{x})+\frac{e^{-x^2}}{\sqrt{\mathit{\pi }}(\mathit{x}-\mathit{z})}\right)-2\mathit{zw}(\mathit{z},\mathit{x}),$(23)

which clearly maintains consistency with Eq. (22)(22) $\frac{\mathrm{d}w(z)}{\mathrm{d}z}=\frac{2\mathit{i}}{\sqrt{\mathit{\pi }}}-2\mathit{zw}(\mathit{z}).$(22) as $\mathit{x}\to \mathrm{\infty }$. This fairly simple connection between the derivative of $\mathit{w}(\mathit{z},\mathit{x})$ and its value can be utilized for efficient sensitivity analyses of the scattering kernel with respect to WMP parameters.

The forthcoming discussion presents some further concepts in the direction of the efficient numerical evaluation of the incomplete Faddeeva function. Going forward, we denote the second integral appearing in Eq. (19)(19) $\mathit{w}(\mathit{z},\mathit{x})=\frac{1}{2}(1+\mathrm{erf}(\mathit{x}))\mathit{w}(\mathit{z})+\frac{i{e}^{-x^2}}{\mathit{\pi }}{\int }_0^{\mathrm{\infty }}\frac{e^{-t^2}{e}^{2\mathit{itz}}\mathrm{d}\mathit{t}}{\mathit{i}(\mathit{x}-\mathit{z})+\mathit{t}}.$(19) as

(24) $\mathit{I}(\mathit{z},\mathit{m})={\int }_0^{\mathrm{\infty }}\frac{e^{-t^2}{e}^{2\mathit{itz}}\mathrm{d}\mathit{t}}{\mathit{im}+\mathit{t}},$(24)

where $\mathbb{R}∍\mathit{m}=\mathit{x}-\mathit{z}$. While it may seem that a half-range Gauss-Hermite quadrature might work well to efficiently approximate Eq. (24(24) $\mathit{I}(\mathit{z},\mathit{m})={\int }_0^{\mathrm{\infty }}\frac{e^{-t^2}{e}^{2\mathit{itz}}\mathrm{d}\mathit{t}}{\mathit{im}+\mathit{t}},$(24) ), this is not the case. First, complex exponentials would have to be calculated at each quadrature point. Second, as the real part of $\mathit{z}$ grows, the integrand oscillates more. In practice, values of $\mathrm{\Re }[\mathit{z}]>10$ are frequently encountered, and low-degree quadratures would not capture the oscillation. Additionally, given that the value of $\mathrm{\Im }[\mathit{z}]$ is small, the denominator becomes nearly singular when $\mathit{x}\approx \mathrm{\Re }[\mathit{z}]$. In fact, when $\mathrm{\Im }[\mathit{z}]=0$, $\mathit{I}(\mathit{z},\mathit{m})$ becomes a discontinuous function in $\mathit{x}$ when interpreted as a principal value integral.

Equation (24)(24) $\mathit{I}(\mathit{z},\mathit{m})={\int }_0^{\mathrm{\infty }}\frac{e^{-t^2}{e}^{2\mathit{itz}}\mathrm{d}\mathit{t}}{\mathit{im}+\mathit{t}},$(24) is equivalent to the incomplete Goodwin-Staton integral, referred to in Ref. [25]. However, to our knowledge, only asymptotic analysis has been performed on this type of integral before, without any development of numerical routines. Recent work in the field of finance^[26] presents results for computing what the authors define as the extended incomplete Goodwin-Staton integral, for which the $\mathit{\nu }=1$ case is of interest in the present discussion. Unfortunately, the authors’ numerical method works for all cases except $\mathit{\nu }=1$, suggesting Eq. (24)(24) $\mathit{I}(\mathit{z},\mathit{m})={\int }_0^{\mathrm{\infty }}\frac{e^{-t^2}{e}^{2\mathit{itz}}\mathrm{d}\mathit{t}}{\mathit{im}+\mathit{t}},$(24) to be of a fundamentally different nature.

In our experience, the difficulty with large $\mathrm{\Re }[\mathit{z}]$ cannot be ameliorated by a stationary phase technique,^[27] as these tend to accentuate the pole behavior and remain of similar difficulty for half-range Gauss-Hermite quadratures.

In the WMP method, $\mathrm{\Im }[\mathit{z}]$ is near zero, as shown in Fig. 4. Therefore, we can expect to frequently encounter nearly singular integrands in Eq. (24)(24) $\mathit{I}(\mathit{z},\mathit{m})={\int }_0^{\mathrm{\infty }}\frac{e^{-t^2}{e}^{2\mathit{itz}}\mathrm{d}\mathit{t}}{\mathit{im}+\mathit{t}},$(24) . An efficient numerical technique that explicitly treats this behavior can be devised by first noticing that $\mathit{I}(\mathit{z},\mathit{m})$ satisfies this differential equation in the complex plane

(25) $\frac{\mathrm{d}I}{\mathrm{d}z}+2\mathit{zI}=\frac{1}{\mathit{x}-\mathit{z}}.$(25)

Fig. 4. Distribution of the imaginary part of the poles in the WMP method. These were collected from every nuclide in the OpenMC regression testing data set based on ENDFVII.1.

This can be used to connect the value $\mathit{I}(z_0,\mathit{x})$ at a point $z_0$ to another point $z_1$. The integrating factor technique shows that

(26) $\mathit{I}(z_1,\mathit{x})={\int }_{z_0}^{z_1}\frac{e^{t^2-z_1^2}\mathrm{d}\mathit{t}}{\mathit{x}-\mathit{t}}+e^{z_0^2-z_1^2}\mathit{I}(z_0,\mathit{x}).$(26)

If the distance between $z_0$ and $z_1$ is small, the exponential in the numerator of Eq. (26)(26) $\mathit{I}(z_1,\mathit{x})={\int }_{z_0}^{z_1}\frac{e^{t^2-z_1^2}\mathrm{d}\mathit{t}}{\mathit{x}-\mathit{t}}+e^{z_0^2-z_1^2}\mathit{I}(z_0,\mathit{x}).$(26) can be Taylor expanded about $z_0$ to yield an efficient numerical scheme.

We have found that the imaginary part of $\mathit{I}(\mathrm{\Re }[\mathit{z}],\mathit{x})$ can be calculated with a closed-form, discontinuous-in-$\mathit{x}$ formula. Because the nearly discontinuous behavior of $\mathit{I}(\mathit{z},\mathit{x})$ in $\mathit{x}$ is the main source of difficulty here, Eq. (26)(26) $\mathit{I}(z_1,\mathit{x})={\int }_{z_0}^{z_1}\frac{e^{t^2-z_1^2}\mathrm{d}\mathit{t}}{\mathit{x}-\mathit{t}}+e^{z_0^2-z_1^2}\mathit{I}(z_0,\mathit{x}).$(26) can be used to resolve this behavior accurately after calculating the value of $\mathit{I}(\mathrm{\Re }[\mathit{z}],\mathit{x})$. The real part of $\mathit{I}(\mathrm{\Re }[\mathit{z}],\mathit{x})$ is continuous in $\mathit{x}$ but is not available in terms of elementary functions; however, it is readily amenable to numerical approximation. In conclusion, the real and imaginary part of $\mathit{I}(\mathrm{\Re }[\mathit{z}],\mathit{x})$ contrast with each other: the former is readily numerically approximated by series or similar methods, whereas the latter is discontinuous and therefore not amenable to series or rational function approximation, and fortunately has an exact formula.

The first goal at hand is to calculate $\mathrm{\Im }[\mathit{I}(\mathit{z},\mathit{m})]$ for real values of $\mathit{z}$. This can be written as

(27) $\mathrm{\Im }[\mathit{I}(\mathit{z},\mathit{m})]={\int }_0^{\mathrm{\infty }}\frac{e^{-t^2}\left(\mathit{t}sin(2\mathit{zt})-\mathit{m}cos(2\mathit{zt})\right)\mathrm{d}\mathit{t}}{t^2+m^2}.$(27)

The linearity of integration can be distributed over both trigonometric terms. Each of the resulting integrals can be found as respective sine and cosine transform integrals, which are listed in Ref. [28]. Combining the results from the sine and cosine transform tables yields

(28) $\mathrm{\Im }[\mathit{I}(\mathit{z},\mathit{m})]=\mathit{\pi }{e}^{-z^2+{(\mathit{m}+\mathit{z})}^2}\left(\frac{1}{2}\mathrm{erf}(\mathit{m}+\mathit{z})-\mathrm{sign}(\mathit{m})\right)\mathit{z}\in \mathbb{R}.$(28)

Next, we can find an approximation for the real part of $\mathit{I}(\mathit{z},\mathit{m})$:

(29) $\mathrm{\Re }[\mathit{I}(\mathit{z},\mathit{m})]={\int }_0^{\mathrm{\infty }}\frac{e^{-t^2}\left(\mathit{t}cos(2\mathit{zt})+\mathit{m}sin(2\mathit{zt})\right)\mathrm{d}\mathit{t}}{t^2+m^2}.$(29)

This expression can neither be written in terms of elementary nor special functions to our knowledge. In order to manipulate it to obtain a numerical expression, consider the auxiliary function

(30) $\mathit{R}(\mathit{z},\mathit{m})={\int }_0^{\mathrm{\infty }}\frac{e^{-t^2}sin(2\mathit{zt})\mathrm{d}\mathit{t}}{t^2+m^2},$(30)

which again, of course, cannot be represented in terms of elementary or special functions. This pinpoints the difficulty in calculating the real part of $\mathit{I}(\mathit{z},\mathit{m})$ because

(31) $\mathrm{\Re }[\mathit{I}(\mathit{z},\mathit{m})]=\mathit{mR}(\mathit{z},\mathit{m})+\frac{1}{2}\frac{\mathrm{\partial R}}{\mathrm{\partial z}}.$(31)

We thus seek a straightforwardly differentiable approximation to $\mathit{R}(\mathit{z},\mathit{m})$. The intuition behind our forthcoming numerical approximation to $\mathit{R}(\mathit{m},\mathit{z})$ comes from the fact that for $\mathit{m}\gg 1$, $\mathit{t}$ is negligible in the denominator compared to $\mathit{m}$ over the range where the $e^{-t^2}$ weighting is large, hence

(32) $\mathit{R}(\mathit{m},\mathit{z})\approx \frac{1}{m^2}{\int }_0^{\mathrm{\infty }}{e}^{-t^2}sin(2\mathit{zt})\mathrm{d}\mathit{t}=\frac{1}{m^2}\mathit{F}(\mathit{z}),$(32)

Where $\mathit{F}(\mathit{z})$ is the Dawson F function.^[24] By standard asymptotic analysis, matching the $\mathit{z}$ derivative at $\mathit{z}=0$ and the $\mathit{z}\gg 1$ asymptote for the Dawson F function results in the improved estimate

(33) $\mathit{R}(\mathit{m},\mathit{z})\approx 1\text{overm}\sqrt{e^{m^2}{E}_1(m^2)}\mathit{F}\left(\mathit{mz}\sqrt{e^{m^2}{E}_1(m^2)}\right),$(33)

which is accurate to within about 5% across the full range of $\mathit{m}$ values. However, this level of accuracy is not appropriate for engineering calculations.

To build more intuition for $\mathit{R}(\mathit{m},\mathit{z})$, its close relation to the Dawson F function is evinced by considering the Maclaurin series in $\mathit{z}$ for both

(34) $\mathit{F}(\mathit{z})=\mathit{z}-\frac{2}{3}{z}^3+\frac{4}{15}{z}^5-\frac{8}{105}{z}^7+\cdots$(34)

and

(35) $\begin{array}{rl}& \mathit{R}(\mathit{m},\mathit{z})=e^{m^2}\left(E_1(m^2)\mathit{z}-\frac{2}{3}{E}_2(m^2)z^3\right.\\ & \left.+\frac{4}{15}{E}_3(m^2)z^5-\frac{8}{105}{E}_4(m^2)z^7+\cdots \right),\end{array}$(35)

which further highlight the utility of maintaining the consistency of $\mathit{R}(\mathit{m},\mathit{z})$ with its asymptotic sister $\mathit{F}(\mathit{z})$. The series are the same, but with the addition of the exponential integral multiplying factors in $\mathit{R}(\mathit{m},\mathit{z})$.

In order to find such a numerical relation that maintains this consistency in the asymptotic case, Eq. (30)(30) $\mathit{R}(\mathit{z},\mathit{m})={\int }_0^{\mathrm{\infty }}\frac{e^{-t^2}sin(2\mathit{zt})\mathrm{d}\mathit{t}}{t^2+m^2},$(30) can be transformed via the interchange of differentiation and integration. Consider the generalized function

(36) $\mathit{R}(\mathit{m},\mathit{z},\mathit{\alpha })=e^{m^2}{\int }_0^{\mathrm{\infty }}\frac{e^{-\alpha (t^2+m^2)}sin(2\mathit{zt})\mathrm{d}\mathit{t}}{t^2+m^2}.$(36)

Differentiating reveals that

(37) $\frac{\mathrm{\partial R}(\mathit{m},\mathit{z},\mathit{\alpha })}{\mathrm{\partial \alpha }}=-e^{m^2}{\int }_0^{\mathrm{\infty }}{e}^{-\alpha (t^2+m^2)}sin(2\mathit{zt})\mathrm{d}\mathit{t}=-\frac{e^{-m^2(\mathit{\alpha }-1)}}{\sqrt{\mathit{\alpha }}}\mathit{F}\left(\frac{z}{\sqrt{\mathit{\alpha }}}\right).$(37)

The fundamental theorem of calculus then applies:

(38) $\mathit{R}(\mathit{m},\mathit{z},\mathrm{\infty })-\mathit{R}(\mathit{m},\mathit{z},1)=e^{m^2}{\int }_1^{\mathrm{\infty }}\frac{1}{\sqrt{\mathit{\alpha }}}{e}^{-m^2\mathit{\alpha }}\mathit{F}\left(\frac{z}{\sqrt{\mathit{\alpha }}}\right)\mathrm{d}\mathit{\alpha }.$(38)

Using the fact that $\mathit{R}(\mathit{m},\mathit{z},\mathrm{\infty })=0$ and doing a change of variables, Eq. (30)(30) $\mathit{R}(\mathit{z},\mathit{m})={\int }_0^{\mathrm{\infty }}\frac{e^{-t^2}sin(2\mathit{zt})\mathrm{d}\mathit{t}}{t^2+m^2},$(30) becomes

(39) $\mathit{R}(\mathit{m},\mathit{z})=2e^{m^2}{\int }_1^{\mathrm{\infty }}{e}^{-m^2{t}^2}\mathit{F}(\mathit{z}/\mathit{t})\mathrm{d}\mathit{t},$(39)

which confirms our hunch about the close relation of $\mathit{F}(\mathit{z})$ to $\mathit{R}(\mathit{m},\mathit{z})$; it is an infinite superposition of stretched and scaled Dawson F functions.

We have found success in inserting an approximation for $\mathit{F}(\mathit{z})$ into Eq. (39)(39) $\mathit{R}(\mathit{m},\mathit{z})=2e^{m^2}{\int }_1^{\mathrm{\infty }}{e}^{-m^2{t}^2}\mathit{F}(\mathit{z}/\mathit{t})\mathrm{d}\mathit{t},$(39) . Well-known approximations to $\mathit{F}(\mathit{z})$ based on rational expressions and other elementary functions^[29,30] do not result in numerically useful expressions. However, noting that

(40) ${\int }_1^{\mathrm{\infty }}{e}^{-m^2{t}^2}(\mathit{z}/\mathit{t}{)}^{(2\mathit{k}+1)}\mathrm{d}\mathit{t}=\frac{1}{2}{z}^{1+2\mathit{k}}{E}_{1+\mathit{k}}(m^2),$(40)

it can be seen that approximations to $\mathit{F}(\mathit{z})$ in the form of a power series can be computationally efficient in light of the recursion relation for exponential integrals:

(41) $\mathit{n}{e}^{m^2}{E}_{\mathit{n}+1}(m^2)=(1-m^2{e}^{m^2}{E}_n(m^2)).$(41)

Careful attention must be paid to the floating point properties of this relation.^[31] The magnification of computational errors grows arbitrarily large, and we later present a specialized numerical algorithm guaranteeing floating point stability.

In order to thus obtain a simple, efficient approximation to $\mathit{R}(\mathit{m},\mathit{z})$, we employ the Chebyshev expansion valid for $\mathit{z}\in [-5,5]$ presented in Ref. [32]. The results of our method could be improved by using a finer piecewise division for the Chebyshev expansion of $\mathit{F}(\mathit{z})$ as in Ref. [33], but we have used the present approach for simplicity of implementation. Thus, if the truncated Chebyshev expansion of $\mathit{F}(\mathit{z})$ is converted to the power series basis

(42) $\mathit{F}(\mathit{z})\approx \sum _{\mathit{i}=0}^n{c}_n(\mathit{z}/5{)}^{2\ast \mathit{i}+1},$(42)

we obtain approximations of the form

(43) $\mathit{R}(\mathit{m},\mathit{z})\approx \sum _{\mathit{i}=10}^n\frac{c_n}{2}{E}_{1+\mathit{i}}(m^2)(\mathit{z}/5{)}^{2\ast \mathit{i}+1}.$(43)

III. NUMERICAL IMPLEMENTATION

III.A. Stable, Efficient Calculation of an $E_n$ Sequence

The magnification of error in the forward recurrence relation for exponential integrals from Ref. [31] is

(44) $|{\rho }_n|=\frac{x^n{E}_1(\mathit{x})}{\mathit{n}!E_{\mathit{n}+1}(\mathit{x})}.$(44)

Notably, the error magnification of the reverse recurrence relation is the reciprocal of this quantity. Moreover, $|{\rho }_n|$ is a function increasing from 1, reaching a maximum and monotonically descending below 1.^[31] As a consequence, a critical index $n^{\ast }$ exists such that iterating outward from it results in a numerically stable recursion algorithm. In terms of evaluating Eq. (43)(43) $\mathit{R}(\mathit{m},\mathit{z})\approx \sum _{\mathit{i}=10}^n\frac{c_n}{2}{E}_{1+\mathit{i}}(m^2)(\mathit{z}/5{)}^{2\ast \mathit{i}+1}.$(43) , this means splitting the polynomial in $\mathit{z}$ into parts above the index $n^{\ast }$ and those below. After computing $e^{m^2}{E}_{n^{\ast }}(m^2)$, Horner’s method is used in the reverse recurring relation down to the term of order $\mathit{z}$, and forward recursion is employed to evaluate the polynomial of degree leading from $2n^{\ast }+1$ up to $2\mathit{n}+1$.

A simple result we have obtained is that the smallest value of $n^{\ast }$, such that

(45) $\frac{x^{n^{\ast }}{E}_1(\mathit{x})}{n^{\ast }!E_{n^{\ast }+1}(\mathit{x})}<1,$(45)

is well approximated by

(46) $n^{\ast }\approx \mathit{ex}-\frac{1}{2}log\mathit{\pi }.$(46)

Appendix B shows how this can be obtained. Figure 5 illustrates the accuracy of Eq. (46)(46) $n^{\ast }\approx \mathit{ex}-\frac{1}{2}log\mathit{\pi }.$(46) . Using this information, Algorithm 1 explicitly states the procedure to calculate $\mathit{R}(\mathit{m},\mathit{z})$ and its derivative with respect to $\mathit{z}$. While the use of a power-basis polynomial is suboptimal, numerically speaking, the main source of numerical error in this scenario originates from the recursive exponential integral formula. The specification of the algorithm assumes that an accurate method for computing $E_n(\mathit{x})e^x$ has been provided, which is well documented in many other works. We have employed a C++ adaptation of the continued fraction approximation employed by the Cephes library,^[34] which is documented in Ref. [24].

Fig. 5. Approximate solution to finding the first value of $\mathit{n}$ such that $\frac{x^n{E}_1(\mathit{x})}{\mathit{n}!E_{\mathit{n}+1}(\mathit{x})}<1$.

Table

Algorithm 1: Stable $\mathit{R}(\mathit{m},\mathit{z})$ Approximation for $|\mathit{z}|\le 5$. Appendix C is used for $|\mathit{z}|>5$.

Input : $\mathit{m}\in \mathbb{R},\mathit{z}\in \mathbb{R}$

Output: R(m,z) from Eq. 30(30) $\mathit{R}(\mathit{z},\mathit{m})={\int }_0^{\mathrm{\infty }}\frac{e^{-t^2}sin(2\mathit{zt})\mathrm{d}\mathit{t}}{t^2+m^2},$(30) and $\frac{\mathrm{\partial }\mathit{R}}{\mathrm{\partial z}}$

Data: cn coefficients of Eq. 42(42) $\mathit{F}(\mathit{z})\approx \sum _{\mathit{i}=0}^n{c}_n(\mathit{z}/5{)}^{2\ast \mathit{i}+1},$(42) , 1≤ n ≤ nmax

n* ← min(max(em^2 − 0.57, 1), nmax);

expnexp ← En* (m^2)e^m2;

expnexp_orig ← expnexp;

result ← 0;

derivative ← 0;

// Backward recurse with Horner scheme for n ← n* – 1 to 1 do

expnexp ← (1– nexpnexp) /m^2;

derivative ← derivative · pow(z/5, 2)+ result;

result ← (result)pow (z/5, 2) + expnexp˙ Cn–1 ;

end

derivative ← 2derivative ·pow(z/5, 2) ;

derivative ← derivative + result ;

result ← result · z/5 ;

// Forward recursing polynomial evaluation

x2 ← pow(z/5, 2n* – 1);

X1← pow(z/5, 2(n* – l));

expnexp ← expnexp_orig;

for n ← n* to nmax do

result ← cn–1·expnexp · x2;

derivative ← cn–1 ·expnexp · xl ·(2n – 1);

x2 ← x2 · pow(z/5, 2);

xl ← xl · pow(z/5, 2);

expnexp ← (1 – m^2 ·expnexp)/n;

end

return (result, derivative);

III.A.1. Jump Integral

After computing the value of $\mathit{I}(\mathrm{\Re }[\mathit{z}],\mathit{m})$, the differential equation of Eq. (25)(25) $\frac{\mathrm{d}I}{\mathrm{d}z}+2\mathit{zI}=\frac{1}{\mathit{x}-\mathit{z}}.$(25) , which $\mathit{I}$ follows, can be used to calculate $\mathit{I}(\mathit{z},\mathit{m})$. We calculate $\mathit{I}(\mathit{z},\mathit{m})$ in this manner due to the nearly discontinuous behavior of $\mathit{I}(\mathit{z},\mathit{m})$; it has a jump discontinuity about $\mathit{m}=0$ if $\mathit{z}\in \mathbb{R}$. Because the imaginary part of $\mathit{z}$ is small in WMP libraries, the resulting behavior is nearly discontinuous, and hence is not captured efficiently by general approximation techniques; finely resolved tables or high polynomial orders would be required. Our approach thus resolves the discontinuous component exactly with the piecewise function [Eq. (28)(28) $\mathrm{\Im }[\mathit{I}(\mathit{z},\mathit{m})]=\mathit{\pi }{e}^{-z^2+{(\mathit{m}+\mathit{z})}^2}\left(\frac{1}{2}\mathrm{erf}(\mathit{m}+\mathit{z})-\mathrm{sign}(\mathit{m})\right)\mathit{z}\in \mathbb{R}.$(28) ]. The nontrivial part of Eq. (26)(26) $\mathit{I}(z_1,\mathit{x})={\int }_{z_0}^{z_1}\frac{e^{t^2-z_1^2}\mathrm{d}\mathit{t}}{\mathit{x}-\mathit{t}}+e^{z_0^2-z_1^2}\mathit{I}(z_0,\mathit{x}).$(26) is the transcendental integral

(47) $\mathit{J}(\mathit{z},\mathit{x})=e^{-\mathrm{\Re }{[\mathit{z}]}^2}{\int }_{\mathrm{\Re }[\mathit{z}]}^z\frac{e^{t^2}\mathrm{d}\mathit{t}}{\mathit{x}-\mathit{t}}.$(47)

We have deemed this term the jump integral because it allows for jumping from values of $\mathit{I}(\mathit{z},\mathit{m})$ on the real line to values above the real line in the complex plane. While it seems that our issue of approximating the transcendental integral $\mathit{w}(\mathit{z},\mathit{x})$ has seemingly not been heretofore ameliorated due to the appearance of yet another transcendental integral [Eq. (47)(47) $\mathit{J}(\mathit{z},\mathit{x})=e^{-\mathrm{\Re }{[\mathit{z}]}^2}{\int }_{\mathrm{\Re }[\mathit{z}]}^z\frac{e^{t^2}\mathrm{d}\mathit{t}}{\mathit{x}-\mathit{t}}.$(47) ], a change of variables puts it into a form suitable for numerical approximation,

(48) $\mathit{J}(\mathit{z},\mathit{x})={\int }_0^{\mathit{i}\mathrm{\Im }[\mathit{z}]}\frac{e^{u^2+2\mathit{u}\mathrm{\Re }[\mathit{z}]}\mathrm{d}\mathit{u}}{\mathit{m}-\mathit{u}},$(48)

where again, $\mathit{m}=\mathrm{\Re }[\mathit{z}]-\mathit{x}$. Because $\mathrm{\Im }[\mathit{z}]$ is small as, shown by Fig. 4, the argument to the exponential term is similarly small. Where this integral is well defined ($\mathit{m}\ne 0$), the exponential term can be expanded in its Maclaurin series and integrated term by term

(49) $e^{u^2+2\mathit{u}\mathrm{\Re }[\mathit{z}]}=\sum _{\mathit{n}=0}^{\mathrm{\infty }}\frac{a_n}{\mathit{n}!}{u}^n.$(49)

It is verified that the coefficients $a_n$ satisfy the two-term recurrence

(50) $a_{\mathit{n}+1}=2\mathrm{\Re }[\mathit{z}]a_n+2(\mathit{n}-1)a_{\mathit{n}-1};a_0=1;a_1=2\mathrm{\Re }[\mathit{z}].$(50)

Next, the term-by-term integrals appear in the form

(51) ${\int }_0^{\mathit{i}\mathrm{\Im }[\mathit{z}]}\frac{u^n\mathrm{d}\mathit{u}}{\mathit{m}-\mathit{u}}=m^n{B}_{\mathit{i}\mathrm{\Im }[\mathit{z}]/m}(1+\mathit{n},0),$(51)

where $B_x(\cdot ,\cdot )$ is the incomplete beta function, defined as

(52) $B_x(\mathit{a},\mathit{b})={\int }_0^z{t}^{\mathit{a}-1}(1-\mathit{t}{)}^{\mathit{b}-1}\mathrm{d}\mathit{t}.$(52)

A recursion formula derived as a special case of formulas in Ref. [24] efficiently calculates these incomplete beta function values of higher $\mathit{n}$ in sequence:

(53) $B_x(\mathit{n}+1,0)=B_x(\mathit{n},0)-\frac{x^n}{\mathit{n}}.$(53)

In combination with the fact that

(54) $B_x(1,0)=-log(1-\mathit{x}),$(54)

this yields an efficient numerical scheme for evaluating an integral of the truncated Maclaurin series of the exponential of Eq. (48)(48) $\mathit{J}(\mathit{z},\mathit{x})={\int }_0^{\mathit{i}\mathrm{\Im }[\mathit{z}]}\frac{e^{u^2+2\mathit{u}\mathrm{\Re }[\mathit{z}]}\mathrm{d}\mathit{u}}{\mathit{m}-\mathit{u}},$(48) .

Algorithm 2 details the combination of all of these facts for an efficient approximation to $\mathit{J}(\mathit{m},\mathit{z})$. This approximation works very well for problems with $|\mathrm{\Re }[\mathit{z}]|\le 5$, which easily covers the range of scattering events where resonances appreciably affect the double differential at temperature. Outside of that range, the integral becomes increasingly oscillatory, so an asymptotic approximation is employed for $|\mathit{z}|>5$. This approximation is documented in Appendix D.

Table

Algorithm 2: Efficient $\mathit{J}(\mathit{m},\mathit{z})$ Approximation for $\mathrm{\Im }[\mathit{z}]<1$ and $|\mathrm{\Re }[\mathit{z}]|<5$. Appendix D describes the approximation for $|\mathrm{\Re }[\mathit{z}]|>5$.

Input : $\mathit{x}\in \mathbb{R},\mathit{z}\in \mathbb{C}$

Output: J(m,z) from Eq. 47(47) $\mathit{J}(\mathit{z},\mathit{x})=e^{-\mathrm{\Re }{[\mathit{z}]}^2}{\int }_{\mathrm{\Re }[\mathit{z}]}^z\frac{e^{t^2}\mathrm{d}\mathit{t}}{\mathit{x}-\mathit{t}}.$(47) result ← 0;

aO ← L0;

al ← 2$\mathrm{\Re }$[z];

result ← bl · a0;

bl ← bl – i$\mathrm{\Im }$[z]/ (x – $\mathrm{\Re }$[x]);

// nmax adjusts the number of truncated series terms for n ← 2 to nmax do

tmp ← 2$\mathrm{\Re }$[z] · al + 2(n – 1) · a0;

a0 ← a1;

al ← tmp;

result ← result + bl ·al $\frac{m^{\mathit{n}}}{\mathit{n}{!}_{\mathit{n}}};$

b1 ← b1 $-\frac{1}{\mathit{n}+1}\left(\frac{\mathit{i}\mathrm{\Im }[\mathit{z}]}{(\mathit{x}-\mathrm{\Re }[\mathit{x}])}\right);$

end

ret1u-n result;

Last, Algorithm 3 gives the overall algorithm to compute $\mathit{w}(\mathit{z},\mathit{x})$ efficiently. It relies on access to some implementation of calculating $\mathit{w}(\mathit{z})$, e.g., the Johnson’s permissively licensed library,^[35] which implements a variety of approximations to achieve high accuracy, or one of the various rational approximations^[23,36,37] when higher error is permitted. Regardless of the chosen $\mathit{w}(\mathit{z})$ implementation, our algorithm maintains asymptotic consistency such that ${lim}_{\mathit{x}\to \mathrm{\infty }}\mathit{w}(\mathit{z},\mathit{x})=\mathit{w}(\mathit{z})$. This work leverages a recent approximation tailored for WMP.^[38]

Table

Algorithm 3: Efficient $\mathit{w}(\mathit{z},\mathit{x})$ Approximation for $\mathrm{\Im }[\mathit{z}]<<1$.

Input : $\mathit{z}\in \mathbb{C},\mathit{x}\in \mathbb{R}$ Output: w(z x) from Eq. 13(13) $\mathit{w}(\mathit{z},\mathit{x})=\frac{i}{\pi }{\int }_{-\mathrm{\infty }}^x\frac{e^{-t^2}}{\mathit{z}-\mathit{t}}\mathrm{d}\mathit{t}.$(13) m ← x – $\mathrm{\Re }$[z];

I I R(m,z) integral and its derivative

rmz, clrmzclz - call(Alg. l);

// l(m,z) integral

imz ← m · rmz + $\frac{1}{2}$drmzdz +$\frac{\mathit{i\pi }}{2}{e}^{-\mathrm{\Re }[\mathit{z}{]}^2+x^2}$(erf(x) – sign(m));

// J (m,z) integral

ji ← call(Alg. 2); imz ← irnz+ ji;

imz ← imz·$e^{\mathrm{\Re }[\mathit{z}{]}^2-z^2}$,

result ← $\frac{i}{\pi }{e}^{-{\mathit{x}}^2}$. imz;

result ← result +$\frac{1}{2}(\mathrm{erf}(\mathit{x})+1)\mathit{w}(\mathit{z});$

return result;

III.B. Pole Sampling Approximation

A key approximation of our technique that enables its computational efficiency is viewing the multipole cross section in the relative speed distribution as a mixture distribution. The theoretical justification is that if poles are present and sufficiently close to the incident neutron energy ($|\mathit{z}|<20$ specifically), the relative speed probability density function (PDF) is well approximated by ignoring the polynomial contribution

(55) $\mathit{P}(\mathit{x}|\mathit{y})\approx e^{-x^2}\mathrm{\Re }\left[\sum _{\mathit{j}\in \mathit{W}({\mathit{\beta }}^{-2}{y}^2)}\frac{\beta r_j}{z_j-\mathit{x}}\right].$(55)

This expression is not employed to actually sample the scattering distribution. Rather, it is best viewed as a mixture distribution in which each pole contributes a probability proportional to

(56) $\mathbb{P}[{\sigma }_s(\mathit{x})=\frac{\beta r_j}{\mathit{z}-\mathit{x}}+{\sigma }_{0,\mathit{j}}+{\sigma }_{1,\mathit{j}}\mathit{x}]\propto \mathrm{\Re }\left[r_j\mathit{w}(z_j)\right],$(56)

which defines a discrete distribution. In order to avoid the need for auxiliary storage and the calculation of a normalizing constant to this distribution, we recommend finding the maximum of $-|r_j\mathit{w}(z_j)|/log({\xi }_j)$, where ${\xi }_j$ are uniform random numbers differing for each pole. The $\mathit{j}$ corresponding to the maximum of this expression follows the desired discrete distribution. We also note that the quantity $\mathrm{\Re }\left[r_j\mathit{w}(z_j)\right]$ is exactly the Doppler-broadened contribution to the integrated scattering cross section, so this sampling procedure incurs no additional Faddeeva function evaluation overhead if this is done in tandem with a WMP cross-section lookup operation.

Finally, we emphasize that the pole sampling approximation is not precisely consistent with the original multipole cross-section representation. Instead, it uses the fact that polynomial contributions to the cross section negligibly affect the scattering kernel, while poles do so substantially.

III.C. Finding the Values of ${\sigma }_0$ and${\sigma }_1$

While it may seem that the polynomial contribution to Eq. (5)(5) $\mathit{\sigma }(\mathit{E})=\frac{1}{E}\mathrm{\Re }\left[\frac{r_j}{p_j-\sqrt{\mathit{E}}}\right]+{\sigma }_0+{\sigma }_1\sqrt{\mathit{E}}.$(5) could come as the first terms from the polynomials defined within the windows, as Ref. [19] used, we have found in practice that this choice is inconsistent with the approximation of the pole sampling technique and leads to negative cross-section values.

To remedy this issue, we use the heuristic that the relative error of the cross section local approximation is minimized by matching polynomial values in the vicinity of the dip of the resonance. The location of the scattering resonance trough is calculated as

(57) ${\sqrt{E}}_{\mathrm{trough}}=\frac{-b+\sqrt{b^2-\mathit{abc}+a^2\mathit{d}}}{\mathit{a}},$(57)

where

(58) $\mathit{a}=-\mathrm{\Im }\left[\left(r_{{\mathit{j}}^{\ast }}-{\bar{r}}_{\mathit{j}}^{\ast }\right)\right],$(58)

(59) $\mathit{b}=-\mathrm{\Im }\left[\left({\bar{r}}_{\mathit{j}}^{\ast }{p}_j^{\ast }-r_j^{\ast }{\bar{p}}_{\mathit{j}}^{\ast }\right)\right],$(59)

(60) $\mathit{c}=-\mathrm{\Re }\left[\left({\bar{p}}_{\mathit{j}}^{\ast }+p_j^{\ast }\right)\right],$(60)

(61) $\mathit{d}=|p_j^{\ast }{|}^2.$(61)

At this point, the window index of ${\sqrt{E}}_{\mathrm{trough}}$ is calculated.^a This may be a different window from the incident neutron energy’s window. The WMP cross section of Eq. (4)(4) $\mathit{\sigma }(\mathit{E})=\frac{1}{E}\left(\mathrm{\Re }\left[\sum _j\frac{r_j}{p_j-\sqrt{\mathit{E}}}\right]+\sum _{\mathit{n}=0}^N{a}_n{E}^{\mathit{n}/2}\right).$(4) is then evaluated at ${\sqrt{E}}_{\mathrm{trough}}$ but excluding the sampled pole $p_j$, i.e.,

(62) ${\sigma }_0=\frac{1}{{\mathit{E}}_{\mathrm{trough}}}\mathrm{\Re }\left[\sum _{j^{\ast }\ne \mathit{j}\in \mathit{W}(E_{\mathrm{trough}})}\frac{r_j}{p_j-\sqrt{E_{\mathrm{trough}}}}\right]+\sum _{\mathit{n}=0}^N{a}_n{E}_{\mathrm{trough}}^{\mathit{n}/2}.$(62)

From there, ${\sigma }_1\approx \frac{\mathrm{\partial }{\sigma }_s(\mathit{E})}{\mathrm{\partial }\sqrt{\mathit{E}}}$ is calculated at the same point, this time only including contributions from the polynomial expansion but not from any poles. This linearization technique has been found to improve the accuracy of our method when considering nuclides with tightly spaced resonances such as ²³⁵U. For complete clarity, the resulting expression is

(63) ${\sigma }_1=\sum _{\mathit{n}=0}^N{a}_n\frac{\mathit{n}-2}{2}{E}_{\mathrm{trough}}^{\mathit{n}/2-2}.$(63)

Finally, a linearization of the cross section in $\sqrt{\mathit{E}}$ space has been obtained. For use with the root finder, the nondimensional variable $\mathit{\beta }(\sqrt{\mathit{E}}-\sqrt{E_{\mathrm{incident}}})$ is preferable, so ${\sigma }_0$ is appropriately shifted and ${\sigma }_1$ appropriately scaled.

III.D. Inverting the Relative Speed CDF

To sample from the relative speed PDF, we employ the CDF inversion technique. A naive attempt at this would be a few bisection root-finding steps followed by a handful of Newton-like iterations. In practice, we have found that five bisection iterations followed by three Halley-Newton iterations resolves the root to within an acceptable tolerance; however, a far more efficient root finder has been developed that takes a maximum of four iterations total, only requiring more work for unusual edge cases.

The bootstrapping step, as we call it, is essential to an efficient implementation of MARS. The bootstrapping step cheaply obtains an initial guess to the solution of the CDF inversion problem, from which a small number of Newton-like iterations improve the solution.

The key to doing so lies in finding a cheap approximation to the inverse of the CDF with general pole parameters. In order to do so, we first move from the root-finding space of $\mathit{x}\in (-\mathrm{\infty },\mathrm{\infty })$ to the nonlinearly mapped variable $\tilde{x}=\frac{1}{2}\left(1+\mathrm{erf}\mathit{x}\right)$. The intuition behind using this modified space is that as the resonances become weak and the incident neutron energy becomes high, it can be shown that the CDF is simply equal to $\tilde{x}$, which ranges between zero and one. Resonances and low-energy free gas effects simply act as perturbations to this linear function, which enables a good starting point for approximating the root location.

The next step in improving the CDF model in the mapped space is to observe that the contribution of collision probability from the resonance largely does not depend on its imaginary part. Increasing the imaginary part of the resonance broadens it and decreases its height. Therefore, the magnitude of the jump in $\mathit{w}(\mathit{z},\mathit{x})$ when $\mathit{x}$ is near $\mathrm{\Re }[\mathit{z}]$ quantifies the probability that the neutron experiences a collision near the peak of the resonance. The jump in $\mathit{w}(\mathit{z},\mathit{x})$ for small $\mathrm{\Im }[\mathit{z}]$ about $\mathrm{\Re }[\mathit{z}]$ is approximately $e^{-\mathrm{\Re }{[\mathit{z}]}^2}$, and therefore the probability contribution due to the resonance is approximately

(64) $p_{\mathrm{jump}}=\mathbb{P}[\mathit{x}\approx \mathrm{\Re }[\mathit{z}]]\approx \mathit{\beta }{r}_{j^{\ast }}\mathit{\pi }{e}^{-\mathrm{\Re }{[\mathit{z}]}^2}/\mathit{C},$(64)

where $\mathit{C}$ is the normalizing constant given by Eq. (12)(12) $\mathit{C}=\mathrm{\Re }\left[\frac{r_j\mathit{\pi }}{\mathit{i}{\beta }^{-1}}\mathit{w}(\mathit{z})\right]+\frac{\sqrt{\mathit{\pi }}}{2{\beta }^2}\left({\sigma }_0(1+2y^2)+{\sigma }_1\mathit{y}\right),$(12) . Because this is an approximation, the probability of Eq. (65)(65) $\mathit{w}(\mathit{z},0)=\frac{1}{2}\mathit{w}(\mathit{z})+\frac{i}{2\mathit{\pi }}{e}^{-z^2}{E}_1(-z^2).$(65) may not be bounded between zero and one, so we threshold it to that range. In practice, the estimate provided here is accurate. We have found that this probability tends to be added into the CDF about $\mathrm{\Re }[\mathit{z}]$ over the interval $[\mathrm{\Re }[\mathit{z}]-\frac{3}{2}\mathrm{\Im }[\mathit{z}],\mathrm{\Re }[\mathit{z}]+\frac{3}{2}\mathrm{\Im }[\mathit{z}]]$. This estimate could obviously be tuned for greater accuracy.

One final tool we employ to bootstrap the root-finding process pertains to the values of the CDF about $\mathit{x}=0$. In this case, numerous instances of functions occurring in its expression, such as $\mathit{erf}(\mathit{x})$ and $e^{-x^2}$, take on easily calculated values. On top of that, the incomplete Faddeeva function has a closed-form expression when $\mathit{x}=0$ of

(65) $\mathit{w}(\mathit{z},0)=\frac{1}{2}\mathit{w}(\mathit{z})+\frac{i}{2\mathit{\pi }}{e}^{-z^2}{E}_1(-z^2).$(65)

Because $e^{-z^2}$ is already computed and cached for the CDF inversion, the calculation of a complex exponential integral is the only difficulty. This is much easier and computationally cheaper to do than the more involved $\mathit{w}(\mathit{z},\mathit{x})$ evaluation, so any off-the-shelf approximation of $E_1(-z^2)$ can be employed here.

Additionally, the derivatives of the CDF with respect to $\mathit{x}$ about $\mathit{x}=0$ are also easily obtainable, which we use to further improve our root-finding guess. So far, we have only incorporated information from the first derivative, which has proven sufficient.

This leaves us with the following pieces of information from which the root estimate is extracted: the probability due to the resonance, its width, the value and slope of the CDF about $\mathit{x}=0$, i.e., $\tilde{x}=1/2$, and the known endpoint values of the CDF at 0 and 1. We therefore construct a function that is piecewise quadratic on the left and right of $\tilde{x}=1/2$. This quadratic interval ranges to either the endpoints $\tilde{x}=0$ or $\tilde{x}=1$, or the resonance’s upper or lower range of probability gain, estimated here as x$\in [\mathrm{\Re }[\mathit{z}]-\frac{3}{2}\mathrm{\Im }[\mathit{z}],\mathrm{\Re }[\mathit{z}]+\frac{3}{2}\mathrm{\Im }[\mathit{z}]]$. Note that this interval has to be mapped to an interval in $\tilde{x}$ space. Because the interval of the resonance is small, the Jacobian of the transformation $\mathit{x}\mapsto \tilde{x}$, which is proportional to $e^{-\mathrm{\Re }{[\mathit{z}]}^2}$ (a quantity already computed), can be used to calculate the range in $\tilde{x}$ space.

With this knowledge, the CDF can be approximated somewhat accurately in $\tilde{x}$ space. Despite the apparent complexity of what was just described, the inversion of the function in the previous paragraph can be done using simple branching logic, and at worst, the solution of a quadratic equation. Because translating the inverse of the previous function into code can take nontrivial effort, C++ code to achieve this has been provided in Appendix E. Figure 6 shows two examples of how this can be a quite satisfactory approximation of the CDF in $\tilde{x}$ space when resonances are influencing the scattering distribution.

Fig. 6. The bootstrapping CDF provides a fairly accurate, easily invertible approximation to the true relative speed CDF to kickstart the root-finding process. The pairs of lines, moving from top to bottom, represent 35.25-eV, 36.25-eV, 38.25-eV, and 66.25-eV incident neutron energies.

IV. RESULTS

IV.A. Calculation of $\mathit{w}(\mathit{z},\mathit{x})$

In order to test the accuracy of Algorithm 3, we have computed reference values of $\mathit{w}(\mathit{z},\mathit{x})$ using the scipy^[39] adaptive quadrature routine, scipy.integrate.quad, to evaluate the integral formulation [Eq. (13)(13) $\mathit{w}(\mathit{z},\mathit{x})=\frac{i}{\pi }{\int }_{-\mathrm{\infty }}^x\frac{e^{-t^2}}{\mathit{z}-\mathit{t}}\mathrm{d}\mathit{t}.$(13) ]. In the approximation of the jump integral [Eq. (47)(47) $\mathit{J}(\mathit{z},\mathit{x})=e^{-\mathrm{\Re }{[\mathit{z}]}^2}{\int }_{\mathrm{\Re }[\mathit{z}]}^z\frac{e^{t^2}\mathrm{d}\mathit{t}}{\mathit{x}-\mathit{t}}.$(47) ], only the first five terms in the series are retained. Where functions such as $log(\mathit{x})$ or $e^x$ appear, C++ standard library implementations have been employed. The implementation of $\mathit{w}(\mathit{z})$ from Ref. [35] has been employed. This results in the error profiles exhibited by Fig. 7, where we have plotted the real part of $\left(w_{\mathrm{approx}}(\mathit{z},\mathit{x})-\mathit{w}(\mathit{z},\mathit{x})\right)/\mathit{w}(\mathit{z})$). Because only the real part is of interest in resonance upscattering calculations, results on the error of the imaginary component are omitted.

Fig. 7. Error of $\mathrm{\Re }[\mathit{w}(\mathit{z},\mathit{x})]$ for a few values of $\mathit{z}$. The legend is the imaginary number added to the real part specified in each figure’s caption. The plotted error $w_{\mathrm{approx}}(\mathit{z},\mathit{x})-\mathit{w}(\mathit{z},\mathit{x})$ is normalized by $\mathrm{\Re }[\mathit{w}(\mathit{z})]$ to match the scaling of Fig. 2. Where the scipy numerical integration routine returned NaNs, the lines disappear.

IV.B. Single Energy Testing

We first present in Fig. 8 the relative speed distribution of ²³⁸U for two different energies and a few temperatures as calculated both by numerical integration and the MARS analytic CDF. The energies correspond to being in the trough and near the peak of a scattering resonance. These plots clearly show the influence of the resonances on the double-differential cross section; a nuclide with a constant cross section has a relative speed distribution that is very nearly an error function at epithermal energies. The relative speed distribution near resonances has a jumping effect that is governed by $\mathit{w}(\mathit{z},\mathit{x})$. They bear a resemblance to the $\mathit{w}(\mathit{z},\mathit{x})$ behavior depicted by Fig. 2.

Fig. 8. The MARS analytic CDF matches the numerically integrated relative speed CDFs. Some errors can be observed for the 1500-K case at 35.25 eV; scattering in the resonance dip is fortunately an extremely rare event.

If the relative speed distribution is correct, the resultant double-differential scattering distribution is also correct. Figure 9 shows this is the case for our method when compared to the RVS method of Ref. [14]. These results were obtained from our modified version of OpenMC.^b It also shows that the pole sampling technique successfully works for ²³⁵U and its tightly spaced resonances.

Fig. 9. Scattering at 1200 K matches results from the RVS method well at two different energies. These energies interact with resonances for both nuclides.

IV.C. Pin Cell Reactivity Feedback

The 2.4%-enriched PWR pin cell example from the OpenMC suite of example problems was used to calculate Doppler reactivity feedback effects with four different models. The first and second used pointwise cross sections that were interpolated between 300 K, 600 K, 900 K, 1200 K, and 2500 K. The model was run at temperatures ranging from 300 to 1800 K in increments of 20 K. Of the two using pointwise cross sections, one used the historical constant cross-section (CXS) free gas scattering approximation and the second used the RVS method. The second two cases both used WMP cross sections, one using RVS and the second MARS. The ENDFB-VII.1 nuclear data set was employed. Figure 10 shows how these cases compare. Two hundred cycles with 10 inactive were employed, using 200 000 particles per cycle. The mean error of k_eff was thus converged to 20 pcm for each case.

Fig. 10. MARS matches the $\mathit{k}$ eigenvalue of the RVS method on a 2.4%-enriched fresh PWR pin cell problem. Line width represents estimated standard deviation of the mean.

It can be seen that the pointwise cross-section representation incurs some interpolation error between 1200 and 1800 K. The MARS method matches the RVS results where multipole cross sections were employed. We can thus conclude that the new method works correctly across the range of energies where resonances influence the double-differential cross section at temperature for nuclides of both strong, distantly spaced resonances (²³⁸U) and closely spaced weak resonances (²³⁵U).

IV.D. Influence on Tracking Rate

Finally, in order to determine the computational efficiency of the new method, tracking rate comparisons were carried out on the same PWR pin cell example problem. The computational performance of both inactive and active cycles was assessed. For the active cycles, a 100 × 100 Cartesian mesh tallied flux, fission rates, and neutron production rates using track-length estimators. In addition, a spatially homogenized energy spectrum tally consisting of 500 equal-lethargy bins was applied.

An Intel Xeon W-2133 with six physical cores carried out the calculations and obtained the results depicted in Table I. The results clearly demonstrate the computational efficiency of MARS compared to the RVS and DBRC methods. While it does not outperform RVS in this scenario, future work will explore its performance on vector computer architectures where we expect it to outperform.

TABLE I Tracking Rate Obtained by the Constant Cross-Section Treatment for Three Resonance Upscattering Models Showing MARS Is Comparable in Speed to Widely Accepted Techniques*

Method	Inactive	Active
CXS	60.5	11.4
DBRC	57.0	11.1
RVS	58.3	11.3
MARS	60.3	11.1

*Track rate in thousands of particles per second.

The computational expense incurred by tallying tends to render the performance impact of our new method particularly negligible. Collision estimators could be used on the mesh tally to improve the tracking rate, but we arbitrarily opted for track-length estimators. Due to subtle hardware-related effects, such as cache utilization or branch prediction, the tracking rates of the three resonance upscattering handling methods have different relative performances when comparing active and inactive cycles. Future work will explore detailed performance results on a variety of architectures.

V. CONCLUSION

Multipole formalism carries a variety of advantages compared to pointwise cross sections. Aside from its potential gains in computational efficiency on modern compute architectures, it enables accurate Doppler broadening without a library-size tradeoff,^[8] offers elegant sensitivity quantification, and narrows the gap between R matrix theory and the cross-section representation.^[21] This work develops yet another advantage to the WMP formalism: closed-form resonance upscattering treatment.

We have demonstrated that the new method matches the results obtained by other resonance upscattering techniques. To achieve this, we derived an expression for the target relative speed distribution and identified a novel special function that universally arises in this application. Novel numerical techniques that balance efficiency and accuracy were derived, implemented, and tested. The overall scheme was shown to achieve the same tracking rate as other resonance upscattering modeling methods.

The new method, MARS, overcomes the storage requirements of the RST method^[7] and avoids rejection sampling as employed by other common approaches. Without the need to access intermediate storage, the accesses to global memory can be reduced on GPU architectures. On top of that, the work discrepancy between threads incurred by rejection sampling on GPUs is similarly overcome. Future work will explore the implementation and optimization of this method on GPUs.

Acknowledgments

Any opinions, findings, conclusions, or recommendations expressed in this paper are those of the authors and do not necessarily reflect the views of the U.S. Department of Energy (DOE) Office of Nuclear Energy.

Disclosure Statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was partially supported by the DOE through the Los Alamos National Laboratory. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of the DOE (contract no. 89233218CNA000001). This paper is also based on work partially supported under an Integrated University Program Graduate Fellowship. This research was also partially supported by the Exascale Computing Project (17-SC-20-SC), a collaborative effort of the DOE Office of Science and the National Nuclear Security Administration.

Notes

^a The term under the square root in Eq. (57)(57) ${\sqrt{E}}_{\mathrm{trough}}=\frac{-b+\sqrt{b^2-\mathit{abc}+a^2\mathit{d}}}{\mathit{a}},$(57) can sometimes be negative in the vicinity of nonphysical poles, which are artifacts of the fitting process, and dealing with imaginary quantities in this case is undesirable. Therefore, our calculation uses a linearization of the nonpole cross section at the incident energy instead in this case.

^b The modified version of OpenMC is available at github.com/gridley/openmc/tree/mars.

APPENDIX A

DERIVATION OF Equation (19)(19) $\mathit{w}(\mathit{z},\mathit{x})=\frac{1}{2}(1+\mathrm{erf}(\mathit{x}))\mathit{w}(\mathit{z})+\frac{i{e}^{-x^2}}{\mathit{\pi }}{\int }_0^{\mathrm{\infty }}\frac{e^{-t^2}{e}^{2\mathit{itz}}\mathrm{d}\mathit{t}}{\mathit{i}(\mathit{x}-\mathit{z})+\mathit{t}}.$(19)

The forthcoming discussion has not been made mathematically rigorous for the sake of brevity and the context of a nuclear engineering journal. We begin by defining the auxiliary complex function F(z),

(A.1) $\mathit{F}(\mathit{z})=e^{z^2}\mathit{w}(\mathit{z},\mathit{x}).$(A.1)

The complex line integral theorem can then be applied when $\mathrm{\Im }[\mathit{z}]>0$:

(A.2) $\mathit{F}(\mathit{z})-\mathit{F}(0)={\int }_0^z\frac{\mathrm{d}F}{\mathrm{d}z}{|}_{z=z^{\prime }}{z}^{\prime }.$(A.2)

Computing $\frac{\mathrm{d}F}{\mathrm{d}z}$ and inserting then reveals

(A.3) $\begin{array}{rl}& e^{z^2}\mathit{w}(\mathit{z},\mathit{x})-\mathit{w}(0,\mathit{x})={\int }_0^z\left(2z^{\prime }{e}^{{z^{\prime }}^2}\mathit{w}(z^{\prime },\mathit{x})\right.\\ & \left.-e^{{z^{\prime }}^2}\frac{i}{\pi }{\int }_{-\mathrm{\infty }}^x\frac{e^{-t^2}\mathrm{d}\mathit{t}}{{(z^{\prime }-\mathit{t})}^2}\right)\mathrm{d}{z}^{\prime },\end{array}$(A.3)

where the linearity of integration has been employed and the interchange of differentiation and integration has also been used. The innermost integrals can now be computed exactly, carrying the $z^{\prime }$ through to the integral defining the incomplete Faddeeva function. This results in

(A.4) $\begin{array}{rl}& e^{z^2}\mathit{w}(\mathit{z},\mathit{x})-\mathit{w}(0,\mathit{x})\\ & =\frac{i}{\pi }\left({\int }_0^z\frac{e^{{z^{\prime }}^2}{e}^{-x^2}}{\mathit{x}-z^{\prime }}\mathrm{d}{z}^{\prime }-\sqrt{\mathit{\pi }}(1+\mathrm{erf}(\mathit{x})){\int }_0^z{e}^{{z^{\prime }}^2}\mathrm{d}{z}^{\prime }\right).\end{array}$(A.4)

Recalling that the Faddeeva function can be defined as

(A.5) $\mathit{w}(\mathit{z})=e^{-z^2}\left(1+\frac{2\mathit{i}}{\sqrt{\mathit{\pi }}}{\int }_0^z{e}^{t^2}\mathrm{d}\mathit{t}\right),$(A.5)

we can identify $\mathit{w}(\mathit{z})$ as the trailing term of Eq. (A.4)(A.4) $\begin{array}{rl}& e^{z^2}\mathit{w}(\mathit{z},\mathit{x})-\mathit{w}(0,\mathit{x})\\ & =\frac{i}{\pi }\left({\int }_0^z\frac{e^{{z^{\prime }}^2}{e}^{-x^2}}{\mathit{x}-z^{\prime }}\mathrm{d}{z}^{\prime }-\sqrt{\mathit{\pi }}(1+\mathrm{erf}(\mathit{x})){\int }_0^z{e}^{{z^{\prime }}^2}\mathrm{d}{z}^{\prime }\right).\end{array}$(A.4) . The term $\mathit{w}(0,\mathit{x})$ must be interpreted in a principal value sense, which results in a contribution in the form of a Heaviside function. The following expression then results:

(A.6) $\begin{array}{rl}& \mathit{w}(\mathit{z},\mathit{x})=e^{-z^2}\left(-\frac{1}{2}(\mathrm{Ei}(-x^2)+\frac{i}{\pi }{e}^{-x^2}{\int }_0^z\frac{e^{{z^{\prime }}^2}}{\mathit{x}-z^{\prime }}\mathrm{d}{z}^{\prime }\right)\\ & +\frac{1}{2}(1+\mathrm{erf}(\mathit{x}))(\mathit{w}(\mathit{z})-e^{-z^2})+\mathit{h}(\mathit{x})e^{-z^2}.\end{array}$(A.6)

This result could perhaps be used for the numerical calculations of $\mathit{w}(\mathit{z},\mathit{x})$. However, it suffers the shortcoming that the exponential integral term goes to infinity for $\mathit{x}=0$, which is canceled out by the integral term. However, $\mathit{w}(\mathit{z},\mathit{x})$ is well defined at $\mathit{x}=0$, and the addition of branching logic to numerical routines to handle this case would be cumbersome. The expression can be made more amenable to numerical approximation with some further simplification.

The integral term goes from 0 to $\mathit{z}$, and the integrand encloses no poles of the following path whenever $\mathit{x}\ne 0$. As such, a change of integration path is employed; the contour of Fig. A.1 results in a convenient cancellation of terms. The top leg of the contour is zero from the $e^{t^2}$ term, resulting in

(A.7) ${\int }_0^z\frac{e^{t^2}}{\mathit{x}-\mathit{t}}\mathrm{d}\mathit{t}={\int }_0^{\mathit{i}\mathrm{\infty }}\frac{e^{t^2}}{\mathit{x}-\mathit{t}}\mathrm{d}\mathit{t}-{\int }_z^{\mathit{i}\mathrm{\infty }}\frac{e^{t^2}}{\mathit{x}-\mathit{t}}\mathrm{d}\mathit{t}.$(A.7)

Fig. A.1. Modified contour used to cancel the exponential integral Eq. (A.6)(A.6) $\begin{array}{rl}& \mathit{w}(\mathit{z},\mathit{x})=e^{-z^2}\left(-\frac{1}{2}(\mathrm{Ei}(-x^2)+\frac{i}{\pi }{e}^{-x^2}{\int }_0^z\frac{e^{{z^{\prime }}^2}}{\mathit{x}-z^{\prime }}\mathrm{d}{z}^{\prime }\right)\\ & +\frac{1}{2}(1+\mathrm{erf}(\mathit{x}))(\mathit{w}(\mathit{z})-e^{-z^2})+\mathit{h}(\mathit{x})e^{-z^2}.\end{array}$(A.6) . The contribution from the top line is zero.

A little bit of algebra shows that

(A.8) ${\int }_0^{\mathit{i}\mathrm{\infty }}\frac{e^{t^2}}{\mathit{x}-\mathit{t}}\mathrm{d}\mathit{t}=\frac{1}{2}{e}^{x^2}\left(-\mathit{i\pi }(\mathrm{erf}(\mathit{x})-\mathrm{sign}(\mathit{x}))+\mathrm{Ei}(-x^2)\right).$(A.8)

The sign function term ends up canceling out the Heaviside term upon substitution of Eq. (A.7)(A.7) ${\int }_0^z\frac{e^{t^2}}{\mathit{x}-\mathit{t}}\mathrm{d}\mathit{t}={\int }_0^{\mathit{i}\mathrm{\infty }}\frac{e^{t^2}}{\mathit{x}-\mathit{t}}\mathrm{d}\mathit{t}-{\int }_z^{\mathit{i}\mathrm{\infty }}\frac{e^{t^2}}{\mathit{x}-\mathit{t}}\mathrm{d}\mathit{t}.$(A.7) back to Eq. (A.6)(A.6) $\begin{array}{rl}& \mathit{w}(\mathit{z},\mathit{x})=e^{-z^2}\left(-\frac{1}{2}(\mathrm{Ei}(-x^2)+\frac{i}{\pi }{e}^{-x^2}{\int }_0^z\frac{e^{{z^{\prime }}^2}}{\mathit{x}-z^{\prime }}\mathrm{d}{z}^{\prime }\right)\\ & +\frac{1}{2}(1+\mathrm{erf}(\mathit{x}))(\mathit{w}(\mathit{z})-e^{-z^2})+\mathit{h}(\mathit{x})e^{-z^2}.\end{array}$(A.6) . The second term in Eq. (A.7)(A.7) ${\int }_0^z\frac{e^{t^2}}{\mathit{x}-\mathit{t}}\mathrm{d}\mathit{t}={\int }_0^{\mathit{i}\mathrm{\infty }}\frac{e^{t^2}}{\mathit{x}-\mathit{t}}\mathrm{d}\mathit{t}-{\int }_z^{\mathit{i}\mathrm{\infty }}\frac{e^{t^2}}{\mathit{x}-\mathit{t}}\mathrm{d}\mathit{t}.$(A.7) easily can be transformed via the change of variables $t^{\prime }=-\mathit{it}$ to the integral in Eq. (19)(19) $\mathit{w}(\mathit{z},\mathit{x})=\frac{1}{2}(1+\mathrm{erf}(\mathit{x}))\mathit{w}(\mathit{z})+\frac{i{e}^{-x^2}}{\mathit{\pi }}{\int }_0^{\mathrm{\infty }}\frac{e^{-t^2}{e}^{2\mathit{itz}}\mathrm{d}\mathit{t}}{\mathit{i}(\mathit{x}-\mathit{z})+\mathit{t}}.$(19) , thus yielding Eq. (19)(19) $\mathit{w}(\mathit{z},\mathit{x})=\frac{1}{2}(1+\mathrm{erf}(\mathit{x}))\mathit{w}(\mathit{z})+\frac{i{e}^{-x^2}}{\mathit{\pi }}{\int }_0^{\mathrm{\infty }}\frac{e^{-t^2}{e}^{2\mathit{itz}}\mathrm{d}\mathit{t}}{\mathit{i}(\mathit{x}-\mathit{z})+\mathit{t}}.$(19) .

APPENDIX B

DERIVATION OF EQUATION (46)(46) $n^{\ast }\approx \mathit{ex}-\frac{1}{2}log\mathit{\pi }.$(46)

The goal is to find the smallest $\mathit{n}$ such that

(B.1) $\frac{x^n{E}_1(\mathit{x})}{\mathit{n}!E_{\mathit{n}+1}(\mathit{x})}<1.$(B.1)

The variation of the left side function of Eq. (B.1)(B.1) $\frac{x^n{E}_1(\mathit{x})}{\mathit{n}!E_{\mathit{n}+1}(\mathit{x})}<1.$(B.1) as $\mathit{n}$ increases, for a constant value of $\mathit{x}$, has been shown to be first increasing from one, reaching a maximum, and monotonically decreasing from that point, never becoming negative.^[31] First, the exponential integrals are replaced with the equivalent upper incomplete gamma function

(B.2) $E_n(\mathit{x})=x^{\mathit{n}-1}\mathrm{\Gamma }(1-\mathit{n},\mathit{x})$(B.2)

so the equation becomes

(B.3) $\mathrm{\Gamma }(0,\mathit{x})=\mathit{n}!\mathrm{\Gamma }(-\mathit{n},\mathit{x}).$(B.3)

Using the asymptotic formula for the upper incomplete gamma function $\mathrm{\Gamma }(\mathit{s},\mathit{x})\to x^{\mathit{s}-1}{e}^{-\mathit{x}}$ gives this approximation to solve

(B.4) $x^n\approx \mathit{n}!,$(B.4)

which can be solved approximately by first inserting Stirling’s formula,

(B.5) $x^n\approx \sqrt{2\mathit{\pi n}}{\left(\frac{n}{e}\right)}^{\mathit{n}}.$(B.5)

Taking the n’th root results in

(B.6) $\mathit{ex}\approx (2\mathit{\pi n}{)}^{\frac{1}{2\mathit{n}}}\mathit{n}.$(B.6)

where the second term on the right can be approximated by expanding the exponential for large $\mathit{n}$:

(B.7) $\mathit{ex}\approx (1+\frac{1}{2\mathit{n}}log(2\mathit{\pi n}))\mathit{n}.$(B.7)

This can be solved exactly in terms of the Lambert W function:

(B.8) $\mathit{n}\approx \frac{1}{2}\mathit{W}\left(\frac{{\mathit{e}}^{2\mathit{ex}}}{\mathit{\pi }}\right).$(B.8)

The asymptotic property of the Lambert W function that $\mathit{W}(\mathit{z})\approx log\mathit{z}-loglog\mathit{z}$ is then used to obtain Eq. (46)(46) $n^{\ast }\approx \mathit{ex}-\frac{1}{2}log\mathit{\pi }.$(46) .

APPENDIX C

ASYMPTOTIC APPROXIMATION FOR $\mathit{R}(\mathit{m},\mathit{z})$ OF EQUATION (30)(30) $\mathit{R}(\mathit{z},\mathit{m})={\int }_0^{\mathrm{\infty }}\frac{e^{-t^2}sin(2\mathit{zt})\mathrm{d}\mathit{t}}{t^2+m^2},$(30)

An efficient approximation can be found by grouping the integrand as

(C.1) $\mathit{R}(\mathit{m},\mathit{z})={\int }_0^{\mathrm{\infty }}{f}^{(0)}(\mathit{t},\mathit{m})sin(2\mathit{zt})\mathrm{d}\mathit{t}=-\frac{1}{2\mathit{z}}{\int }_0^{\mathrm{\infty }}{f}^{(0)}(\mathit{t},\mathit{m})\frac{\mathrm{\partial }}{\mathrm{\partial t}}\left[cos(2\mathit{zt})\right]\mathrm{d}\mathit{t},$(C.1)

where

$f^{(0)}(\mathit{t},\mathit{m})=\frac{e^{-t^2}}{t^2+m^2}.$

Integrating the rightmost expression in Eq. (C.1)(C.1) $\mathit{R}(\mathit{m},\mathit{z})={\int }_0^{\mathrm{\infty }}{f}^{(0)}(\mathit{t},\mathit{m})sin(2\mathit{zt})\mathrm{d}\mathit{t}=-\frac{1}{2\mathit{z}}{\int }_0^{\mathrm{\infty }}{f}^{(0)}(\mathit{t},\mathit{m})\frac{\mathrm{\partial }}{\mathrm{\partial t}}\left[cos(2\mathit{zt})\right]\mathrm{d}\mathit{t},$(C.1) by parts repeatedly results in a divergent series approximation of the form

(C.2) $\mathit{R}(\mathit{m},\mathit{z})=\frac{{\mathit{f}}^{(0)}(0,\mathit{m})}{2\mathit{z}}+\frac{{\mathit{f}}^{(2)}(0,\mathit{m})}{8z^3}+\frac{{\mathit{f}}^{(4)}(0,\mathit{m})}{32z^5}+\frac{{\mathit{f}}^{(6)}(0,\mathit{m})}{128z^7}+\cdots ,$(C.2)

where $f^{(\mathit{n})}(\mathit{t},\mathit{m})$ denotes the n’th derivative of $f^{(0)}(\mathit{m},\mathit{t})$ with respect to $\mathit{t}$. Some of the subsequent values evaluated about $\mathit{t}=0$ are

(C.3) $f^{(2)}(0,\mathit{m})=-\frac{2(1+m^2)}{m^4},$(C.3)

(C.4) $f^{(4)}(0,\mathit{m})=\frac{12(2+2m^2+m^4)}{m^6},$(C.4)

(C.5) $f^{(6)}(0,\mathit{m})=-\frac{120(6+6m^2+3m^4+m^6)}{m^8},$(C.5)

(C.6) $f^{(8)}(0,\mathit{m})=\frac{1680(24+24m^2+12m^4+4m^6+m^8)}{m^{10}}.$(C.6)

While seemingly progressing without a clear pattern, after a considerable amount of staring at these expressions a simple recursive formula can be obtained to compute these values: prime for computer implementation. Consider the sequences $a_n,c_n\in \mathbb{R}$ defined by

(C.7) $a_0=m^{-2};c_0=2$(C.7)

and

(C.8) $a_{\mathit{n}+1}=\frac{2n(2n-1)a_n+c_n}{m^2};c_{\mathit{n}+1}=(4\mathit{n}+2)c_n.$(C.8)

Using this, one can show that $f^{(2\mathit{n})}(0,\mathit{m})=a_n$. This allows for easy evaluation of the asymptotic series of Eq. (C.2)(C.2) $\mathit{R}(\mathit{m},\mathit{z})=\frac{{\mathit{f}}^{(0)}(0,\mathit{m})}{2\mathit{z}}+\frac{{\mathit{f}}^{(2)}(0,\mathit{m})}{8z^3}+\frac{{\mathit{f}}^{(4)}(0,\mathit{m})}{32z^5}+\frac{{\mathit{f}}^{(6)}(0,\mathit{m})}{128z^7}+\cdots ,$(C.2) . Numerical experimentation has shown this to be an excellent approximation with a maximum error around ${10}^{-6}$ when $|\mathit{z}|=5$ and $|\mathit{m}|=1$, retaining only five terms. The error rapidly falls from there as $|\mathit{z}|\to \mathrm{\infty }$ and $|\mathit{m}|\to \mathrm{\infty }$.

APPENDIX D

ASYMPTOTIC APPROXIMATION FOR $\mathit{J}(\mathit{z},\mathit{x})$ OF EQUATION (47)(47) $\mathit{J}(\mathit{z},\mathit{x})=e^{-\mathrm{\Re }{[\mathit{z}]}^2}{\int }_{\mathrm{\Re }[\mathit{z}]}^z\frac{e^{t^2}\mathrm{d}\mathit{t}}{\mathit{x}-\mathit{t}}.$(47)

Compared to the asymptotic approximation for the $\mathit{R}(\mathit{m},\mathit{z})$ integral, a clean expression for simple computer code is not available to our knowledge. Obtaining an asymptotic expression thus relies on access to a computer algebra system. Finding this starts by applying a simple change of variables to Eq. (47)(47) $\mathit{J}(\mathit{z},\mathit{x})=e^{-\mathrm{\Re }{[\mathit{z}]}^2}{\int }_{\mathrm{\Re }[\mathit{z}]}^z\frac{e^{t^2}\mathrm{d}\mathit{t}}{\mathit{x}-\mathit{t}}.$(47) to find

(D.1) $\mathit{J}(\mathit{z},\mathit{x})={\int }_0^{\mathrm{\Im }[\mathit{z}]}\frac{e^{-t^2}{e}^{2\mathit{i}\mathrm{\Re }[\mathit{z}]\mathit{t}}\mathrm{d}\mathit{t}}{\mathit{i}(\mathrm{\Re }[\mathit{z}]-\mathit{x})-\mathit{t}},$(D.1)

where it becomes clear that the numerical difficulty from expanding the exponential term originates from the $e^{2\mathit{i}\mathrm{\Re }[\mathit{z}]\mathit{t}}$ modulation. This is the term to isolate to obtain the correct asymptotic behavior as the integrand becomes increasingly oscillatory. The standard repeated integration by parts procedure can then be applied. This is pure tedium, so we simply report the C++ code that evaluates five terms as follows:

APPENDIX E

ROOT-FINDING BOOTSTRAP FUNCTION

Given the approximate CDF of the relative speed at x = 0, its derivative, the probability associated with the resonance, and the resonance parameters $\mathit{z}$, this function quickly approximates the relative speed quantile function.