The GeV-class gamma-rays from astronomical objects such as blazers and pulsars are observed by satellites (e.g. Fermi-LAT) or telescope arrays to study cosmic-ray acceleration mechanism, evolution of galaxies, etc. On the other hand, the cosmic ray heavy ions interacting with the nuclei in the atmosphere produce various secondary particles such as nucleons, nuclear clusters, fragments, mesons, and prompt gamma-rays. The gamma-rays from energetic projectile fragments are boosted by Doppler Effect and observed as GeV-range gamma-rays in the laboratory frame. The energetic gamma-rays from cosmic-ray heavy ions should be considered in the experiments.
Three different calculation steps, dynamic phase, evaporation and prompt gamma-decay, are necessary to simulate production of such high energy gamma-rays (Fig.1). It is also vital to simulate reactions on event-by-event basis to deduce the excitation energy of each fragment. A number of simulation codes have been developed for prediction of cosmic ray reactions; however, very few can simulate event-by-event prompt gamma-ray emission[1-4] because it is necessary to determine the isotopic species produced in each reaction and to calculate the internal transition intensity between discrete levels based on the nuclear structure data [5-7].
The general purpose radiation transport code PHITS (Particle and Heavy Ion Transport code System) incorporates a heavy ion reaction model JQMD (JAERI Quantum Molecular Dynamics Model) [8,9], statistical decay model GEM (Generalized Evaporation Model) , and prompt gamma-ray production model EBITEM (ENSDF-Based Isomeric Transition and isomEr production Model)  and therefore PHITS is capable of simulating prompt gamma-rays from heavy ion fragments.
Prompt gamma-ray emission was simulated in typical conditions such as 10 AGeV Oxygen interactions with nitrogen and 1 AGeV boron interactions with nitrogen. The calculated gamma-ray energy spectra showed sharp peaks corresponding to particular internal transition paths.
In this study, cosmic ray transport and interaction simulation is performed by using PHITS. PHITS is a general-purpose radiation transport simulation code which can simulate transport and reactions of nearly all kinds of particles from thermal energies to 1 TeV. In PHITS, the 3 reaction phases (dynamic phase, evaporation and prompt gamma-decay) were simulated by the JQMD, GEM, and EBITEM, respectively. Because all these models are event generators, secondary particle emission was calculated in event-by-event basis. The three models are briefly explained below.
Dynamic phase model JQMD
JQMD is one of the reaction event generators based on the quantum molecular dynamics (QMD) approach , which is often used for simulation of nucleus-nucleus reaction dynamic phase. In comparison with other QMD models, JQMD can simulate various perspectives, such as secondary particle production and fragment yields, by using fewer parameters.
In JQMD, nucleons follow a Lorentz-covariant equation of motion. The long range interaction between nucleons is described by the potential,
where A is one of the Skyrme force constant (-219.4MeV), ρ is the saturation density (0.168 fm ), τ is another Skyrme force constant (4/3), B is the other Skyrme force constant (165.3 MeV), ρ is the local nucleon density at the position of i-th particle, is the spatial coordinate of the centroid of the i-th particle, c is charge of i-th particle, e is elementary charge, L is the width of wave packet (2 fm ), C is the symmetry force constant (25 MeV) and ρ is the wave function overlap of i-th particle and j-th particle. Strong repulsion attributed to the hard core is not included in this potential. Instead, the nucleon pairs close to each other undergo elastic or inelastic scattering based on their cross sections. Production of baryon resonance particles such as Δand N* as well as absorption of pions by nucleons are also considered.
After 150fm/c of time evolution, nucleons close to each other are bound to form clusters. Mass, charge, kinetic energy, and excitation energy are then calculated by summing the contribution from the constituent nucleons. More details on JQMD are available elsewhere [7,8].
Evaporation model GEM
GEM is one of the event generators to simulate statistical decay of excited nuclei. Nuclei undergo sequential particle emission based on the Weisskopf-Ewing approach [13,14] in GEM. Particle emission probability Γis calculated by,
where subscript i denotes initial nucleus, subscript j denotes ejectile particle, and subscript d denotes daughter nucleus. S is ejectile spin, m is ejectile mass,ρ is level density of daughter nucleus calculated by Fermi gas model or by the model proposed by Gilbert and Cameron , Q is Q value, V is Coulomb potential between the projectile and the ejectile, R , c , c , and b are model parameters given in , respectively.
The initial state of statistical decay is specified by the clustering phase of JQMD. The mass, charge, kinetic energy, and excitation energy are passed to GEM.
Prompt gamma-decay model EBITEM
A theoretical model to simulate gamma de-excitation of excited nuclei, EBITEM (ENSDF-Based Isomeric Transition and isomEr production Model), is based on the Evaluated Nuclear Structure Data File (ENSDF), and theories . When excitation energy became smaller than the nucleon separation energy, evaporation by GEM was stopped and the residue, characterized by mass, charge, angular momentum and excitation energy, starts gamma de-excitation. If the excitation energy was smaller than 3 MeV or the nuclear mass was smaller than 40 amu, the level structure taken from ENSDF was employed. Otherwise, the level density was calculated by the Gilbert-Cameron formula expressed in Eq. (3),
where U is excitation energy, a is level density parameter, J is angular momentum, and σ is spin cut-off parameter.
The probability of transition from one state to another state was taken from ENSDF if the nucleus was at an excitation state given in ENSDF; otherwise it was calculated by the single-particle model. The transition probability was calculated as a function of total angular momentum transfer and gamma-ray energy. By sampling the de-excitation path at random, the next level characterized by excitation energy and total angular momentum was determined. By repeating this process corresponding to de-excitation, the nuclei reach either the ground state or the metastable states.
The detailed description of the theory and validation of EBITEM is provided elsewhere .
3. Results and Discussion
Some case studies on the application of PHITS are discussed.
Fig. 2 shows particle trajectories after an interaction of C with N at 100 AGeV in the atmosphere. The air density averaged from the ground level to 10 km in the sky was adopted as the air density in this calculation. Various particles such as muons, pions, photons, neutrons, and neutrinos are seen as particle tracks. High energy gamma-rays are also observed in this event. In the laboratory frame, two photons with energies 410 MeV and 646 MeV were detected in coincidence and their emission angle was shifted by 0.17 degrees from each other. This gamma emission event is attributed to internal transition of C emitting 3.12 MeV (7.56MeV - 4.44MeV) and 4.44 MeV (4.44MeV- Ground) gamma-rays. Gamma-ray energy in the laboratory frame is shifted by Doppler effect taking into account for the projectile energy (100 AGeV) and the angle between gamma emission and projectile momentum.
In this way, PHITS can be used to reproduce one event starting from one heavy ion interacting with the target nucleus.
Fig. 3 shows energy spectrum of gamma-rays from the fragments produced by 10 AGeV B interacting with nitrogen. The gamma-ray energies are calculated in the rest frame of the projectile.
Gamma-ray peaks attributed to various fragment species are seen in the spectrum. It is possible to simulate correlation between the fragment and its gamma-ray. In reality, this peak structure is blurred because the kinetic energy of fragments is distributed and gamma-ray energies are shifted by Doppler Effect.
In this study, it is shown that PHITS is capable of simulating interactions of cosmic-ray heavy ions. In particular, for predicting prompt gamma-rays from the fragments of cosmic-ray heavy ions, combination of a non-equilibrium reaction model, a statistical decay model and a prompt-gamma decay model is essential.
In addition, needful functions and data such as transformation between the frames, event-by-event correlation and nuclear level structure are implemented to help users to obtain desired quantities in useful conditions.
These facts indicate that PHITS is one of the most robust tools to plan and analyze cosmic ray experiments. In particular, PHITS offers a unique capability to simulate prompt gamma-rays from cosmic ray heavy ion fragments.
We express our gratitude to the operation team of the Center for Computational Science and E-systems (CCSE) at the Japan Atomic Energy Agency. The simulations reported in this paper were executed on a PC-cluster system located at CCSE. This work was supported by JSPS Grants-in-Aid for Scientific Research (Number 26790072).