I was investigating the current status of the Argon neutron cross sections. On the way I also stumbled upon what I think is a bug in the skip_missing_isotopes flag.
For context, I work on the LEGEND experiment, which searches for neutrinoless double-beta decay. We run high purity germanium detectors, enriched in Ge-76, bare in a liquid argon cryostat. One of our backgrounds is the production of unstable isotopes by neutrons from muon-induced shower. Argon has a relatively high muon-induced neutron yield, so we care about how well their transport and capture are modeled.
Our previous studies used G4.10.5 and I now want to move to a newer version. In the process I stumbled upon the argon changes in the [release notes] of 11.2.2, also discussed in this [thread]. I understand that the Ar-36 and Ar-38 cross-sections from JEFF-3.3 are imprecise, and that this is why G4NDL 4.7.1 removed them. I expected argon to then behave close to pure Ar-40 as the two isotopes are only 0.4 % of natural argon. A thermal neutron should capture on Ar-40 and give Ar-41, and the two removed isotopes should no longer contribute (using skip_missing_isotopes).
That is not what seems to happen when I ran a few tests based on Hadr04. I started 2000 thermal neutrons in a 3 m cube of G4_lAr and counted the nuclides produced. I repeated the test with three libraries: G4NDL 4.7 (before the removal), G4NDL 4.7.1 (after it) and ENDF/B-VIII.0 (as a reference). Each library ran with /process/had/particle_hp/skip_missing_isotopes set to true and to false. The results are in the table below.
| library | skip_missing_isotopes |
substitutions | Ar-41 | Ar-39 | Ar-37 | S-33 | Cl-36 |
|---|---|---|---|---|---|---|---|
| G4NDL 4.7.1 | true | 6 | 163 | 78 | 3 | 1691 | 65 |
| G4NDL 4.7.1 | false | 6 | 163 | 78 | 3 | 1691 | 65 |
| G4NDL 4.7 | true | 0 | 1921 | 2 | 57 | 0 | 0 |
| G4NDL 4.7 | false | 0 | 1921 | 2 | 57 | 0 | 0 |
| ENDF/B-VIII.0 | true | 0 | 1902 | 2 | 50 | 0 | 0 |
| ENDF/B-VIII.0 | false | 0 | 1901 | 2 | 50 | 0 | 0 |
G4NDL 4.7 and ENDF/B-VIII.0 give what I would expect. Almost every neutron captures on Ar-40 and gives Ar-41. Ar-36 adds a small Ar-37 contribution, and Ar-38 a negligible Ar-39 one.
G4NDL 4.7.1 gives only 163 Ar-41 out of 2000 neutrons. At the same time the Ar-39 production is about forty times too high, and the Ar-37 production about twenty times too low. This is what you would get if Ar-38 takes the cross section of Ar-39, and Ar-36 that of Ar-37 as per substitution which the run also warns about regarding these isotopes. Also notice the production of S-33 and Cl-36 via (n, alpha) and (n, p) reactions on Ar-36 even with kinetic energies of 25 meV.
What I did not expect is that skip_missing_isotopes=true changes nothing. I assumed it suppresses the substitution and sets those cross sections to zero. The Ar-39 and Ar-37 production should then drop. Instead the numbers are identical to the false case.
In the G4 source code the NDL libraries seem to be read at two separate points. During transport the cross section comes from G4CrossSectionHP, which builds the data per element in Initialise(Z). There it applies the condition18 == Z && 40 != A -> continue [link], so Ar-36 and Ar-38 are never considered.
Once an interaction is decided, the model picks the target isotope and the final state through G4ParticleHPChannel. This class builds its own data with G4ParticleHPIsoData::Init and G4ParticleHPNames::GetName, and that path has no such condition. GetName falls back to the nearest available mass number, A+1 first. Ar-36 therefore gets the Ar-37 files and Ar-38 the Ar-39 files [link].
The skip flag therefore does not work as I understand it. This is the messenger output
=== G4ParticleHPMessenger CHANGED PARAMETER SkipMissingIsotopes TO 1 ===
-> Use only exact isotope data files, instead of allowing nearby isotope files to be used:
if the exact file is not available, the cross section will be set to zero !
and this is the description in the Application Developer Guide
/process/had/particle_hp/skip_missing_isotopes true
This UI commands sets to zero the cross section of the isotopes which are not present in the neutron library. If Geant4 doesn't find an isotope, then it looks for the natural composition data of that element. Only if the element is not found then the cross section is set to zero. On the contrary, if this variable is not defined, Geant4 looks then for the neutron data of another isotope close in Z and A, which will have completely different nuclear properties and lead to incorrect results (highly recommended).
G4NDL holds no 18_nat_Argon file. By that description argon should therefore get a zero cross section, not a substitution.
In GetName the branch that sets the cross section to zero has the condition fManager->GetSkipMissingIsotopes() && (Z != result.GetZ() || !result.IsThisNaturalAbundance()) [link]. For a substitution inside the same element the first part of the bracket is false. The second part is false as well. nat in G4ParticleHPDataUsed is initialised to true, and the only function that can change it, SetNaturalAbundanceFlag(), also sets it to true. IsThisNaturalAbundance() therefore always returns true. The condition seems to me only be true when the fallback crosses to a different element, which for argon it seems to never do.
So my first question is whether this is the expected behavior of skip_missing_isotopes. Is it meant to suppress a substitution within the same element as well? If so, should I report it in bugzilla?
The second question is about the data libraries. I understand from [this thread] that the removal of the cross sections does not imply a recommendation on which library to use. Has any other user tested this and potentially has some advice? (I initially thought about adding this to the above thread, but with the bug I thought it better to be a separate post.)
For my part, I ran the three libraries above on the LEGEND setup and compared the isotope production. Using MUSUN, I started muons around our setup and simulated them using Shielding. All simulations share the same geometry and exposure. The Ar-41 production with G4NDL 4.7.1 is about a factor of eight below the other two, while other isotopes produced via inelastic scattering (e.g. Cl-36, S-33) are significantly boosted. In comparison, G4NDL 4.7 and ENDF/B-VIII.0 give very similar results. If there is interest to repeat it with other datasets, I can deliver them in a matter of a day approximately.
Finally, I am aware of [Capture cross-section issues using Shielding in G4.11.4] and the corresponding [bugzilla 2759]. I ran the tests above with G4.11.3.2 to stay clear of that one. I get identical numbers with G4.11.4, so I do not think the two are the same problem.
I attached the setup files to reproduce the numbers in the table. They are three macros and a README.
README.txt (1.6 KB)
mac_lar_skip_true.txt (271 Bytes)
mac_lar_skip_false.txt (272 Bytes)
_Geant4 Version: 11.3.2 and 11.4.0
_Operating System: Ubuntu 24.04.3 LTS
_Compiler/Version: 13.3.0
_CMake Version: 3.28.3