A software toolkit to simulate activation background for high energy detectors onboard satellites
G. Galgoczi, J. Ripa, G. Dilillo, M. Ohno, R. Campana, N. Werner
AA software toolkit to simulate activation background for highenergy detectors onboard satellites
G´abor Galg´oczi a,b , Jakub ˇR´ıpa a,c,d , Giuseppe Dilillo e,f , Masanori Ohno a , Riccardo Campana g,h ,Norbert Werner c,a,d , and the HERMES-SP collaboration a Institute of Physics, E¨otv¨os Lor´and University, Budapest, Hungary b Wigner Research Centre, Budapest, Hungary c Department of Theoretical Physics and Astrophysics, Faculty of Science, Masaryk University,Brno, Czech Rep. d Astronomical Institute of Charles University, Prague, Czech Rep. e University of Udine, Udine, Italy f INAF/Astronomical Observatory of Trieste, Trieste, Italy g INAF/Astrophysical and Space Science Observatory (OAS), Bologna, Italy h INFN-Sezione di Bologna, Bologna, Italy
ABSTRACT
A software toolkit for the simulation of activation background for high energy detectors onboard satellites is pre-sented on behalf of the HERMES-SP collaboration. The framework employs direct Monte Carlo and analyticalcalculations allowing computations two orders of magnitude faster and more precise than a direct Monte Carlosimulation. The framework was developed in a way that the model of the satellite can be replaced easily. There-fore the framework can be used for different satellite missions. As an example, the proton induced activationbackground of the HERMES CubeSat is quantified. Keywords: simulation, CubeSat, activation, background
1. INTRODUCTION
The High Energy Rapid Modular Ensemble of Satellites (HERMES) mission consists of a fleet of 3U CubeSatshosting broadband X and gamma-ray detectors. The aim of the mission is to detect transient sources with highsky coverage and good angular resolution. The CubeSats will utilize triangulation to locate the investigatedX-ray sources. Therefore proton induced activation of these satellites have to be taken into account since itincreases background in the long term significantly. An important factor to evaluate for the HERMES missionis the activation background produced by the satellite passages inside the regions with high fluxes of trappedprotons, such as the South Atlantic Anomaly (SAA). Different orbits will have different characteristics (depthand duration) and therefore different irradiation doses. Moreover, decay times of all the radio-activated isotopesshall be evaluated. We have developed a framework similar to Odaka et al. which is much more efficient from acomputational point of view with respect to a full-scale Monte Carlo simulation. This methods consists of threesteps.1. A simulation to determine the production rate of each long lived unstable nuclear isotope.2. Analytical calculation of the decay rates and chains of each isotope at any given time.3. A second simulation involving only the decaying isotopes and their products. Further author information: (Send correspondence to G. Galg´oczi)G. Galg´oczi: E-mail: [email protected] , Telephone: +36 30 497 84 25 a r X i v : . [ a s t r o - ph . I M ] J a n . DETERMINING THE NUMBER OF PRIMARY NUCLEAR ISOTOPES The framework which is presented in this paper is based on Geant4. In the first step a mass model is usedto simulate an irradiation of the satellite with monochromatic proton fluxes. One can think of this step asthe determination of response to the proton irradiation similarly to the Green function technique. 16 energiesbetween 4 MeV and 700 MeV were simulated. Energies lower than 4 MeV does not create activation in ourcase and energies above 700 MeV have typically very low fluxes. In Figure 1 the number of different identifiedisotopes is shown, for a simulation with the HERMES mass model and in Figure 2 the simulation is shown.Higher proton energy results in more possible reactions.
Figure 1. Number of identified species of isotopes created by proton irradiation on the HERMES satellite for eachirradiation energy. Different volumes implemented in the mass model are treated independently, but for this plot theyare summed together.
In the “SteppingAction” Geant4 class radioactive decays are detected. Then to prevent recording a decaychain several times, the decayed isotope and the simulation of its secondary particles is stopped. For each volumeof the geometry the name, energy, and excitation level of these isotopes are stored with their excitation energy.Each of the “Physical” volumes in the mass model have to be treated independently. The reason for this isdue to two factors. Each volume consists of different materials, therefore different radioactive isotopes will becreated. The other reason is that the chance that a decaying isotope inside different volumes will be actuallydetected can be very different: e.g., an isotope decaying in the detector will be certainly detected, but one in thesolar panels is very unlikely to be recorded.Different energy irradiations yield in the creation of very different isotopes as nuclear cross sections varysignificantly by energy. Also, the abundance of the created nuclear isotopes is very different. For example,simulating a proton fluence of 100 000 cm − with energy of 30 MeV yields mostly to the production of Tb155and Ge68. For this fluence, about 5000 isotopes of each type are created. For an irradiation at 700 MeV, onthe other hand, the most abundant produced isotope is O15, with a number of 170 000 generated nuclei. Thesecond most abundant isotope in this case is Ga68 with 82 000 nuclei. Also, at 700 MeV the number differentcreated isotopes is 20 times higher than at 30 MeV.
3. CALCULATION OF ACTIVITY FOR TIME OF INTEREST
The input spectrum used was determined by the AP9 model of SPENVIS Monte Carlo mode with 90% confidencelevel for an orbit with height of 600 km and inclination of 40 ◦ . It should be emphasized that this represents,both for the orbit type and for radiation model used, an absolute worst case for the HERMES scenario. In this igure 2. The mass model of the HERMES satellite being irradiated by protons. case, most of the protons reach the satellite when it is passing through the SAA. In order to obtain a meanSAA flux, the proton spectrum was averaged simulating an orbit through 30 days sampled at a step of 10 s.Then the fluence of the protons inside one single pass through SAA was determined by multiplying the averagespectrum by the transit duration 60 ×
90 s (there is one SAA transit each 90 minutes; therefore the fluence in asingle passing corresponds to 60 ×
90 s of the average spectrum assuming that the proton fluence is zero outsideof SAA). In this case, an integral proton fluence of 2.5 × protons/cm at energies above 4 MeV was reached. Figure 3. Irradiation profile used to calculate activity for 700 MeV protons.
The contribution of each SAA passage is different since the created radioactive isotopes have different timeelapsed to decaying. Also each isotope has a different decay time, therefore the contribution of each SAA willvary by time. We can solve the Bateman equation for each 1 second timestep inside each SAA transit. Thereforewe choose to group the decaying times into ”time slices”. The seven SAA passages closest to the time of interestwere sampled in 21 points. The following 9 SAA transits were treated as a single irradiation. The rest of thepassages were grouped into 10 irradiation with a logarithmical integration window. In Figure 3 the describedirradiation profile can be seen.The first seven SAA transits were simulated with a Gaussian time profile sampled in 21 points. The standardeviation of the Gaussian was chosen to be 100 seconds. In order to keep computation time reasonable but stillkeep precision the 9 passes were assumed to be single irradiation. The rest of the passes were divided into 9irradiations in time. As the decays have a logarithmic nature, these 9 irradiations were grouped logarithmicallyin time as can be seen in Figure 2. Figure 4. The activation for each irradiation timeslice from 100 MeV protons for all of the volumes and for only thescintillator. Roughly half of the total satellite activation is due the scintillator.
The X-axis of Figure 4 shows time between the time of interest and the given irradiation. For each of thesepoints the decays chains are calculated and then summed up to yield the total number of isotopes. The explicitformula to determine the amount of the N th isotope in a decay chain after and elapsed time t time is: N n ( t ) = n (cid:88) i =1 N i (0) × n − (cid:89) j = i λ j × n (cid:88) j = i (cid:32) e − λ j t (cid:81) p = i,p (cid:54) = j n ( λ p − λ j (cid:33) (1)The algorithm identifies the decaying isotopes created by the proton-induced activation in the first simulation.Most of them are in ground state. Afterwards the algorithm identifies the decay chains and the decay rates ofthe corresponding daughter nuclei before the deletion of very short lived particles.The first step in this algorithm is the construction of the decay chains, which describe the decay of the isotopesidentified in the first simulation. This is implemented in C++, employing the built-in function of Geant4 toaccess the decay products of the element in question. The decay tables are taken from the ESNDF datasetincluded in the Geant4 distribution. Excited states and floating levels (same excitation energy but differentangular momentum) of the nuclei are treated individually. (An example for the importance of this: for Ta178,half life can be an order of magnitude different in case of different floating level of the nuclei.) The decay chainswere created by identifying the decay modes of the original isotopes and deducing what the daughter isotope(s)can be. The implemented algorithm looks recursively for the unstable daughters of the isotope in question. Thiswas repeated until the last element has no unstable daughters. In order to keep the number of decay chainsreasonable, daughter isotopes with branching ratios lower than a given threshold are neglected.An example of such a decay chain is:Hf156 → Yb152 → Tm152[482.320] → Tm152 → Er152[1715.400] → Er152 → Ho152[179.400] → Ho152 → Dy152[3500.000] → Dy152 → Tb152[256.930] → Tb152 → Gd152[2880.670] → Gd152 → Sm148 → Nd144.n the brackets is shown the excitation energy in keV. The number of decay chains increases by input protonenergy, all together there are more than 1.5 million chains. The Bateman equation has to be solved for each ofthese chains for every time in the irradiation profile. Moreover, the decay chains are treated independently ineach volume. Still on a standard PC (Intel Core i7-8700K CPU @ 3.70GHz ×
12, 32 GB of RAM) the algorithmruns in a few hours by setting a lower threshold for the cutoff of branching ratios to 5%.The second part of the algorithm is the exclusion of isotopes (mostly isomers) with very short half life fromthe decay chains. This is needed since to obtain the amount of each isotope after a given time, one needs to solvethe Bateman equation. In this equation, for very large decay rates ( λ ) even the double-precision floating-pointformat is not sufficient. Therefore the isotopes with a decay rate larger than 10 are erased from the decaychains: this corresponds to a half life of 21 ms. This includes only very short lived isomers (like Ta181[6.237] orW181[365.550]) which would not be abundant after the time in question passes. Therefore it does not affect theresults of the simulation. Figure 5. The contribution of different energy bands of primary protons to the activation of the satellite.
The third part of this step is to solve the Bateman equation for these decay chains. The solution of theBateman equation does not have a closed form, although it was possible to implement the solution recursivelyin C++. In order to be able to solve the equation, the decay chains were created in linear fashion. In order toobtain the activity of each radioactive isotope, their amount needs to be determined. The software determinesthis normalization factor from the provided integral photon flux.As can be seen in Figure 5 about half of the activation due to 100 MeV protons is created in the GadoliniumAluminium Gallium Garnet (GAGG) scintillator. For other energies this is different as different isotopes havedifferent energy dependence on proton activation cross section. Also what can be seen is that the activity ofthe volumes changes by time as the isotopes in the scintillator decay ”faster” than those on the other satellitestructures. The proton energy that contributes most to the overall activity in the HERMES case is about100 MeV, since protons around this energy have intermediate flux and also quite high activation cross-sectioncompared to lower energies.For different timescales different isotopes dominate the activity. For example let’s take the decay chain:Tb151 → Gd151[839.320] → Gd151 → Eu151 → Pm147 → Sm147Tb151 decays with a half-life of 18 hours. After about a day Gd151 becomes more important than Tb151 ascan be seen in Figure 6. For the simulated orbit, the activity of the whole satellite is 851 Bq after half a year inspace assuming the discussed worst-case irradiation profile (AP9 90% c.l. model, 600 km altitude, inclination of40 ◦ ). igure 6. The activity of Tb151 and Gd151 in different timeslices with respect to the time of interest. The previousdecays into the latter. For long timescales the latter is more important for shorter ones, the previous.
4. SIMULATION OF DETECTOR RESPONSE
The final result of the last algorithm is a list of isotopes sorted by activity in each volume of the mass model. Inthis simulation the most active isotopes are put in the respective volume where they were created and their decayis simulated. The energy deposited inside the Silicon Drift Detector (SDD) and the scintillator is registered. Thehistogram of the energy depositions are normalised to the isotope in question. In the end these histograms aremerged into one single histogram which represents what we would actually measure in space due to activation(assuming a perfect energy resolution). Figure 7 shows the merged histogram obtained when assuming half ayear spent in space.
Figure 7. Energy spectrum deposited inside the GAGG scintillator by the proton induced nuclear decays after half ayear in space. he deposited energy inside the scintillator can be seen in Figure 7. One can notice the β + decay peak at511 keV and also the characteristic electron capture and γ peaks. Assuming a lower detector thresholdof 20 keV the measured background would be 60 counts per second after spending half a year inspace.
It is interesting to compare this to the 800 Bq activity of the satellite. The deposited energy inside theSDD can be seen in Figure 11. The number of detections originating in the SDD is negligible compared to theones in the scintillator as the SDD is much thinner and consists of material with lower atomic number.
Figure 8. Deposited energy spectrum of the SDD by nuclear decays from activation.
5. VALIDATION OF THE FRAMEWORK
In order to ensure the validity of the results of the implemented simulation framework both direct Geant4simulations and activation measurements from literature were used. In the direct simulations, Geant4 did all thecalculations regarding the decaying nuclei. Energy depositions between 10692 minutes and 10908 minutes afterthe irradiation were recorded in this case to be compared to the results of our framework and in the experimentalresults of Sakano et al. A GAGG crystal with dimensions of 5 . × . ×
40 mm was irradiated by a protonbeam with an energy of 150 MeV. A photomultiplier tube was used to read out the activation signal. Figure 9. Deposited energy in a single GAGG crystal after 10800 minutes of irradiation by protons with an energy of150 MeV. The direct simulation performed by us matches the experimental results of Sakano et al. o be able to run the direct simulations in a reasonable timeframe, only a single GAGG crystal was irradiatedby 10 protons with an energy of 150 MeV. In Figure 9 the energy spectrum obtained in direct Geant4 simulationsdue to proton induced activation is compared to measurements obtained by Sakano et al. The results of thesimulation were convolved by a Gaussian function with changing width to include the statistical process of photondetection in the measurements. Since the used measurements only span the low energy regime (up to 200 keV)we decided to validate our framework with direct Geant4 simulations. Figure 10. Deposited energy in a single GAGG crystal after 10800 minutes of irradiation by protons with an energy of150 MeV. The direct simulation performed by us matches the results of our custom framework that we developed.
6. DISCUSSION
Within the HERMES collaboration a Geant4 based framework was developed to understand and quantify theeffect of proton induced activation. The framework utilizes direct Monte Carlo and analytical calculations tobe able to run two orders of magnitude faster and more accurate simulation, with respect to a direct MonteCarlo approach. The framework was developed in a way that the geometry of the satellite can be replaced easily.We have applied the framework to our satellite and determined energy deposition spectrum in the HERMESscintillators and in the SDD detectors as well after half a year in space. Assuming a low energy threshold of20 keV the background due to proton activation would be 60 cps which is within the required specifications. PPENDIX A. THE CONSTRUCTION OF THE DECAY CHAINS
Figure 11. The state diagram chart of the decay chain construction from the decay scheme of the mother nuclei.
ACKNOWLEDGMENTS
This project has received funding from the European Union Horizon 2020 Research and Innovation FrameworkProgramme under grant agreement HERMES-Scientific Pathfinder n. 821896 and from ASI-INAF AccordoAttuativo HERMES Technologic Pathfinder n. 2018-10-HH.0. The research has been supported by the EuropeanUnion, co-financed by the European Social Fund (Research and development activities at the E¨otv¨os Lor´andUniversity’s Campus in Szombathely, EFOP-3.6.1-16-2016-00023).
REFERENCES [1] Fiore, F. et al., “The hermes-technologic and scientific pathfinder,”
Society of Photo-Optical InstrumentationEngineers (SPIE) Conference Series , 11444–166 (2020).[2] Evangelista, Y., “The scientific payload on-board the HERMES-TP and HERMES-SP CubeSat missions,”
Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series , 11444–168 (2020).[3] Odaka, H. et al., “Modeling of proton-induced radioactivation background in hard X-ray telescopes: Geant4-based simulation and its demonstration by Hitomi’s measurement in a low Earth orbit,”
Nuclear Instrumentsand Methods in Physics Research Section A (vol. 891) , 92–105 (2018).[4] Allison, J. et al., “Recent developments in GEANT4,”
Nuclear Instruments and Methods in Physics ResearchSection A (vol. 835) , 186–225 (2016).[5] Bateman, H., “The solution of a system of differential equations occurring in the theory of radioactivetransformations,”
Proc. Cambridge Philos. Soc (Vol. 15) , 423–427 (1910).6] Sakano, M. et al., “Estimating the radiative activation characteristics of a Gd Al Ga O :Ce scintillator inlow earth orbit,” Journal of Instrumentation (vol. 9) (2014).[7] Campana, R., “The HERMES background and response simulations,”