Techniques for neutron ray-tracing

The neutron ray-tracing packages have a very similar philosophy in the way simulations are described and performed, and in the way the neutrons are represented. We here describe this in some detail.

Representing the neutrons in simulations

In simulations, a neutron is represented semiclasically by simultaneously well-defined position, \({\bf r}\), velocity, \({\bf v}\), and and all three components of the neutron spin vector, \({\bf s}\). Formally, this violates the laws of quantum mechanics, in particular the Heisenberg uncertainty relations, given for the position/momentum coordinates by

\begin{equation}\label{dummy2120409709} \delta x \delta p \geq \dfrac{\hbar}{2} \end{equation}

However, as we shall illustrate in the problem Validity of the semiclassical approximation, the semiclassical approximation is very good for describing instruments that use "typical" neutrons with velocities of the order 0.1 - 10 km/s. One important region where the semiclassical approximations breaks down is for very slow (or "ultra-cold") neutrons, where quantum effects become prominent in the optical properties. Ultra-cold neutrons will not be discussed in this version of the notes.

The validity of the semiclassical approach is discussed in detail in Ref. .

The neutrons are simulated in "rays", by which we mean the neutron trajectory; position \({\bf r}\), as a function of time \(t\). Therefore, at any point in the simulation of one ray, \(t\) is a necessary parameter. This is of particular importance when simulating pulsed neutron sources.

The neutron weight factor

To represent realistic values for neutron numbers, a neutron ray in general represents more than a single physical neutron. As a consequence, the ray contains an additional parameter, the weight factor, \(p\). This generally has the unit of neutrons per second. When the ray begins at the source, \(p\) has a typical initial value of thousands or millions neutrons per second.

To improve simulation speed, the weight factor can be manipulated through the simulation. For example, when some physical neutrons are "lost" due to e.g. finite reflectivity or absorption, the simulated ray will in general continue in the simulations, while \(p\) is adjusted to reflect the correct average physical behaviour. When (if) the neutron ray reaches the detector, \(p\) may be only a fraction of a neutron per second. This weight factor adjustment is a very efficient implementation of importance sampling, as presented in the section Introduction to the Monte Carlo technique.

To quantify the method, we consider a certain position (a surface perpendicular to the beam) in the simulated instrument. Here, the neutron intensity is given by the sum of all rays reaching this surface (again with units neutrons per second):

\begin{equation}\label{dummy945804845} I_j = \displaystyle\sum_{i=1}^N p_{i,j-1} \end{equation}

where \(i\) is the ray index, \(N\) is the total number of rays, and \(j\) is the index of the given component. The index \(j-1\) on \(p\) indicate that we consider the intensity being emitted by the previous component. If the ray does not reach this point, we take \(p_{i,j-1}=0\). The weight of the neutrons emitted from this point in the simulations, i.e. after interacting with component \(j\), is expressed by

\begin{equation}\label{dummy653998044} p_j = w_j p_{j-1} ,\, \end{equation}

where the ray index, \(i\), is omitted for simplicity. The weight multiplier of the \(j\)'th component, \(w_j\), is calculated by the probability rule

\begin{equation}\label{eq:weight_master} f_{\rm MC, b} w_j = P_{\rm b} , \, \end{equation}

where \(P_{\rm b}\) is the physical probability for the event "b", and \(f_{\rm MC,b}\) is the probability that the Monte Carlo simulation selects this event.

Often, there is only one allowed event, giving \(f_{\rm MC}=1\), whence \(w_j=P\). This may, e.g., be the case for neutrons being attenuated when passing through absorbing materials. When a Monte Carlo branch point is reached (selection between several events), we have \(f_{\rm MC, b} < 1\) for each branch, b. However, since \(f_{\rm MC}\) is a probability function, we must have

\begin{equation}\label{dummy540191185} \displaystyle\sum_b f_{\rm MC, b} = 1 . \, \end{equation}

Estimates of simulation uncertainty

In a stochastical simulation, it is important to be able to estimate the uncertainty, in the same way as for experiments. We here present a simple derivation of the uncertainty in simulations with weight factors.

As a simple start, let us imagine \(M\) rays, all with weight factor \(p\). Each ray is imagined to have an overall probability, \(P_{\rm d}\), of reaching the detector. The distribution of observed rays will be binominally distributed with a mean value \(\bar{N} = MP_{\rm d}\) and variance \(\sigma^2(N) = M P_{\rm d} (1-P_{\rm d})\). Very often, we will have small values of \(P_{\rm d}\) and large values of \(M\), so that \(M P_{\rm d} \gg 1\). Here, it is valid to approximate the distribution as a Gaussian with standard deviation

\begin{equation}\label{dummy730577740} \sigma(N) = \sqrt{MP_{\rm d}} = \sqrt{N} , \end{equation}

where the observed value, \(N\), is used as the best estimate for the average, \(\bar{N}\). Recalling that the rays have weight factors, the total simulated intensity will in this case be

\begin{equation}\label{dummy185225129} I=Np, \, \end{equation}

with a variance

\begin{equation}\label{eq:MCerror1} \sigma^2(I) = Np^2 . \, \end{equation}

Imagine now that the rays can have a number of (discrete) different weight factors, \(p_i\). The simulated number of rays in each class is denoted \(n_i\) (standard deviation \(\sqrt{n_i}\)). The total simulated result is now

\begin{equation}\label{dummy2068296395} I = \displaystyle\sum_i n_i p_i . \, \end{equation}

Using that the variance of a sum is the sum of variances, we reach the statistical error of the simulated result:

\begin{equation} \label{eq:MCerror2} \sigma^2(I) = \displaystyle\sum_i n_i p_i^2 . \end{equation}

Now, let the simulations occur with rays of arbitrary weight factors, \(p_j\), different for each ray, i.e. \(n_i \equiv 1\), corresponding to a continuous distribution of \(p_i\). If the distribution of these weight factors is reasonably well behaved, we can generalize the equations above to reach

\begin{equation}\label{dummy2014930404} I = \displaystyle\sum_i p_i , \, \end{equation}

with a statistical error:

\begin{equation} \label{eq:MCerror_final} \sigma^2(I) = \displaystyle\sum_i p_i^2 , \end{equation} which is consistent with \eqref{eq:MCerror1} and \eqref{eq:MCerror2}. This is the way uncertainties of simulation results are calculated in all detectors (sometimes denoted monitors) in the McStas package.

Scattering from a sample

The place where most Monte Carlo choices are made is when the neutron ray interacts with a sample. First, it must be decided whether the ray scatters or not. If scattering takes place, the algorithm must select the point in the sample for the scattering event, the scattered direction, and (for inelastic scattering) the final energy. Finally, there is potentially an issue of multiple (repeated) scattering.

To simplify the description, let us just study the scattering direction. Assume we have a sample that only scatters elastically and isotropically with volume specific cross section \(\Sigma\) and has no absorption. Then, the attenuation factor is \(\mu = \Sigma\), and the physical probability for scattering is

\begin{equation}\label{dummy547165093} P_{\rm scatt} = 1-\exp(-\Sigma L) , \, \end{equation}

where \(L\) is the length of the particular ray within the sample. However, we must also consider the outgoing ray. The probability for scattering into the solid angle \(d\Omega\) is

\begin{equation}\label{dummy1409330943} P(\Omega) d\Omega = P_{\rm scatt.} \dfrac{d\Omega}{4\pi} . \end{equation}

Let us require that the ray must always scatter, and let us select the outgoing ray direction with uniform probability (the "physical" scattering). Then we have

\begin{equation}\label{dummy469886538} f_{\rm MC\, phys.} d\Omega = \dfrac{d\Omega}{4\pi} , \end{equation}

and equation \eqref{eq:weight_master} gives

\begin{equation} \label{eq:w_simple} w_{\rm phys.} = \dfrac{P(\Omega)d\Omega}{f_{\rm MC\, phys.}d\Omega} = 1-\exp(- \Sigma L) \approx \Sigma L , \end{equation}

where the rightmost approximation is valid only for "thin" samples, \(\Sigma L \ll 1\). Note that we have not taken into account the attenuation of the outgoing neutron ray.

Focusing in sample scattering

Some neutron instruments have a geometry such that only rays scattered in certain "interesting directions" have any chance of being detected. In such cases, one will employ the technique of focusing to improve simulation efficiency. This implies a need to modify \eqref{eq:w_simple}, as we discuss below.

In focusing, the ray will be emitted only within a certain solid angle, \(\Delta \Omega\). This solid angle must contain all the directions contributing to the detector counts, otherwise the focusing will give wrong results. Assuming uniform selection within \(\Delta \Omega\), the Monte Carlo probability will be given by

\begin{equation}\label{dummy1611418157} f_{\rm MC\, focus} d\Omega = \dfrac{d\Omega}{\Delta \Omega} , \end{equation}

and equation \eqref{eq:weight_master} gives

\begin{equation}\label{dummy686215173} w_{\rm focus} = \dfrac{P_{\rm scatt.}}{f_{\rm MC\, focus}} = \dfrac{\Delta\Omega}{4\pi} \big(1-\exp(- \Sigma L) \big) \approx \dfrac{\Delta\Omega}{4\pi}\Sigma L . \end{equation}

Comparing to the physical case, the focusing method gives smaller weight factors per ray, but a larger number of rays traveling towards the detector. This gives the same final result, but with a smaller statistical error. The focusing technique is an important example of variance reduction through importance sampling.