Incorrect photon mean free path when Rayleigh and photonuclear used together

[ftessier]

Description

The photon mean free path is incorrect when both Rayleigh scattering = On and Photonuclear attenuation = On. Given a density $\rho$, and a cross section $\sigma_0 \equiv (\sigma_\text{pair} + \sigma_\text{Compton} + \sigma_\text{photoelectric})$, the current implementation of Rayleigh and Photonuclear corrections implies that the photon mean free path (gmfp in the code) is set to

$$ \lambda' = \frac{1}{\rho\sigma_0} \left(\frac{\sigma_0}{\sigma_0 + \sigma_\text{Rayleigh}}\right) \left(\frac{\sigma_0}{\sigma_0 + \sigma_\text{photonuclear}}\right). $$

Expected behaviour

Strictly speaking, the mean free path of the photon should be set to

$$ \lambda = \frac{1}{\rho} \left(\frac{1}{\sigma_0 + \sigma_\text{Rayleigh} + \sigma_\text{photonuclear}}\right). $$

Impact

• This bug, in principle, affects any simulation where both effects are turned on, but in practice it is unlikely to induce a discernable bias in typical simulations, for example calculations of deposited dose.

• It has low impact because Rayleigh scattering is significant for photon energies in the 0.01-0.1 MeV range, while photonuclear attenuation is relevant above about 10 MeV. They don't overlap. Hence, at least one bracket in $\lambda'$ is very close to unity for any given photon energy.

• It was introduced when photonuclear attenuation was added to EGSnrc, around 2012.

Root cause

• The Rayleigh interaction was originally implemented as a correction factor to the mean free path $\lambda_0 \equiv 1/(\rho\sigma_0)$, correctly, as in

$$ \lambda_\text{R} = \lambda_0\left(\frac{\sigma_0}{\sigma_0 + \sigma_\text{Rayleigh}}\right) = \frac{1}{\rho}\left(\frac{1}{\sigma_0 + \sigma_\text{Rayleigh}}\right). $$

• The photonuclear attenuation was then modelled on the Rayleigh one, using a similar correction factor, to yield:

$$ \lambda' = \lambda_\text{R}\left(\frac{\sigma_0}{\sigma_0 + \sigma_\text{photonuclear}}\right). $$

Proposed solution

• It would be possible to adjust the photonuclear correction factor, e.g., $(\sigma_0 + \sigma_\text{Rayleigh})/(\sigma_0 + \sigma_\text{Rayleigh} + \sigma_\text{photonuclear})$. But then a new variable is needed for branching, interactions become interrelated, and the logic is convoluted further.

• A more consistent solution is to remove the Rayleigh and photonuclear correction factors altogether, and instead add branching ratios for these interactions. All photon interactions are then handled in a the same way.

• A disadvantage of refactoring the Rayleigh and photonuclear logic is that the meaning of the cohfac and photonucfac will change: they will become cumulative branching ratios (similar to gbr1 and gbr2), instead of gmfp correction factors. Hence all code depending on these factors, and more generally on the cohe or photonuc data arrays, has to be carefully reviewed.

• Notably, cross section enhancement is implemented in EGSnrc applications by redefining the $RAYLEIGH-CORRECTION macro. This macro should become a no-op, and a new macro $CROSS-SECTION-ENHANCEMENT ought to be added for this purpose. Again, any code calling the former macro has to be reviewed and adjusted.