High energy limit should be linear combination of be the free atom cross section of isotopes

[marquezj]

Currently materials in NCrystal are represented as a system of atoms. At high energies, the cross section goes to the free atom cross section that can be computed from the bound atom cross section of the element, and the average mass of the element. The problem is the high energy limit should actually be the abundance weighted average of the free atom cross sections for each isotope. This is caused by the interaction of high energy neutrons with independent nuclei, which can only be independent nuclides, not elements.

The approximation of using elemental cross section is mostly fine elements for which one isotope is dominating, and also for heavier elements, because the difference in mass between the atoms is small. But the differences for light elements are dramatic.

An example of this is the calculation of the free atom cross section of mixtures of light and heavy water:

import numpy as np
import NCrystal as NC
import matplotlib.pyplot as plt

# Digitized from Fig.2, https://doi.org/10.1103/PhysRevLett.90.105302
txt = """0.000 1.493
0.150 1.322
0.300 1.148
0.400 1.033
0.500 0.920
0.900 0.466
1.000 0.350
"""
exp_blostein = np.fromstring(txt,sep=' ')
exp_blostein.shape = (-1,2)

fracD = np.linspace(0.001,0.999,100)
free_XS_O = NC.atomDB(Z=8, A=16).freeScatteringXS()
free_XS_H = NC.atomDB(Z=1, A=1).freeScatteringXS()
free_XS_D = NC.atomDB(Z=1, A=2).freeScatteringXS()
xs = []
xs_linear = []
rho = 3.0*0.6022/18.015 # at/b*cm or at/Aa^3
#
# Compute the cross section at 10 eV for different H2O/D2O mixtures
#
for x in fracD:
    scat = NC.createScatter(f'freegas::H2O/{rho}perAa3/H_is_{1-x}_H1_{x}_H2/O_is_O16')
    xs.append(rho*scat.crossSectionIsotropic(10.0))
    xs_linear.append(rho*1.0/3.0*(free_XS_O+2*((1-x)*free_XS_H+x*free_XS_D)))
plt.plot(fracD, xs, label='elemental FG')
plt.plot(fracD, xs_linear, label='isotopic combination')
plt.plot(exp_blostein[:,0], exp_blostein[:,1], 'o', label='Blostein (2003)')
plt.legend()
plt.xlabel('D fraction')
plt.ylabel('Macro XS [cm^-1]')

As expected, the epithermal total cross section is the linear combination of the free atom cross sections for the independent isotopes:

A possible solution for this is that the inelastic cross sections would be computed independently for each isotope of the material, and their contribution would be added using the natural abundance or an user provided isotopic distribution.

[tkittel]

Yes, there is clearly an issue here of how to transfer into the high energy region. Of course, one option (given that the abundancy is somehow available), is to not use a single-mass free-gas model for the highE extrapolation, but rather add one free-gas SAB contribution for each isotope. Still, for consistency one would actually need to have the lowE parts separated out as well.

However, I am surprised that you can get NCrystal to not give you the correct limit. Try something higher than 10eV.

[tkittel]

Are the blostein values actually measured at 10eV?

[marquezj]

I chose 10 eV as an energy where the value of the XS already is close to the free atom value, but that is not relevant. You can calculate the same using the value from freeScatteringXS().

import numpy as np
import NCrystal as NC
import matplotlib.pyplot as plt

# Digitized from Fig.2, https://doi.org/10.1103/PhysRevLett.90.105302
txt = """0.000 1.493
0.150 1.322
0.300 1.148
0.400 1.033
0.500 0.920
0.900 0.466
1.000 0.350
"""
exp_blostein = np.fromstring(txt,sep=' ')
exp_blostein.shape = (-1,2)

fracD = np.linspace(0.001,0.999,100)
free_XS_O = NC.atomDB(Z=8, A=16).freeScatteringXS()
free_XS_H = NC.atomDB(Z=1, A=1).freeScatteringXS()
free_XS_D = NC.atomDB(Z=1, A=2).freeScatteringXS()
xs = []
xs_linear = []
rho = 3.0*0.6022/18.015 # at/b*cm or at/Aa^3
#
# Compute the cross section at 10 eV for different H2O/D2O mixtures
#
for x in fracD:
    xsf = 0.0
    info_HD = NC.createInfo(f'freegas::H2O/{rho}perAa3/H_is_{1-x}_H1_{x}_H2/O_is_O16')
    for di in info_HD.dyninfos:
        xsf = xsf + di.fraction*di.atomData.freeScatteringXS()
    xs.append(rho*xsf)
    xs_linear.append(rho*1.0/3.0*(free_XS_O+2*((1-x)*free_XS_H+x*free_XS_D)))
plt.plot(fracD, xs, label='elemental FG')
plt.plot(fracD, xs_linear, label='isotopic combination')
plt.plot(exp_blostein[:,0], exp_blostein[:,1], 'o', label='Blostein (2003)')
plt.legend()
plt.xlabel('D fraction')
plt.ylabel('Macro XS [cm^-1]')

Also: it is not important either that I used a free gas model. The same thing happens with a solid: the problem comes from considering a single scatterer for the element with the average mass:

c_HD = NC.NCMATComposer()
c_HD.set_dyninfo_msd('H',msd=0.022,temperature=50,fraction=1.0) # msd is not real, just an example
c_HD.set_density(0.1,'g/cm3') # density is not real, just an example
c_HD.set_composition(f'H is 0.5 H1 0.5 H2')
c_HD.register_as('hd.ncmat')
info_HD = NC.createInfo(f'hd.ncmat')
print(info_HD.dyninfos[0].atomData)

# H=H{50%H1+50%H2}(cohSL=1.4652fm cohXS=0.269776barn incXS=44.5655barn absXS=0.166559barn mass=1.51096amu Z=1)

[marquezj]

Yes, there is clearly an issue here of how to transfer into the high energy region. Of course, one option (given that the abundancy is somehow available), is to not use a single-mass free-gas model for the highE extrapolation, but rather add one free-gas SAB contribution for each isotope. Still, for consistency one would actually need to have the lowE parts separated out as well.

I agree.

[ramic-k]

Hi Thomas,

For SCALE integration of NCrystal this issue would be beneficial as well to connect TSL cross sections to the higher energy cross sections.