Geant4 Version: 11.4.0
Operating System: Windows 11 (64-bit)
Compiler/Version: MSVC 19.50 (Visual Studio 2026)
CMake Version: 4.3.0
Hi,
I am trying to simulate a Ra–Be neutron source in Geant4 with the following setup:
I generate Ra-226 ions at the origin using G4IonTable
I am using G4RadioactiveDecayPhysics to model the decay
The geometry consists of:
A Beryllium spherical shell (0–5 mm)
A detector shell (5–6 mm) to record neutron energy
I record neutron energies when they enter the detector
I also included the command:
/process/had/rdm/thresholdForVeryLongDecayTime 1.0e+60 s
Problem:
Even after running 1,000,000 events , I get only a few neutron entries (around 6–8 values). Most events produce no neutrons.
Example output:
2.44
8.54
0.60
4.54
...
My understanding:
It seems that Ra-226 is not decaying frequently enough during the simulation due to its long half-life (~1600 years), so very few alpha particles are produced, and hence very few (α,n) reactions occur.
Question:
Is this expected behavior when simulating real radioactive decay of long-lived isotopes?
What is the correct way to simulate a Ra–Be source efficiently in Geant4?
Is there any paper/articles which have reported the correct simulation?
You could test if Ra226 is decaying properly by changing everything in your simulation to G4_Galactic (Vacuum) and making your detector detect Alphas. Then the alpha count should be equal to the primary count, since the daughter isotope will have a momentum and leave the simulation.
With 1e6 simulated events, you would expect around 400 neutrons, considering a neutron emmision ratio of ~1.5e4 neutrons/(s*mCi) (- this number comes from an old Amersham catalogue i am having around).
Ra-226 has a decay chain of 5 alpha decays until it is stable, so is this decay chain is not processing, then your neutron count is of by a factor of 5, which is still not enough to be close to the ‘theoretical’ value i would expect…
What physics list and cross section table are you using for the production of neutrons, i.e. the Be(a,xn)Y-Reaction?
Thanks for your detailed response and suggestions.
I performed the test you suggested by setting the entire geometry to vacuum (G4_Galactic) and modifying the stepping action to detect alpha particles instead of neutrons. I also printed the first few events to verify decay.
However, I did not observe any alpha particles even for the first few events. This suggests that Ra-226 is not decaying within the simulation, most likely due to its long half-life (~1600 years) and the lack of sufficient time evolution in the simulation.
So it seems that the issue is not with the (α,n) interaction, but rather that the radioactive decay itself is not occurring.
Regarding your earlier questions:
I am using the following physics list:
G4HadronElasticPhysicsHP
G4HadronPhysicsQGSP_BIC_HP
G4IonPhysics
G4IonElasticPhysics
G4StoppingPhysics
G4EmStandardPhysics_option3
G4DecayPhysics
G4RadioactiveDecayPhysics
For neutron production, I am relying on the default Geant4 hadronic models. I have not included any external cross-section libraries (e.g., TENDL or JENDL) specifically for the Be(α,n) reaction.
From this, I suspect that:
The radioactive decay is not being triggered due to the long half-life.
Even if decay is enabled, the (α,n) reaction modeling for Be may be limited in default Geant4.
I would appreciate your guidance on the following:
What is the correct way to handle long-lived isotopes like Ra-226 in Geant4?
Should I use biasing or time scaling to force decay?
For simulating a Ra–Be neutron source, is it standard practice to:
Replace the Ra decay chain with an equivalent alpha source?
Or include evaluated cross-section libraries (like TENDL/JENDL) for Be(α,n)?
Are there any recommended examples or references for correctly modeling (α,n) neutron sources in Geant4?
Thanks again for your help, and I will also go through the AmBe thread you shared.
Do you set the energy of the Ra226 to be 0 MeV?
If you dont set any energy, Geant4 will create primaries with 1 MeV by default, i think.
Then your Ra226 will simply leave your simulation (world-)volume before it decays.
What Max has said about the particles needing to decay first is correct, and they correctly linked to the AmBe simulation.
You have indeed hit the problem of any neutron generator source in Geant4 that use the 9Be(a,n) reaction: being strongly anisotropic, it strictly depends on differential cross section. Geant4 does not have these.
I wrote a simulation for an AmBe source: Simulation and analysis of an AmBe source - ScienceDirect (open access) and made it open source. On the relative GitHub repository GitHub - bhamnuclear/AmBeSimulation · GitHub you can find examples for PuBe and CmBe. There, you can create your input isotope card and your Eloss equations. I would suggest using Lise++ for doing this, as it is one of the few codes that can do energy loss above natural uranium masses.
You will need to edit the PrimaryGenerator to load the created isotope card (see lines 165-168), and also PrimaryGeneratorAction::Eloss to implement your energy loss equations.
You can follow the readme in the github repo.
If you have any questions, please let me know.
Max - you were absolutely right about the primary energy. I had not set it initially, and the Ra-226 ions were leaving the world volume before decaying. After setting the energy to 0 MeV, the decay chain started working correctly.
I also investigated the decay chain behavior further. By using:
/process/had/rdm/analogueMC false
I am now able to observe the full decay cascade, and I can see multiple alpha particles per event (around 15–20), corresponding to different branches and excited states.
Filippo - thanks a lot for pointing out the limitations of Geant4 for the 9Be(a,n) reaction. I understand now that the lack of differential cross-section data means the neutron energy and angular distributions will not be accurate if I rely purely on Geant4 physics.
For now, I am focusing on correctly simulating the Ra-226 decay and alpha spectrum. As a next step, I plan to:
Implement a basic a + Be setup in Geant4 to understand the interaction.
Then move towards a more realistic approach using external data (as you suggested), possibly using your open-source repository and tools like LISE++.
I will go through your GitHub repository and paper in detail. If I have questions while implementing the isotope card or energy loss model, I will reach out.
Thanks again for your guidance - it clarified several important aspects of the simulation.
1- The 5 alpha created in the decay chain of Ra226 are in [`~4 MeV, ~8 MeV]. See and run Ra226.mac.txt (349 Bytes)
2- such alpha have a nuclear interaction length (aka MeanFreePath) in Beryllium of the order of 10 - 20 cm. See and run hadronic03.mac.txt (256 Bytes)
( in PhysicsList.cc, use G4G4IonPhysicsXS or G4IonPhysicsPHP)
3- however, those alpha have a stopping range of the order of 20 - 50 um (eg. energy lost by ionisation). See and run Em0.mac.txt (229 Bytes)
It means that alpha, most of the time, loose their kinetic energy in few micrometers and stop instead of interact with Beryllium nuclei
4- to illustrate this, hadronic06.mac.txt (704 Bytes) is a simulation similar to yours: alpha source in a sphere of Beryllium.
In hadronic06.out.txt (1.8 KB) , I started directly with alpha 5 MeV : 100 000 alpha.
As you can see, there are only 12 alphaInelastic interactions, among them 6 (C12 + n) produced. The rest of alpha are stopped by ionisation. ( in PhysicsList.cc, I used G4IonPhysicsXS)
ps. You can start from Ra226. Results are the same, but more complicated to read.