© F. Ballarini et al., Published by EDP Sciences, 2024

1. Introduction

FLUKA [1] is a Monte Carlo code capable of simulating the transport of all elementary particles and ions in complex geometries and materials. Possible beams include all flavourless and strange neutral and charged hadrons that do not undergo strong decays, ions as light as deuterons and as heavy as uranium, muons, electrons, positrons and photons. Hadronic resonances can be produced and can decay during nuclear interactions. The energy range supported by FLUKA is from 1 keV (keV/ n for ions) to 10²¹ eV for most particles. There are two notable exceptions: photons and neutrons can be transported down to 100 eV and 10⁻⁵ eV respectively.

FLUKA is used in various fields, including medical physics, cosmic ray studies, shielding, dosimetry, radiation protection, calorimetry, detector simulations, accelerator-driven systems and more. The user community has several thousand members worldwide and is continuing to grow.

This paper is structured to give a global overview of the code, focusing on the developments made over the last decade. For the most part, the reader is referred to the relevant literature for features that have been stable for a long time. Some of these features are recalled in this paper when their description is only available in hard-to-access references (typically hardback conference proceedings with no electronic counterpart). New additions and recent developments are described in some detail below. The history of FLUKA is presented in Section 2, followed by a general review of the fundamental interaction models contained in the code in Section 3. Recent developments in hadronic models are summarized in Section 4. Section 5 is dedicated to new point-wise fully-correlated neutron transport, Section 6 deals with recent developments in the EM part, while Section 7 covers novelties about charged particle transport. Examples of recent and advanced FLUKA applications are presented in Section 8. Finally, Section 9 briefly describes some technical improvements.

2. History of FLUKA

The “first” FLUKA was born in 1970, when Johannes Ranft, then working at RHEL (UK), developed the first analogue hadron cascade code, named “FLUKA” after “FLUktuierende KAskade”. In the same year, J. Ranft got a position at the Karl-Marx-Universität, Leipzig (GDR), and until the dissolution of the GDR, the development of FLUKA, in particular the development of the hadronic models, took place mainly in Leipzig.

*8pc

Starting from these early pioneering efforts, three different generations of “FLUKA” codes can be distinguished over the years, which can be roughly identified as the FLUKA of the 1970s (main authors J. Ranft and J. Routti), the FLUKA of the 1980s (P. Aarnio, A. Fassò, H.-J. Moehring, J. Ranft, G.R. Stevenson), and, from 1989 to 2019, the modern FLUKA (A. Fassò, A. Ferrari, J. Ranft and P.R. Sala). This paper describes the fourth generation of the FLUKA code which is a further evolution of the modern FLUKA, and which is still authored by A. Fassò, A. Ferrari, and P.R. Sala. Johannes Ranft passed away in 2018. This paper is dedicated to his memory.

Each new “generation” was not just an improvement on the existing program. It was a quantum leap in code, physics, design and objectives. The same name “FLUKA” has been retained as a reminder of this historical development – mainly as a tribute to J. Ranft, who was involved as an author and mentor from the beginning until his death – but the third and fourth generation codes are completely different and much more powerful than the previous generations. In fact, the modern FLUKA is a completely new model, with many capabilities that have not been present in any of the previous FLUKA generations.

The modern FLUKA [1–6] was born at the INFN in 1989. It was mainly the work of Alfredo Ferrari and Paola Sala, with important contributions from Alberto Fassò and Johannes Ranft. Its development was carried out entirely within the INFN until 2003 when CERN joined it. This was done under an agreement between the two institutions and the authors. This agreement came to an end in 2019. The development and maintenance of FLUKA is continued by the surviving legitimate authors, Alberto Fassò, Alfredo Ferrari, Paola Sala, and their collaborators in the framework of the FLUKA Collaboration. The code is publicly available at http://www.fluka.org. At present there are worldwide more than six thousand active users for the last FLUKA generation.

3. The FLUKA models: a general description

The following is a summary of the various models used in FLUKA to describe nuclear and atomic interactions and transport. The nuclear interaction models are by far the most complex part of the code, a scheme showing the relations between incident particles/energies and the various algorithms used for that purpose is shown in Figure 1.

3.1. Hadron-Nucleon interactions

An essential prerequisite for a successful description of hadron–nucleus elastic and inelastic scattering is an accurate hadron–nucleon interaction model. Details on the different approaches for different projectiles and energies used in FLUKA for the description of hadron–nucleon interactions can be found in [7–9]. For nucleons, pions, kaons and anti-nucleons elastic, charge exchange, strangeness exchange and annihilation on nucleons, a mixture of fits to available experimental data and phase shift analyses is used up to a few GeV. Non-elastic (particle production) hadron–nucleon interactions up to a few GeV’s are based on the production and subsequent decay of hadronic resonances through resonant and non-resonant channels. The relative probabilities of the different channels are calculated according to available experimental and prescriptive data and then extended to unmeasured channels using isospin symmetry considerations [8].

At higher energies, FLUKA relies on the hadronization of hadronic strings to describe the particle-producing hadron–nucleon interactions. In fact, in QCD, the colour field acting between the quarks is carried by gluons, the vector bosons of the strong interaction, which are themselves “coloured”. Thus, the

characteristic feature of gluons (and QCD) is their strong self-interaction, which unfortunately prevents the use of the QCD

Lagrangian in the non-perturbative, soft p_T regime, which dominates the bulk of hadron interactions. If we imagine that the quarks are held together by colour lines of force, the gluon-gluon

interaction will pull them together in the form of a tube or a string. Due to quark confinement, the energy required to “stretch” such a string becomes increasingly large, until it is sufficient to

materialise a quark-antiquark pair from the vacuum, and the string breaks into two shorter ones, with still (anti)quarks at either end.

Several quark-string models exist; in FLUKA, the hadron–nucleon interactions are described within the framework of the Dual Parton Model (DPM) [8, 10]. The Dual Parton Model is a special quark/parton string model. It gives reliable results up to tens of TeV. In the DPM, hadron–hadron interactions lead to the production of two or more QCD colour strings. From these, hadrons must be produced. These hadronic strings are then hadronized using a native FLUKA hadronization algorithm.

Some examples of the performances of the FLUKA DPM implementation can be found in the next sections. In particular, Section 4.3 contains a few comparisons with experimental data, while describing recent improvements to the hadronization algorithm.

3.2. Hadron–Nucleus interactions: PEANUT

In the second FLUKA generation, the intermediate (E <  5 GeV) energy model was based on the production and decay of one or more resonances to model the primary hadron–nucleon interactions in conjunction with a parameterized intra-nuclear cascade. Similarly, the high-energy hadron–nucleus event generator was based on parameterizations for the number of primary collisions, as well as for the intra-nuclear cascade. Both models did not include any evaporation stage and were unable to produce low-energy particles and/or residual nuclei. Despite many improvements in the parameterizations, in the kinematics and the treatment of the nuclear effects [7, 13, 14], those models became increasingly inaccurate, particularly at energies below 1 GeV. Therefore, since 1991, these two models have been gradually replaced by the PEANUT model [2, 8, 9, 15, 16] (PreEquilibrium Approach to NUclear Thermalization). At that time, most Monte Carlo codes still lacked a correct treatment of the nuclear thermalisation process, including a pre-equilibrium phase, with the exception of LAHET, which included a specialised algorithm [17], and GEANT3, which implemented an older version of PEANUT [15].

PEANUT was originally designed to simulate the interactions of nucleons, pions and γ-rays with nuclei from approx. 2 GeV down to the reaction threshold (or 20 MeV for neutrons). As far as photo-nuclear reactions are concerned, details can be found in [18].

Since then, the model has undergone a long evolution and is now able to treat all stable¹ hadrons and anti-hadrons, and real and virtual photons inelastic interactions with nuclei, from the threshold up to 10–20 TeV, laboratory energy. Starting with Fluka2023.3, the old model, which was still used in some niche situations (for photon-nucleus interactions in the intermediate energy range), is no longer distributed.

As explained in Section 3.1, up to a few GeV, inelastic hadron–nucleon collisions are described by resonance production and decay. A model based on the Dual Parton Model [10] (DPM) takes over at higher energies. It is possible to extend DPM to hadron–nucleus collisions [10, 19] using the Glauber-Gribov approach [20–23]. Furthermore, DPM provides a theoretical framework for the description of hadron diffractive scattering in both hadron–hadron and hadron–nucleus collisions.

The hadron–nucleus reaction mechanism is therefore modelled in PEANUT by the Glauber-Gribov cascade, followed by the Generalised Intra-Nuclear Cascade (GINC), smoothly coupled to statistical (exciton) pre-equilibrium emission [24, 25]. The pre-equilibrium stage is in turn followed by evaporation or fragmentation of the excited remainder and finally by gamma de-excitation.

Some of these steps are omitted, depending on the primary hadron energy. The onset of the Glauber-Gribov cascade with the resulting multiple primary collisions occurs at a few GeV. Lower energy hadrons go directly to the GINC stage or, in the case of protons and neutrons below a few tens of MeV, to the pre-equilibrium stage.

An obvious requirement for all models involving an intra-nuclear cascade is that the wavelength associated with the motion of the hadron must be much shorter than the mean free path of the hadron inside the target nucleus, and also much shorter than the mean distance between two neighbouring nucleons. It is easy to show that such an assumption breaks down for energies below a few hundred MeV (see [26]) unless special care is taken in the introduction of quantum effects and in the treatment of the effect of the nuclear mean field. In any case, at energies well below 100 MeV INC models are definitely ruled out and one has to resort to pre-equilibrium models to properly describe particle emission. Indeed, PEANUT includes a similar model as a second stage after the INC part for sufficiently large projectile energies (E >  50 MeV) or directly as the first stage at lower energies.

At the end of the INC stage, three further steps are performed before completing the reaction:

• Pre-equilibrium [24, 25] stage, whenever all excited nucleons are below a given energy threshold (typically 30–50 MeV);

• Evaporation stage, whenever the pre-equilibrium stage is finished, and the system can be assumed to be equilibrated;

• Final de-excitation stage when the excitation energy is below the threshold for particle emission and it is dissipated through photon emission.

A description of the physics involved in all stages can be found in [2, 8, 9, 16], together with several examples of its performances when compared with experimental data. One example is presented in Figure 9 where the computed double differential spectra for the reactions ⁹⁰Zr(p,xn) and ⁹⁰Zr(p,xp) at 160 MeV are compared with experimental data [27, 28]. A summary of the main steps of our approach is provided in the following sections, giving more emphasis on recent and/or unpublished developments.

3.2.1. Glauber-Gribov cascade

The Glauber formalism [41, 42] allows computing the scattering amplitude as well as all relevant cross-sections for hadron–nucleus interactions, using the knowledge of elementary hadron–nucleon scattering and of the nuclear ground state only. The equations used in PEANUT can be found in references [8, 16]. Here we only remind that the Glauber theory allows us to compute the total, elastic, non-elastic (aka absorption or particle production), and inelastic (incoherent elastic or quasi-elastic) scattering cross-sections starting from the hadron–nucleon scattering amplitude and the target nucleus density distribution.

Expressing the hadron–nucleon scattering amplitude by means of the S-matrix approach in impact parameter space:

$$\begin{array}{c}\hfill S_{\mathit{hN}}(\mathrm{b⃗},s)={\mathrm{e}}^{\mathrm{i}{\chi }_{\mathit{hN}}(\mathrm{b⃗},s)}={\eta }_{\mathit{hN}}(\mathrm{b⃗},s){\mathrm{e}}^{2\mathrm{i}{\delta }_{\mathit{hN}}(\mathrm{b⃗},s)}\end{array}$$

the Glauber cross-sections for inelastic and elastic scattering can be written as (Ψ_i is the initial state global wave function of the target nucleus):

$$\begin{array}{cc}\hfill {\sigma }_{hA\mathrm{\Sigma }f}(s)& \equiv \sum _f{\sigma }_{hAfi}(s)\hfill \\ \hfill & =\int {\mathrm{d}}^2\mathrm{b⃗}\int {\mathrm{d}}^3\mathrm{u⃗}{|{\mathrm{\Psi }}_i(\mathrm{u⃗})|}^2{\left|[1-\prod _{j=1}^A{S}_{\mathit{hN}}(\mathrm{b⃗}-{\mathrm{r⃗}}_{j\perp },s)]\right|}^2\hfill \end{array}$$(1) $$\begin{array}{cc}\hfill {\sigma }_{hAel}(s)& =\int {\mathrm{d}}^2\mathrm{b⃗}{|\int {\mathrm{d}}^3\mathrm{u⃗}{|{\mathrm{\Psi }}_i(\mathrm{u⃗})|}^2[1-\prod _{j=1}^A{S}_{\mathit{hN}}(\mathrm{b⃗}-{\mathrm{r⃗}}_{j\perp },s)]|}^2\hfill \end{array}$$(2)

σ_{hA el} and σ_{hA Σf} are the cross-sections for elastic scattering and inelastic scattering to excited states. They do not include contributions from individual hN non-elastic (particle production) interactions.

The total cross-section can be expressed as:

$$\begin{array}{cc}\hfill {\sigma }_{hAT}(s)& =2\int {\mathrm{d}}^2\mathrm{b⃗}\int {\mathrm{d}}^3\mathrm{u⃗}\hfill \\ \hfill & |{\mathrm{\Psi }}_i{(\mathrm{u⃗})|}^2[1-\prod _{j=1}^A\mathrm{Re}{S}_{\mathit{hN}}(\mathrm{b⃗}-{\mathrm{r⃗}}_{j\perp },s)].\hfill \end{array}$$(3)

The quasi-elastic cross-section (incoherent–elastic) is:

$$\begin{array}{c}\hfill {\sigma }_{hAqe}(s)\equiv {\sigma }_{hA\mathrm{\Sigma }f}(s)-{\sigma }_{hAel}(s).\end{array}$$(4)

The absorption (particle production) cross-section can be written as:

$$\begin{array}{cc}\hfill {\sigma }_{hAabs}(s)& \equiv {\sigma }_{hAT}(s)-{\sigma }_{hAel}(s)-{\sigma }_{hAqe}(s)=\hfill \\ \hfill & =\int {\mathrm{d}}^2\mathrm{b⃗}\int {\mathrm{d}}^3\mathrm{u⃗}{|{\mathrm{\Psi }}_i(\mathrm{u⃗})|}^2\hfill \\ \hfill & \{1-\prod _{j=1}^A\{1-\left[1-|S_{\mathit{hN}}(\mathrm{b⃗}-{\mathrm{r⃗}}_{j\perp },s){|}^2\right]\}\}\hfill \\ \hfill & \equiv \int {\mathrm{d}}^2\mathrm{b⃗}{\mu }_{hAabs}(\mathrm{b⃗},s).\hfill \end{array}$$(5)

An example of the application of the Glauber-Gribov cascade plus the Generalized Intranuclear Cascade (next section) is presented in Figure 4 for 31 GeV/c protons on carbon, compared with experimental data from the NA61 experiment [12]. Other examples can be found in Section 4.

3.2.2. Generalized IntraNuclear cascade

Hadron–nucleus non–elastic interactions are described in our model in the framework of the IntraNuclear Cascade (INC) model, possibly applied after a Glauber-Gribov cascade if the projectile energy is in excess of a few GeV. The concept of the INC model was developed at the very beginning (the original ideas go back to the end of the 50s [43, 44]) of the history of energetic nuclear interaction modelling but it is still valid and in some energy ranges it is the only available choice. The model is intrinsically a Monte Carlo model, well suited for numerical applications, while no closed analytical expression can be derived without severe approximations. In the energy range going from the pion production threshold (approx. 290 MeV for a free nucleon, down to 200 MeV for nucleons in nuclei because of the Fermi motion) to high energies, INC models are widespread tools to describe hadron–nucleus interactions. At lower energies, a variety of pre-equilibrium models work very well, with foundations in physics that become increasingly more robust as the energy decreases, in contrast to those of INC.

In our INC stage, all particles are transported along semi-classical trajectories that are subject to curvature in the Coulomb and nuclear potentials. In this way, refraction and reflection at the nuclear surface are taken into account, and the Coulomb effects are properly described. The nuclear potential depths and shapes used for neutral and charged hadrons have been found to have a strong influence on the results: more details can be found in [2, 8]. A discussion of the assumptions underlying the INC approach and related references can be found in [26].

The basic assumptions of the FLUKA INC model can be summarized as follows:

1. The target nucleus is modelled as a few concentric spheres of different densities, adapted to reproduce the experimental density distribution of nuclei;

2. Nucleons in the target nucleus move according to a density-dependent potential derived from the local Fermi energy;

3. The projectile impact parameter is selected with a constant probability over the geometrical cross-section area;

4. Hadrons propagate like free particles in the nuclear medium, with interaction probability per unit length given by free space cross-sections, properly averaged over the Fermi motion of the target nucleons, times the local nuclear density; an exception is represented by photons and pions where the effect of nuclear medium on the free cross-section in the Δ region are taken into account according to [45–47];

5. The particle motion is formulated in a semi-classical way. It is subject to a mean nuclear potential, and to the Coulomb potential, which must be added to the free kinetic energy of the particle as it travels through the nucleus. The depth of the nuclear field and the radial profiles depend on the particle and the energy;

6. The effect of the nuclear mean field on the particle motion produces curved trajectories in a semi-classical approach, according to energy and momentum conservation, depending on the model. The curvature effects induced by the nuclear mean field are usually referred to as refraction and reflection effects;

7. Interactions occur like in free space in the centre of mass System of the two colliding hadrons. Because of Fermi motion, obviously, the laboratory frame will not coincide with the frame where the target nucleon is at rest, but suitable Lorentz boosts are applied to transform the secondary particles back to the laboratory frame;

8. Interactions occur in a completely incoherent and uncorrelated way. No coherence or diffractive effects among successive interactions are taken into account. No multi-body or cluster processes are included, with the exception of muon, pion, and photon absorption;

9. Quantum effects include Pauli blocking, nucleon–nucleon correlations, fermion anti-symmetrization and coherence length (for elastic and charge exchange elementary interactions), and formation zone (for inelastic ones) effects;

10. Secondaries are treated exactly like primary particles, with the only difference being that they start their trajectory already inside the nucleus;

11. Binding energies are calculated from recent nuclear mass databases [48]. They evolve with the evolution of the reaction and account at every emission stage for the proper Q.

It should be emphasized that the typical problems of INC codes with binding energies and reaction Q’s are completely solved in our approach, due to the use of dynamically changing binding energies.

In the following, we will summarize the implementation of the assumptions introduced in the FLUKA model, such as the required quantum effects and some important features related to nuclear physics, such as coalescence. We will also briefly discuss some special treatment of pions, kaons, antinucleons and negative muon capture. Other stages of the interaction that occur after the INC: pre-equilibrium, nuclear evaporation, fission and finally gamma de-excitation will also be described. These aspects are essential for a reliable prediction of the abundance of residual nuclei after a reaction.

3.2.3. Nuclear density and Fermi motion

In both stages, INC and pre-equilibrium, the nucleus is modelled as a sphere with density given by a symmetrized Woods-Saxon [49] shape for A >  16, and by a harmonic oscillator shell model for light isotopes (see [50]). A standard local density approximation Fermi momentum distribution is implemented in PEANUT in order to compute the nucleon mean field:

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathit{dk}}=\frac{{|k|}^2}{2{\pi }^2}\end{array}$$(6)

for k up to a local Fermi momentum k_F(r) given by

$$\begin{array}{c}\hfill k_{\mathrm{F}}(r)={\left(\frac{3{\pi }^2}{2}\rho (r)\right)}^{\frac{1}{3}}.\end{array}$$(7)

Densities, and therefore Fermi momenta, are evaluated separately for protons and neutrons.

3.2.4. Quantum effects

The main quantum effects considered in the INC and pre-equilibrium stages of FLUKA include Formation Zone, Coherence Length and Nucleon–Nucleon anti-symmetrization. They are essential in reducing hadron re-interaction rates with respect to naive expectations based on free hadron–nucleon cross-sections.

3.2.4.1. Pauli blocking.

This is a well-known effect due to the fermionic nature of nucleons: for each energy level there can be at most four nucleons, namely two protons and two neutrons, each pair with opposite spin. Thus, interactions that would result in a nucleon occupying an already-filled level are not allowed.

3.2.4.2. Formation zone.

A simple treatment of the INC at high energies leads to an overestimation of the computed particle yields whenever the incident energy is higher than a few GeV. The experimental evidence suggests that a mechanism exists that limits the number of re-interactions of energetic (fast) particles. The physical mechanism at work is the so-called “formation zone”. This concept has a sound theoretical foundation [51, 52] and is analogous to the Landau–Pomeranchuk–Migdal effect [53, 54], which has been carefully verified experimentally. However, an exact derivation would require knowledge of the generated hadron wave functions. In any case, it is not clear whether it should be applied to the resonances produced by the hadronisation rather than to their final decay products.

An approximate formula can be derived if one considers that hadrons are composite objects and that the typical time of the strong interaction is of the order of 1 fm/c. After emerging from an inelastic interaction, it takes some time for hadrons to “materialise” and be able to undergo further interactions. This time can be estimated qualitatively by considering a secondary particle with mass M and momentum components p_∥ and p_T with respect to the original projectile direction. After carrying out a Lorentz transformation to the frame in which the longitudinal momentum is zero, the uncertainty principle provides:

$$\begin{array}{c}\hfill \overline{t}=\mathrm{\Delta }t\approx \frac{ħ}{E_T}=\frac{ħ}{\sqrt{p_T^2+M^2}}.\end{array}$$(8)

This time interval can be translated into the particle proper time, τ, and the lab frame time, t, taking care of the Lorentz dilation among the various frames:

$$\begin{array}{cc}\hfill \tau & =\frac{M}{E_T}\overline{t}=\frac{ħM}{p_T^2+M^2}\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill t_{\mathit{lab}}& =\frac{E_{\mathrm{lab}}}{E_T}\overline{t}=\frac{E_{\mathrm{lab}}}{M}\tau =\frac{ħE_{\mathrm{lab}}}{p_T^2+M^2}.\hfill \end{array}$$(10)

The time interval in the lab frame can be also expressed as a function of the particle rapidity, y, by:

$$\begin{array}{c}\hfill t_{\mathrm{lab}}=\overline{t}coshy=\frac{ħ}{\sqrt{p_T^2+M^2}}coshy.\end{array}$$(11)

If such an interaction has to take place inside a nucleus, the condition for having (possible) re-interactions of our secondary particle can be expressed by:

$$\begin{array}{c}\hfill v·t\le R_A\approx r_0{A}^{\frac{1}{3}}.\end{array}$$(12)

Such an equation can be used to define a critical rapidity above which the particles have no chance of interacting, and in fact “materialise” outside the nucleus. Inserting typical values for pions into equation (11), E_T ≈ 0.3 GeV, E ≈ 10 GeV, exemplifies that the formation zone can easily exceed the radii of heavy nuclei. All this derivation should be considered as pragmatic theory-inspired and semi-quantitative approach. It does not attempt to be a solid theoretical derivation.

Among all the other uncertainties, it is not clear whether the particle entering the formulae must be the final one, or rather the “mother” resonance. Nevertheless, the concept is very powerful, and excellent agreement with experimental data has been obtained for hadron–nucleus and nucleus-nucleus interactions using a universal scale factor for the formation time as a free parameter [16, 52, 55]. Indeed, in PEANUT, the adopted universal parameter in front of equation (11) is just 1.

3.2.4.3. Coherence length.

Another critical issue is the “coherence” length after elastic or charge exchange hadron–nucleon reactions. In analogy to the formation zone concept, such interactions cannot be localized better than the positional uncertainty associated with the four-momentum transfer of the collision. Re-interactions occurring at distances shorter than the coherence length would interfere with each other and, in any case, cannot be treated as independent interactions on individual nucleons. The coherence length concept is applied in PEANUT to the secondaries generated in inelastic, charge exchange, or, more generally, quasi-two-body interactions, with the following recipe: given a two-body interaction with four-momentum transfer q = p_1i − p_1f, between a particle 1i and a nucleon 2i with final particles (which could be resonances) 1f and 2f the energy transfer seen in a frame where the particle 2 is at rest is given by:

$$\begin{array}{c}\hfill \mathrm{\Delta }{E}_2={\nu }_2=\frac{q·p_{2i}}{m_2}.\end{array}$$(13)

From the uncertainty principle this ΔE₂ corresponds to an indetermination in the proper time given by: Δτ ⋅ ΔE₂ = ℏ, which boosted to the lab frame gives a coherence length:

$$\begin{array}{c}\hfill \mathrm{\Delta }{x}_{\mathrm{lab}}=\frac{p_{2\mathrm{lab}}}{m_2}·\mathrm{\Delta }\tau =\frac{p_{2\mathrm{lab}}}{m_2}\frac{ħ}{{\nu }_2}.\end{array}$$(14)

3.2.4.4. Nucleon–Nucleon anti-symmetrization.

An important quantity which can be computed in the Fermi gas model is the so-called correlation function. Due to the anti-symmetrization of the fermion’s wave function, given a nucleon at position $\overrightarrow{r}$ in a nucleus with density ρ₀, the probability of finding another alike nucleon at position ${\overrightarrow{r}}^{\prime }$ is decreased for small values of the distance $r=|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|$ by a factor

$$\begin{array}{c}\hfill g(x)=1-\frac{1}{2}{\left[\frac{3}{x^2}(\frac{sinx}{x}-cosx)\right]}^2\end{array}$$(15)

where x = K_F ⋅ r, K_F is the (local) Fermi momentum, and the factor $\frac{1}{2}$ in front of the parenthesis accounts for the two possible spin orientations.

3.2.5. Other nuclear effects

3.2.5.1. Coalescence.

Emission of energetic light fragments, d, t, ³He, and α through the coalescence mechanism is included all along the PEANUT reaction chain. Coalescence allows to reproduce the high-energy tail of the light fragment spectra [56, 57]. An example of the deuteron emission spectrum from protons reactions on Aluminum at 14.6 GeV/c is shown in Figure 2, together with pions, protons, and kaons differential production cross-sections, all compared with experimental data [11]. Double differential spectra for positive pion production for the same reactions are shown in Figure 3.

3.2.5.2. Hard core.

The strong nuclear interaction between nucleons includes a short-range repulsive component, known as hard core. The radius of this repulsive core is of the order of 0.5–1 fm. As a result, a hadron that underwent an h–N interaction cannot re-interact with another nucleon before exiting from the hardcore zone.

3.2.6. Special cases addressed by the PEANUT model

3.2.6.1. Pions.

The modelling of π nucleus interactions in PEANUT [2, 7, 8] in the Δ resonance region follows the theoretical approach of Oset and collaborators [45], taking into account:

• The resonant nature of the π–N interaction, mostly dominated by the Δ(1232);

• The effect of the nuclear medium on the π–N interaction;

• The possibility of absorption (both s-wave and p-wave) on two or more nucleons;

• The resonant nature of the pion-nucleus potential, which is rapidly varying with the pion energy.

3.2.6.2. Kaons.

All nuclear effects already present in PEANUT are exploited in the FLUKA implementation. The initial kaon–nucleon interaction is described by means of phase shift analysis [58–60]. Isospin relations are used to link different charge states. Mass differences are taken into account. K⁺ and K⁰ interact weakly with nucleons, only the elastic and charge-exchange channels are open up to 800 MeV/c (π production threshold). Thus their interactions are easily described and they are a good probe of the nuclear medium, in particular the Fermi distribution. Examples of quasi-elastic K⁺ scattering from nuclei can be found in [5]. Conversely, a variety of hyperons can be produced in K⁻ and ${\overline{K}}^0$–N interactions, with the Σπ and Λπ channels already open at rest.

3.2.6.3. Negative muon capture.

Negative muons coming to rest are captured on atomic orbitals and can either decay or be eventually absorbed by the nucleus, via the weak interaction process μ⁻ + p → ν_μ + n embedded however in the nuclear environment.

Since this process is slow, it is assumed to occur after the muon has cascaded down to the (muonic) atomic K-shell. The muonic X-ray de-excitation in this phase is taken into account in FLUKA. The slight reduction in the muon decay rate due to binding is also taken into account.

Muon capture rates are computed following the Goulard–Primakoff [61] formula. This formula involves a calculated nuclear effective charge, and Pauli blocking terms that are determined by fits to experimental data. The parameters adopted in FLUKA are those suggested in [62]. In addition to the single-nucleon reaction described above, absorption on a nucleon pair is also considered, as suggested for example in [63]. Single-nucleon absorption alone would lead to a dramatic underestimation of both the high-energy tail of the emitted nucleon spectra and the emitted neutron multiplicity. The two-body mechanism is in fact the first approximation to a more coherent absorption description that allows for different energy-sharing mechanisms. Spectra with different relative probabilities for two-nucleon absorption are shown in Figure 5 for the μ⁻ absorption on lead, together with the experimental data. In FLUKA, a probability of 20% for two-body absorption is assumed. Work is in progress for a better estimation of the two-body probability as a function of the nuclear mass. Kinematics for two-nucleon absorption is then calculated according to the phase space. The absorption position in the nucleus is sampled as a function of the nuclear density for single-nucleon interactions and as a function of the square of the nuclear density for the two-nucleon mechanism. The model adopted does not attempt to describe any additional features that are driven by details of the nuclear structure, as this would go well beyond the scope of a general-purpose interaction and transport code. The results in terms of the capture ratio and the spectra of the emitted particles are more than satisfactory, as shown in Figure 5 for lead and calcium.

3.2.6.4. Antinucleons.

Low-energy anti-nucleons colliding with nucleons undergo elastic, charge exchange, and annihilation processes. Bendiscioli and Kharzeev [33] wrote a complete review of available experimental data on anti-nucleon interactions on nucleons and nuclei. They also provide parametrizations for the different cross-sections depending on the energy range. These parametrizations can be joined together for a smooth Monte Carlo sampling. The result for the antiproton-proton system is shown in Figure 6 together with experimental data from [33]. Particle multiplicities and spectra following annihilation are calculated through the production and decay of two or more resonant intermediate states whose branchings are adjusted to reproduce experimental multiplicities.

Anti-nucleon interactions on nuclei occur preferentially on peripheral nucleons. Reaction cross-sections are roughly proportional to the square of the nuclear radius. In FLUKA, the annihilation depth is chosen as a function of the atomic number, the interaction being more and more peripheral as the nuclear mass increases. Annihilation is the dominant process up to about 1 GeV/c. Examples of experimental [33] and simulated reaction cross-sections are shown in Figure 7.

Final particle spectra and multiplicity are influenced by final state interactions, in particular by pion absorption and charge exchange. For instance, the average charged pion multiplicity following $\overline{p}$ annihilation at rest decreases from 3.1 on deuterium to 2.4 on uranium. Examples of PEANUT results are shown in Figure 8. In the case of compounds, the relative annihilation probabilities are calculated following [64] for hydrogenated compounds, and [65] for other compounds. Satisfactory agreement between FLUKA results and experimental data has been published by the AEGIS collaboration [66].

3.2.7. Pre-equilibrium stage

The INC step goes on until all nucleons are below 50 MeV (with the further specifications discussed before) and all particles but nucleons (typically pions) have been emitted or absorbed. At the end of the INC stage a few particles may have been emitted leaving a so-called “compound” nucleus (actually the nucleus is not yet at all in an equilibrated state and the term “compound” is somewhat incorrect, but it is in common use). The nuclear configuration is characterized by:

• the total number of protons, Z_pre⁰, and neutrons, N_pre⁰, of the compound nucleus;

• the number of particle type excitons (nucleons excited above the Fermi level), n_p (n_p = n_pro + n_neu);

• the number of hole-like excitons (holes created in the Fermi sea by the INC interactions), n_h;

• the compound nucleus excitation energy U ($U=\sqrt{s}-M_{N_{\mathrm{pre}}^0,Z_{\mathrm{pre}}^0}$, where $\sqrt{s}$ is the centre of mass energy of the system);

• the compound nucleus momentum components, p_{i comp}.

All the above quantities can be derived by properly counting what occurred during the INC stage and they represent the input configuration for the pre-equilibrium stage.

Already in 1966 Griffin [24] described the spectra following nucleon-induced reactions in terms of a “preequilibrium” model, that is, a transition between the first step of the reaction and the final thermalization.

Since then, many models have been developed (see [25] for an exhaustive review). The two leading approaches (with many different implementations) are the quantum-mechanical multi-step model [67], which has a solid theoretical background, but is complex and entails some difficulties for the description of multiple nucleon emissions, and the exciton model [24, 25, 68], which relies on statistical assumptions.

The exciton formalism employed in PEANUT follows that of M. Blann and coworkers [69–72], called the Geometry Dependent Hybrid Model (GDH).

The pre-equilibrium process in the exciton model is described as a chain of steps, each step corresponding to a particular configuration of “excitons”, where an exciton can be either a particle above the Fermi surface or a hole below the Fermi surface. The statistical assumption underlying the exciton model is that any partitioning of the excitation energy E among n, n = n_h + n_p, excitons has the same probability of occurring. The nucleus proceeds in the chain by nucleon–nucleon collisions, which increase the number of excitons by two units. The chain stops and equilibrium is reached when either the exciton number n is sufficiently high or the excitation energy is below any emission threshold.

The initial number of excitons depends on the reaction type and the cascade history.

At each step, there is a definite probability P_x, n(ϵ) of emitting a nucleon of type x and energy ϵ in the continuum. This probability depends on the exciton configuration, the inverse cross-section, and the re-interaction probability.

The GDH [70, 72] formulation includes an impact parameter dependence on all the position-dependent quantities entering the model, and it constrains the exciton state density to take into account the variation of Fermi energy with position. The FLUKA approach (see [8] for details) is similar but with a few differences. First of all, geometry-dependent quantities are no longer averaged as a function of the impact parameter, but localized, event-specific, values are used, since the position of the first interaction is known from INC. This locality is progressively released with an increasing number of excitons by averaging over increasingly larger fractions of the nucleus.

For the exciton re-interaction rate, nucleon–nucleon cross-section corrected [73] for Fermi motion and Pauli principle, have been used. Further corrections connected to nucleon correlations and to coherence considerations after scattering events have been introduced when computing the re-interaction rate, for consistency with what is done in the INC. These corrections proved to be very useful and prevented the need for somewhat arbitrary reduction factors of the nucleon–nucleon cross-sections, which were often required in similar models to match the experimental data.

In the exciton model, the angular dependence has to be somehow added, since it is not intrinsic to formulation. A non-isotropic angular distribution has been implemented in PEANUT, following the fast particle approximation [74], as implemented by Akkermans et al. [75].

At the end of the pre-equilibrium stage, a true compound nucleus is left with Z_res and N_res, moving with ${\overrightarrow{p}}_{\mathrm{res}}$, and with excitation energy U. The evaporation/fission/fragmentation stage is then simulated starting from this configuration.

3.2.8. Evaporation/Fragmentation/Fission

Evaporation in FLUKA is based on the standard Weisskopf–Ewing formalism [76], which is an application of the detailed balance principle. A description of the basic formalism can be found in [8, 77]. Over the years, the algorithm has been vastly improved with the addition of

• Improved state density $\rho \propto exp(2\sqrt{\mathit{aU}})/U^{5/4}$ where U is the excitation energy and a the level density constant

• Analytic integration of the emission widths

• Competition with γ emission

• Sub-barrier emission, obtained by modification of the inverse cross-section:

  $$\begin{array}{c}\hfill {\sigma }_{\mathrm{inv}}^x={(R+\mathrm{\lambda ̅})}^2\frac{ħ{\omega }_x}{2E}ln(1+exp(\frac{2\pi (E-V_{\mathrm{C}})}{ħ{\omega }_x}))\end{array}$$

  where V_C is the Coulomb barrier, $\overline{\lambda }$ the particle-reduced De Broglie wavelength, and ω_x a penetrability factor that ranges between 2 and 6 MeV, depending on the particle type

• Emission of heavy fragments up to A ≤ 24

For low-mass residual nuclei (A ≤ 16) the evaporation is substituted by a phase space Fermi Break-up model[8, 78, 79].

Fission is also in competition with evaporation. The algorithm in FLUKA [8, 77] follows the approach of Atchison [80], with several recent improvements that will be described in Section 4.13.

The evaporation/fission stage is of utmost importance for the correct calculation of residual nuclei. Examples of the results obtained with FLUKA are shown in Figures 10 and 11.

3.2.9. Residual nucleus de-excitation

The evaporation stage ends when the nuclear excitation energy becomes lower than all separation energies for nucleons and fragments. This residual excitation energy is then dissipated through a cascade of consecutive photon emissions until the ground state is reached. The gamma de-excitation algorithm in FLUKA proceeds through statistical emission, rotational emission, and transitions among experimentally known levels, as described in [81].

More recently, the full set of data on nuclear levels and gammas has been retrieved from the RIPL [82] database and is exploited in all stages of the reaction chain.

Examples of applications to prompt photon emission as a monitoring tool for hadron therapy are published in [6, 83].

3.3. Nucleus–nucleus interactions

Reactions initiated by ions are dealt with by different event generators, depending on the projectile energy. The code switches gradually between event generators at the threshold energies.

At the highest energies (up to 10¹¹ GeV/n) and down to approximately 12.5 GeV/n, FLUKA calls DPMJET-III [84–87]. From Fluka2021 on, the latest version, 19.3, of DPMJET-III [88] is available through an improved interface.

For ions in the few GeV/n energy range and down to approximately 0.1 GeV/n, FLUKA uses an interface to a modified version of rQMD–2.4 [89, 90]. rQMD is a relativistic quantum molecular dynamics model that can also be run in intra-nuclear cascade mode. Examples of FLUKA results compared with experimental data when running with the modified RQMD–2.4 model can be found in [86, 91].

The Boltzmann Master Equation (BME [92]) model has been implemented into FLUKA to deal with the lowest energies, below about 150 MeV/n. BME describes the thermalization of a composite nucleus, created in the complete or incomplete fusion of two ions. The FLUKA implementation [93] samples from the results of the numerical integration of the BME. Complete fusion covers the lowest impact parameter interval. For more peripheral collisions both a three-body description of the reaction and a single nucleon break-up/transfer are modelled. For all event generators, the final de-excitation of the remaining equilibrated nucleus is handled by the FLUKA evaporation/fission/fragmentation module. Moreover, in recent years both rQMD and BME have been interfaced with the PEANUT pre-equilibrium module [83].

3.4. Transport and interaction of low-energy neutrons

In most FLUKA applications, neutrons are a very important radiation component. For neutron energies above 20 MeV, FLUKA handles neutron interactions, as for all other hadronic interactions, with microscopic models. Below this threshold, the complexity of nuclear physics requires specialised algorithms. The standard method for low-energy neutron transport in FLUKA is a multi-group algorithm based on publicly available evaluated nuclear data files. Multi-group transport is a perfect tool for many applications because it is fast and reliable. Comparisons of the FLUKA multi-group transport with experimental data have been published already in ’90 [94, 95]. However, it requires an estimated spectrum for the group averaging of the cross-sections and, by definition, does not allow a direct link between the exact neutron energy and the corresponding cross-section/reaction products. Furthermore, due to the inclusive distribution provided by the data sets, energy conservation is only ensured on average and not at each interaction. A major step forward has been realized in FLUKA and made available to the public with the FLUKA 2021 version: the point-wise fully correlated neutron transport, described in Section 5.2.

3.4.1. Group-wise cross-sections

The multi-group technique, which is widely used in low-energy neutron transport codes, splits the energy range of interest into a given number of intervals (“energy groups”). Elastic and inelastic reactions are not simulated as exclusive processes, but by transfer probabilities from group to group, forming a so-called “down-scattering matrix” (up-scattering is also possible in the thermal energy range). In the FLUKA neutron cross-section library, the energy range up to 20 MeV is divided into 260 energy groups of approximately equal logarithmic width. 31 of these groups are thermal. The angular probabilities for neutron scattering are obtained by discretising a P5 Legendre polynomial expansion of the actual scattering distribution, preserving its first 6 moments. A multi-group scheme is also used to treat gamma production by low-energy neutrons (but not gamma transport). There are 42 gamma energy groups in the FLUKA library, covering the range from 1 keV up to 50 MeV. The actual energy of the produced photon is randomly sampled in the energy interval corresponding to its group, except for a few important gamma lines such as the 2.2 MeV capture transition of deuterium and the 478 keV photon from the ¹⁰B(n,α) reaction. These photons are transported in the same way as all other photons in FLUKA, using continuous cross-sections and an explicit and detailed description of all their interactions with matter. This allows the production of electrons, positrons, and even secondary particles from photo-nuclear reactions. The recoil protons from the interactions on the hydrogen and the protons from the ¹⁴N(n,p) reaction as well as the α particles from ¹⁰B(n,α), are explicitly produced and transported. No other charged secondary particles, including fission fragments, are transported but their kinetic energy is deposited at the point of interaction (kerma approximation). For most materials, FLUKA libraries contain group-dependent information on residual nuclei produced by low-energy neutron interactions, which can be used for activation studies. The standard FLUKA neutron cross-section library includes more than 250 different materials (natural elements or single isotopes), selected for their interest in physics, dosimetry, and accelerator technology. The preparation of the library involves the use of the NJOY code [96, 97] and several ad hoc programs written to tailor the output to the particular FLUKA structure [98]. Hydrogen and carbon cross-sections, which are of particular importance for the neutron slowdown, are also available for different types of molecular binding: free gas, H₂O, CH₂, and ice for hydrogen, free gas, and reactor graphite for carbon. For some important elements (Al, Fe, Cu, Au, W, Pb, Bi) cross-section sets with different degrees of self-shielding have been included in the FLUKA libraries. Some cross-sections are available in the library at two or three different temperatures, mainly for simulations of calorimeters with cryogenic media. Doppler broadening is taken into account.

Point transport was and is available as a user option for a few isotopes because of their relevance: ¹H, ²H, ³He, ⁴He and ¹²C. For the reactions ¹⁰B(n,α), ⁷Li and ¹⁴N(n,p)¹⁴C, the α and the proton respectively are explicitly generated.

3.5. EMF

Electro-Magnetic Fluka (EMF) deals with the transport, interaction and scattering of electrons from 1 keV to 10⁸ TeV and photons from 100 eV to 10⁸ TeV. Photon atomic interaction cross-section are taken from evaluated databases (currently EPICS2017 [99], see Sect. 6), or from specially developed FLUKA algorithms as it is the case for Compton [100], bremsstrahlung, and for the Landau–Pomeranchuk–Migdal effect for both bremsstrahlung and pair production [4].

The treatment of the photoelectric effect takes into account the angular distribution of the photo-electrons in a fully relativistic manner [101]. In addition, fluorescence and Auger electron emission are simulated for six K and L single sub-shells [102, 103]. Compton scattering is described with full consideration of orbital electron motion and binding, with a shell-by-shell, element-by-element description of atomic orbitals (see [100] for details).

Photon polarization is considered when modelling photoelectric, Compton and coherent (Rayleigh) interactions. Accurate sampling of the electron–positron energy and angular distribution is accounted for when modelling pair production.

The code inherently distinguishes between electrons and positrons. As a result, it is able to reproduce specific differences in their interactions with matter. These include stopping power, pair production and bremsstrahlung. In FLUKA, the latter is treated by combining low-energy results from Berger and Seltzer [104] for the differential cross-section data with the analytical high-energy DBMO theory including Coulomb and Elwert corrections at energies above 50 MeV. The treatment takes into account the Landau–Pomeranchuk–Migdal (LPM) [53, 54] and Ter–Mikaelyan [105] (soft photon suppression polarization) effects.

Photon–nuclear interactions, including giant resonance, quasi-deuteron, and resonance interactions as well as virtual photon–nuclear interactions originating from electrons, positrons, high-energy hadrons, and muons, are implemented from the threshold up to approximately 100 EeV [102]. More details on the implementation of the photon–nuclear interactions are given in the following sections. These interactions have recently been thoroughly revised and improved.

Optical photons can be generated and transported. Optical properties of materials have to be provided by the user, as well as emission spectra and intensities of scintillation light. Cherenkov photons are generated according to the standard emission law. Transport includes coherent scattering, absorption, reflection, and refraction on boundaries. Transport time is calculated and recorded. Examples of FLUKA performances in optical photon applications can be found in [106, 107].

Electromagnetic products of nuclear interactions are also explicitly treated by EMF since it is fully coupled to the hadronic interaction models.

3.6. Charged particle transport

The transport of charged particles includes continuous energy loss, energy loss straggling, δ-ray production and multiple Coulomb scattering [4, 83, 108], fully integrated with transport in a magnetic field when requested.

3.6.1. Ionization losses

The most important atomic processes that charged particles undergo during their passage through a medium are related to Coulomb scattering, both with the nuclei and with the electrons. Interactions with the latter give rise to electronic stopping power, which is the dominant process for particle energy loss down to very low energies, resulting in the characteristic shape of the depth-dose profile for charged particles heavier than electrons with the Bragg peak at the end. The energy loss of charged particles is usually expressed as the mean energy loss per unit path length and the corresponding fluctuations around the mean. In FLUKA the electronic stopping power for heavy particles is calculated based on an improved Bethe-Bloch formalism including higher-order corrections [109–113], starting from 1 keV energies up to the maximum energies supported by FLUKA. Approximate transport of particles below the energy cut-off is also performed according to their residual range.

A description of the implementation in FLUKA of the electronic stopping power and associated fluctuations for particles heavier than electrons can be found in [83, 114]. Here a few reminders are given, together with a focus on the practical implementation of the Mott correction, which is important for medium-heavy nuclei.

The formula for the average, unrestricted, energy loss of particles much heavier than electrons and with charge z, can be expressed by:

$$\begin{array}{cc}\hfill {\left(\frac{\mathrm{d}E}{\mathrm{d}x}\right)}_0& =\frac{2\pi n_{\mathrm{e}}{r}_{\mathrm{e}}^2{m}_{\mathrm{e}}{c}^2{z}_{\mathrm{eff}}^2}{{\beta }^2}[ln\left(\frac{2m_{\mathrm{e}}{c}^2{\beta }^2{T}_{\max }}{I^2(1-{\beta }^2)}\right)-2{\beta }^2\hfill \\ \hfill & +2z{L}_1\left(\beta \right)+2z^2{L}_2\left(\beta \right)+M_C(z,\beta )\hfill \\ \hfill & -2\frac{C\left(\beta \right)}{Z}-\delta \left(\beta \right)]\hfill \end{array}$$(16)

for spin 0 particles and by a corresponding expression for spin $\frac{1}{2}$ particles.

β is the projectile velocity relative to the speed of light, n_e is the target material electron density ($n_{\mathrm{e}}=\frac{\rho N_{\mathit{Av}}Z}{A}$ for a given element), I its mean excitation energy, M is the projectile mass, $\gamma =\frac{1}{\sqrt{1-{\beta }^2}}$, T_i = (γ−1)Mc², p_i = βγMc², are the projectile kinetic energy and momentum, and T_max is the maximum energy transfer to a stationary electron, which is dictated by kinematics and given by ($\sqrt{s}$ is the centre-of-mass energy):

$$\begin{array}{c}\hfill T_{\max }=\frac{2m_{\mathrm{e}}{c}^4{p}_i^2}{s}=\frac{2m_{\mathrm{e}}{c}^2{\beta }^2{\gamma }^2}{1+2\gamma \frac{m_{\mathrm{e}}}{M}+{\left(\frac{m_{\mathrm{e}}}{M}\right)}^2}.\end{array}$$(17)

The terms δ, C/Z, L₁, L₂, and M_C are all corrections to the Bethe-Bloch formalism. The “density correction” δ, and the shell correction C, which are important at high and low energy respectively, will not be discussed further. The reader can refer to the extensive literature about these corrections, and to references [83, 114] for details about the FLUKA implementation.

z_eff is the projectile “effective charge” which takes into account the partial neutralization of the projectile charge z when its velocity is not much larger than those of the atomic electrons, the FLUKA implementation is briefly discussed in [4, 83], as well as the implementation of the z³, Barkas [115], and z⁴, Bloch [112], corrections (indicated by L₁ and L₂).

The expression for the restricted energy loss of particles much heavier than electrons and charge z, with energy transfers to atomic electrons restricted at T_δ, is given by:

$$\begin{array}{cc}\hfill {\left(\frac{\mathrm{d}E}{\mathrm{d}x}\right)}_{0T_{\delta }}& =\frac{2\pi n_{\mathrm{e}}{r}_{\mathrm{e}}^2{m}_{\mathrm{e}}{c}^2{z}_{\mathrm{eff}}^2}{{\beta }^2}[ln\left(\frac{2m_{\mathrm{e}}{c}^2{\beta }^2{T}_{\delta }}{I^2(1-{\beta }^2)}\right)-{\beta }^2(1+\frac{T_{\delta }}{T_{\max }})\hfill \\ \hfill & +2z{L}_1\left(\beta \right)+2z^2{L}_2\left(\beta \right)+M_{\mathrm{C}}{(z,\beta )}_{T_{\delta }}-2\frac{C\left(\beta \right)}{Z}-\delta \left(\beta \right)]\hfill \end{array}$$(18)

for spin 0 particles and by a corresponding expression for spin 1/2 particles.

The Mott correction is applied as a parameterized factor in FLUKA based on offline calculations [114], to yield the appropriate stopping power factor based on the relationship between Mott and Bethe differential cross-section ($\frac{\mathrm{d}{\sigma }_M}{\mathrm{d}ϵ}$ and $\frac{\mathrm{d}{\sigma }_B}{\mathrm{d}ϵ}$, respectively)

$$\begin{array}{c}\hfill M_C(z,\beta )\propto \int (\frac{\mathrm{d}{\sigma }_M}{\mathrm{d}ϵ}-\frac{\mathrm{d}{\sigma }_B}{\mathrm{d}ϵ})ϵ\mathrm{d}ϵ\end{array}$$(19)

where ϵ is the transferred energy. Secondary electron production spectra and energy loss fluctuations are corrected as well. Overall, straggling in FLUKA simulations, using the Mott cross-section, was verified to be within the intervals of measured variances for various combinations of relativistic ions and targets as measured in [116].

The contributions of radiative losses for positrons and electrons due to photons below the threshold for explicit production are computed consistently with the bremsstrahlung data, accounting for the Ter-Mikaelyan and LPM corrections when required. This contribution is added to the continuous losses due to ionization.

3.6.2. Ionization fluctuations

Ionization fluctuations are reproduced with a FLUKA-specific algorithm, which represents an alternative implementation with respect to the Landau-Vavilov distribution [117]. The FLUKA algorithm exploits statistical properties of cumulants of distributions, particularly the cumulants of the distribution of Poisson distributed variables, applied to elementary charged hadrons, ions and leptons [4].

For a given particle and known energy, step length, and secondary electrons, δ, explicit production threshold combination, the cumulants of the energy loss distribution can be computed starting from those of the δ production cross-section. The problem is then reduced to sampling from a distribution with given cumulants/moments. In FLUKA, this is solved via the expansions given in [118], for transforming a Gaussian random variate into a variate of given cumulants, up to the 6^th order. This model can be applied for discrete losses below a predetermined δ-ray production threshold, taking into account the latter production cross-sections as well as the effect of Mott corrections, matching the average restricted stopping power. Distant collision fluctuations [119] can be included in the same formalism [4]

3.6.3. Multiple Coulomb Scattering

The transport of charged particles in FLUKA features an original implementation of multiple scattering, based on the Molière theory of Multiple Coulomb Scattering [108]. It is characterized by a “condensed history” approach which optimizes computation time since it does not sample all interactions individually and also accounts for correlations between the final step angular distribution, lateral deflections (average/second moment), path-length corrections and its variance. This methodology greatly enhanced the accuracy of energy deposition/track-length scoring over a single step, further complementing the automatic energy deposition apportioning algorithm of FLUKA.

A rejection technique is used to apply nuclear form factors and spin-relativistic corrections at the 1^st or 2^nd Born approximation level [108]. Additionally, it includes the Fano correction for charged hadrons and muon multiple scattering [122].

Optionally, a Single Scattering method can be activated albeit with a CPU penalty [3]. It consists of an algorithm based on the Mott formula with a screening factor consistent with Molière’s theory implementation in FLUKA. It is also capable of reproducing electron back-scattering and energy deposition in very thin material layers where Molière theory does not apply.

3.7. Neutrino interactions

FLUKA has its own neutrino interaction generator, that embeds the basic neutrino–nucleon interaction in the PEANUT nuclear environment. Quasi-elastic neutrino interactions are sampled according to the formalism in [133], accounting for the mass of the produced lepton. Validation of the algorithm and of the Fermi motion as implemented in PEANUT is described in an experimental paper about neutrino interactions in liquid argon [134]. At higher energies, the resonant and deep inelastic scattering (DIS) channels open. The production of Δ resonances is modelled following Rein–Sehgal [135]. No other resonance is considered. Care is taken to avoid double-counting with the non-resonant DIS channel, with a linear decrease (increase) of the resonance (DIS) probability as a function of the mass of the intermediate hadronic system.

The (anti)neutrino-nucleon DIS generator, called NUNDIS, is described in references [136, 137]. The NUNDIS package treats neutral current (NC) and charged current (CC) DIS interactions for neutrinos and antineutrinos incident on protons and neutrons, based on standard Parton Distribution Function (PDF) sets, in the energy range from threshold to at least 10 TeV. It samples the energy of the outgoing lepton, its angle and (for CC) the polarization. It also returns a set of variables to the PEANUT model, including the squared four-momentum transfer, Q², the squared mass of the hadronic system, W², and the flavour of the partons entering and leaving the interaction. Based on these variables, the FLUKA hadronization then proceeds to the creation of the hadrons, and then the standard PEANUT mechanisms take care of their propagation and, possibly, of the final state interactions in the target nucleus.

At present, NUNDIS adopts the GRV98 [138] PDF set in the NLO (DIS) version. The Extrapolation from Q² = 1 GeV² to Q² = 0 has recently been improved by adopting the formulation of Bertini et al. [139].

NUNDIS predictions in terms of total cross-section agree very well with available experimental data, even at the highest energies where there is no longer a linear dependence with the neutrino energy, as shown in Figure 12.

The output from the sampling of DIS events, at parton level, is passed to the hadronization package of FLUKA, to be treated in the same way as for hadron–hadron interactions. Thus, neutrino events profit from the original FLUKA treatment of low-mass chains, which allows for accurately reproducing low-DIS events (see [6]), of particular interest for the next generation neutrino oscillation experiments.

Initial state effects (essentially Fermi motion) and final state effects are simulated in the PEANUT framework.

NUNDIS was able to accurately reproduce event rates and topologies in the ICARUS experiment at Gran Sasso as in [140]. The unique capability of FLUKA to describe gamma de-excitation after neutrino interaction is demonstrated in the analysis of low-energy events in the liquid argon detector ArgoNeuT [141].

3.8. Online evolution of activity and dose

Both PEANUT and the low-energy neutron algorithm provide information about residual nuclei. FLUKA allows the user to calculate online the time evolution of radioactivity and to propagate radiation products at selected times after the irradiation. “Online” refers here to the same run used for the prompt simulation. The advantage with respect to commonly used two-step methods, where the residual nuclei distribution is stored and subsequently transmitted to an evolution code, is both the easiness of use and the fact the exact spatial distribution of the residual nuclei is retained.

The time evolution of radioactive isotopes is calculated analytically with an algebraic solution of the Bateman [142] equations that incorporate decay and build-up during and after irradiation [143]. Up to seven different decay branching ratios for each isotope are included in the decay tree. 51 different decay modes are supported, including all those with delayed particles, like β⁻n. Known isomeric states are included in the evolution.

The generation and transport of decay radiation is possible. A dedicated database of decay emissions has been created using mostly information from NNDC [144], sometimes supplemented by additional data, and checked for consistency. Emitted and transported products include γ, α, β^±, X-ray conversion, and Auger electrons. Decay products are not correlated, meaning that, for instance, an α transition to level i in the daughter nucleus could be followed by γ rays competing to level j. However, in all cases, the primary transition is correctly sampled and never doubly counted. Particular care is taken in the simulation of β spectra, with the inclusion of Coulomb, recoil, and screening corrections.

The user can specify arbitrarily complex irradiation profiles with multiple intensities and time intervals, as well as multiple cooling times at which the various quantities are to be evaluated. The user can even choose to assign different material compositions to selected regions for immediate and delayed radiation transport.

4. Recent developments in the hadronic interaction models

4.1. Colour transparency

The observed hadron–nucleon cross-section, σ_hN, is just the average of the cross-sections corresponding to all possible (quark) configurations of the incoming hadron.

Considering the incoming hadron as a colour dipole (q-qbar for mesons, q-qq for baryons), and with some hypotheses on the distribution of spatial configurations, P(ρ²), with the constraint ∫dρ²P(ρ²) = σ_hN, a fluctuating σ can be used within the Glauber formalism.

Recalling equation (5) in Section 3.2.1, in the Glauber formalism, the absorption (particle production) cross-section can be written as ${\sigma }_{hAabs}(s)\equiv \int {\mathrm{d}}^2\overrightarrow{b}{\mu }_{hAabs}(\overrightarrow{b},s)$.

Let us introduce the thickness function for non-elastic reactions, T_rj, where ρ_j is the density distribution for the j^th target nucleon, $T_{\mathit{rj}}(\overrightarrow{b},s)$ is the amount of nuclear matter seen by the incident hadron traveling along the impact parameter $\overrightarrow{b}$, when folded with its profile function, ${\kappa }_{\mathit{hN}}(\overrightarrow{b},s)$.

$$\begin{array}{cc}\hfill T_{\mathit{rj}}(b)& ={\int }_{-\infty }^{+\infty }\mathrm{d}z\int {\mathrm{d}}^2\mathrm{s⃗}{\rho }_j(\sqrt{z^2+s^2}){\kappa }_{hNj}(\sqrt{b^2+s^2-2bscos\varphi })\hfill \\ \hfill & ={\int }_{-\infty }^{+\infty }\mathrm{d}z\int {\mathrm{d}}^2{\mathrm{s⃗}}^{\prime }{\rho }_j(\sqrt{z^2+b^2+s^{\prime 2}-2b{s}^{\prime }cos\varphi }){\kappa }_{hNj}(s^{\prime }).\hfill \end{array}$$

Assuming that ρ_j does not depend on the nucleon index and that σ_{hn r} = σ_{hp r} ≡ σ_{hN r}, the μ factor defined in equation (5) becomes:

$$\begin{array}{cc}\hfill {\mu }_{hAabs}(b)& =1-{[1-{\sigma }_{hNr}{T}_r(b)]}^A\hfill \\ \hfill & =\sum _{\nu =1}^A\left(\begin{array}{c}A\\ \nu \end{array}\right){\left[{\sigma }_{hNr}{T}_r(b)\right]}^{\nu }{[1-{\sigma }_{\mathit{hN}}{T}_r(b)]}^{A-\nu }.\hfill \end{array}$$

Since σ_hNr T_r(b) is the probability of getting one specific nucleon hit and there are A possible trials, this is exactly the binomial distribution for getting ν out of A trials.

The Glauber formalism should be further developed including the so-called inelastic screening correction [145, 146]. This correction arises because the Glauber expression for the scattering amplitude contains only contributions from elastic hadron–nucleon scattering. There are additional contributions in second and higher-order diagrams from intermediate states in which the incoming hadron is diffractively excited on one target nucleon, and reverts to its ground state in a later interaction.

Introducing cross-section fluctuations into the Glauber formalism corresponds to substituting:

$$\begin{array}{cc}\hfill \langle {\mu }_{hAabs}(b)\rangle & \approx 1-{[1-⟨{\sigma }_{hNr}⟩T_r(b)]}^A\hfill \\ & =\sum _{\nu =1}^A\left(\begin{array}{c}A\\ \nu \end{array}\right){\left[⟨{\sigma }_{hNr}⟩T_r(b)\right]}^{\nu }\hfill \\ & {[1-⟨{\sigma }_{\mathit{hNr}}⟩T_r(b)]}^{A-\nu }\hfill \end{array}$$

with:

$$\begin{array}{cc}\hfill \langle {\mu }_{hAabs}(b)\rangle & =\langle 1-{[1-{\sigma }_{hNr}{T}_r(b)]}^A\rangle \hfill \\ & =\langle \sum _{\nu =1}^A\left(\begin{array}{c}A\\ \nu \end{array}\right){\left[{\sigma }_{hNr}{T}_r(b)\right]}^{\nu }\hfill \\ & {[1-{\sigma }_{\mathit{hNr}}{T}_r(b)]}^{A-\nu }\rangle \hfill \end{array}$$

where now the proper expression averaging is used instead of using the average (observed) cross-section, ⟨σ_{hN r}⟩, into the μ_abs expression.

This approach is often referred to in the literature as colour transparency. It automatically takes into account the first-order inelastic screening correction to the Glauber formalism, thus allowing the calculation of reliable cross-sections from 1 GeV up to cosmic ray energies. The practical implementation in PEANUT is rather complex, since the configuration changes on a timescale of approximately 1 fm/s proper time, i.e. up to 10–20 GeV, for heavy nuclei, σ can fluctuate within the width of the target nucleus.

As well as being essential for the accurate evaluation of hadron–nucleus cross-sections, colour transparency also affects the particle multiplicity distributions. In particular, the spread of the multiplicity distributions becomes significantly larger, and longer tails appear at large multiplicities, especially for heavy targets.

Two examples of the effect of the cross-section fluctuation implementation are shown in Figures 13 and 14. In the former, the measured negative particle multiplicity distributions for 250 GeV K⁺ on aluminium and gold [40] are compared with PEANUT predictions with and without the colour transparency algorithm enabled. The latter shows the calculated fast and negative particle multiplicity distributions for 7 TeV protons on tungsten, again with and without the colour transparency algorithm.

4.2. Improved hadron–nucleus cross-sections

Many problems related to cosmic ray propagation and interactions have been studied with FLUKA, among them [147–151].

To handle very high-energy cosmic rays, the code must include accurate cross-sections for nuclear interactions up to energies of the order of 10²⁰ eV.

Significant extensions to the Gribov-Glauber high-energy collision model in PEANUT now allow first-principles calculation of quasi-elastic (QE, see Eq. (4)) and particle-production hadron–nucleus cross-sections, once the hadron–nucleon scattering amplitudes are known. These extensions include the colour transparency algorithm described in the previous section and other refinements necessary to make the model applicable up to the highest cosmic range energies. FLUKA (PEANUT) predictions for proton–carbon and proton–lead cross-sections are compared with experimental data in Figures 15 and 16. Calculated particle-production cross-sections for proton-air are shown together with experimental data from CR experiments up to 10¹⁹ eV in Figure 17. The agreement with the most recent measurements is particularly good.

From Fluka2023 onwards, the old high-energy parameterizations developed by the Leipzig group in the 1980s have been replaced by the PEANUT-calculated QE and particle production cross-sections for all hadron–nucleus combinations.

4.3. New hadronization model

Several experiments, including some very recent ones, have measured surprisingly large yields of ρ⁰ mesons in both pion–nucleon and pion–nucleus interactions, especially in the forward region (see [153] for a review and discussion).

In addition, the ratio of the ρ⁰/ω/π⁰ mesons that share the same quark composition is measured to be strongly dependent on the velocity/Feynman x_F. Such a result is in contradiction to a critical assumption of all hadronization models that have been used so far in hadronic interaction codes, the so-called chain universality. This result is not only surprising, but also important for very energetic cosmic rays, since it implies a deep reassessment of hadronization assumptions for hadron projectiles. Indeed, a leading ρ⁰ instead of π⁰ can slow down the EM fraction increase in air showers (the ρ⁰ decays into two charged pions instead of two photons) and thus increase the hadronic and muon content of the resulting showers. This process is also important for calorimeters because it can affect the e/h ratio. This is a critical parameter in determining the compensation characteristics of a detector and ultimately its energy resolution.

The observation that, contrary to expectations, the universality of jet hadronization does not hold in soft physics is also evident at much higher energies, such as those probed by the LHC experiments.

Most of the assumptions about hadronization have their origin in the observation of hadron production in e⁺e⁻ colliders. However, a fundamental difference in hadronic collisions is that the projectile and the target (anti)quarks are valence quarks, rather than sea quarks, and as such they carry quantum numbers that are different from those in the vacuum, in particular the isospin. Therefore, it is not unreasonable to assume that the hadronization of chain ends containing valence quarks may be somewhat different from that of pure sea quarks, especially for hadrons that differ in their isospin, such as ρ⁰ and ω.

Figures 18 and 19 show the effect of the modifications for 250 GeV/c positive pions on hydrogen and 158 GeV/c negative pions on carbon. The right-hand plots show the “old” hadronization results, where the “standard” assumption of fixed ratios between ρ⁰/ω/π⁰ mesons is used. In the left-hand plots, these ratios are assumed to be different when valence quarks are involved, while still preserving exactly all isospin symmetries. The improvement in the ρ⁰ and ω production compared to experimental data [152, 153] is striking.

This improvement can still be seen in the inclusive positive particle, negative particle, and π⁰ distributions after resonance decays as shown in Figure 20.

Finally, an example of the present performances of the FLUKA hadronic models at higher energies is presented in Figures 21 and 22 for proton-proton collisions at $\sqrt{s}=31\mathrm{GeV}(E_{\mathrm{lab}}\approx 512\mathrm{GeV})$ and 62 GeV (E_lab ≈ 2050 GeV).

4.4. Improvements at intermediate energies

A number of refinements have been introduced into the hadron–nucleon resonance model, which is used to describe the interactions above the pion production threshold and below the 3–5 GeV/c limit. In particular, the width of the different resonances is no longer a fixed value but rather depends on the actual mass and spin of the resonance. This leads, especially for energies close to the pion production threshold, to significant improvements in the description of pion production in the energy range of the resonance model. It is also crucial for the correct description of photo-nuclear interactions, where the actual shape of the ρ⁰ resonance (a rather broad resonance) plays a key role. An example of these improvements can be seen in Figures 23 and 24, which show the positive pion production for protons on proton and carbon, respectively, at 585 MeV.

4.5. High-energy: DPMJET-III-19.3

The DPMJET event generator comprises three stages. For E_lab >  5 GeV, hadron–nucleon interactions are managed by the PHOJET event generator [165, 166], making DPMJET events and cross-sections identical to PHOJET’s in these instances. The hadronization of colour strings is handled by the Lund String Fragmentation model from PYTHIA 6 [167]. The DPMJET framework employs these two codes for modelling nuclear phenomenology, which shares many common points with PEANUT, such as the Glauber–Gribov approximation or the formation zone intranuclear cascade. DPMJET, along with PEANUT, is one of the few event generators capable of simulating photo-nuclear interactions at high energies [168].

The main differences between PEANUT and DPMJET for hadron–hadron interactions are the models of minijets and multi-parton interactions (MPI), which become important at $\sqrt{s}$ > a few hundred GeV (FLUKA switches to DPMJET at E_lab >  15 TeV for hadron–hadron and hadron–nucleus collisions). While DPMJET explicitly simulates INC, it does so in a simplified way, without the extra ingredients that distinguish PEANUT GINC from classical INC, in particular, multibody processes, refraction and reflection, coherence zone, nucleon anti-correlations, and radial and energy-dependent nuclear densities and potentials. These limitations and the absence of a pre-equilibrium step are not critical at energies above a few tens of GeV due to limited intra-nuclear cascading. More importantly, it also lacks a framework to dissipate the excitation energy of the nuclear remnants. Therefore, the FLUKA evaporation modules handle this final de-excitation stage.

Following the availability of LHC minimum bias measurements, DPMJET-III version 3.0-6, distributed with FLUKA, was significantly updated and re-tuned to incorporate new LHC data [84, 169]. Key changes include:

1. Internal cross-section tables are now computed for each projectile-nucleon combination on demand. The fundamental combinations with available measurements are proton-(anti-)proton, neutron-proton, pion-proton, and kaon-proton, including their charge conjugates. For less abundant projectiles, the cross-section tables are identical to one of these combinations; however, the valence quark content is swapped.

2. The charged particle multiplicity distributions were underestimated at LHC energies and corrected by extending the probability distributions for MPI. This change necessitated adjusting the minijet cross-section’s energy dependence and moderately retuning the model parameters.

3. The parton distribution functions (PDFs), from which sea flavours are sampled in the minijet model and which are used for the calculation of the hard cross-section, have been updated from GRV98 [138] to CT14 LO [170]. Along with the multiplicity correction, this change required a readjustment of the perturbative cutoff parameter $p_{\perp ,\mathrm{cut}}(\sqrt{s})$ and motivated further parameter retuning.

4. Minor changes included code cleanup, the migration from PYTHIA 6.1 to 6.4 [167] for hadronization.

A more detailed view of these changes is given in the following sections.

4.5.1. Updated interaction cross-sections

The hadron interaction cross-sections are calculated within the PHOJET event generator, and DPMJET accesses those for further calculations of the nuclear cross-sections using the Glauber-Gribov approximation. The cross-section fit closely follows the procedure described in Section 5 and Appendix A of reference [165]. DPMJET defines amplitudes for exchanges of a single soft Pomeron, a Reggeon, a hard Pomeron, and for high-mass diffraction, a Triple-Pomeron and a Loop-Pomeron. Low-mass diffraction is implemented using a two-channel Good-Walker paradigm [171] with one generic excited state. These amplitudes are unitarized following the Eikonal approximation. The rise of the cross-sections as a function of s is attributed to the slope of the super-critical Pomeron Δ_P = α(t = 0)−1 >  0, which is universal across different hadron types. Therefore, the high-energy behaviour can be derived from the proton–proton and proton-antiproton data. The fit is shown in Figure 25. Whereas the old DPMJET had a change of slope around Tevatron energies, the new DPMJET demonstrates a clean power-law dependence, constrained by data from the LHC. The hadron–proton and hadron–neutron cross-sections, shown in Figure 26, share most parameters with the proton projectiles. Fitting just the pomeron-particle coupling g_𝒫, 0 and the reggeon-particle coupling g_ℛ, 0 on the non-proton side plus the Reggeon intercept α_ℛ turned out to be sufficient in most cases. We noticed that the elastic slope b_ela at LHC energies cannot be well fitted together irrespective of parameter choices.

4.5.2. Changes to the hard cross-section and the multiplicity distributions

In early comparisons with LHC data, it became clear that events with very high multiplicities were missing in DPMJET. After a lengthy technical investigation, the main reason was identified as an artificial eikonal series truncation in multiple cut tables within the cross-section routines of PHOJET. The extension of the tables resulted in the contrary problem of too high multiplicity of the present cross-section model at $\sqrt{s}$ above 1 TeV and required adjustments of the hard cross-section and some minor parameter tuning. Figure 27 shows the new energy dependence of soft and hard pomeron cuts, as well as the impact of the choice of the p_⊥-cutoff on the charged particle multiplicity distribution, which enters the integration boundaries for the hard cross-section calculation, expressed as the sum of lowest-order perturbative QCD parton-parton cross-sections:

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{QCD}}& =\sum _{i,j,k,l}\frac{1}{1+{\delta }_{\mathit{kl}}}\int \mathrm{d}{x}_1\mathrm{d}{x}_2\underset{Q_{\min }\propto p_{\perp }^{\mathrm{cutoff}}}{\int }\mathrm{d}{Q}^2\hfill \\ \hfill & \times f_i(x_1,Q^2)f_j(x_2,Q^2)\frac{\mathrm{d}{\sigma }_{i,j\to k,l}}{\mathrm{d}{Q}^2},\hfill \end{array}$$(20)

where x₁, x₂ are the longitudinal momentum fractions of the incoming partons i, j, Q the virtuality or the momentum transfer of the process, σ_{i, j → k, l} the leading-order QCD matrix elements, and f_i, f_j the density (PDFs) of flavour i partons in the incoming particle on side 1 and j for the other particle, respectively.

As they were no longer updated with the latest data, the default GRV98 [138] parton distribution functions were replaced with more recent CT14 LO [170] PDFs. Note that in simulations of minimum-bias events or particle cascades, the exact choice of PDFs has little effect compared to simulations of high-p_⊥ jets, for example. What matters more is consistent small x behaviour, i.e. the parton densities must not become negative or diverge, and good performance at low and high values of Q. Since the hard cross-section is defined by the value of the integral in equation (20) above a specific cutoff, there are no additional free parameters. The differences between different PDFs or cutoff behaviours are, therefore, absorbed by the other free parameters of the model in the cross-section fit discussed in the previous subsection.

To improve the compatibility with LHC data and preserve the convergence of the cross-section model, we have chosen a dipole-model-inspired energy dependence [172] of the cutoff parameter:

$$\begin{array}{c}\hfill p_T^{\mathrm{cutoff}}(s)=2\mathrm{GeV}/\mathrm{c}{\left(\frac{\sqrt{s}+400\mathrm{GeV}}{200\mathrm{GeV}}\right)}^{0.19}.\end{array}$$(21)

Using this setup, the cross-section calculation is numerically stable for $\sqrt{s}\simeq 200\mathrm{TeV}$, enabling the calculation for Ultra-High-Energy Cosmic Ray showers with energies up to 100 EeV. Different combinations of parameters were tested against the majority of minimum-bias data available until 2016, and a less steep increase of the hard cross-section was favored compared to the older version of DPMJET.

A comparison with ATLAS multiplicity distributions measured in pp collisions at 900 GeV and 7 TeV is shown in the upper panels of Figure 28 for different phase-space regions. While no data for energies above 7 TeV have been used for tuning, the good description is retained at 8 and 13 TeV (see lower panels in Fig. 28), particularly for lower and medium event multiplicities.

For proton–nucleus and nucleus–nucleus collisions at LHC energies, the majority of the available data is taken at mid-rapidities, which has a moderate impact on calculations of particle cascades or energy deposition profiles. Since the proton–proton cross-section acts as input to the Glauber calculations, it will result in a small change to the production cross-sections. As shown in the proton-lead pseudorapidity density in Figure 29, the recent DPMJET has a better description of ALICE data.

4.6. Deuteron pre-formation

The coalescence mechanism is very effective in reproducing the emission of light fragments at projectile energies in excess of a few tens of MeV. However, at lower energies, coalescence fails to reproduce the observed (p/n,d) reaction cross-sections, particularly on light targets. This is not surprising since the probability of emitting sequentially a proton and a neutron which could coalesce in a deuteron when close to threshold is very reduced. A very effective mechanism for overcoming this issue is to assume that when a proton–neutron first interaction occurs, there is a definite probability that the two outgoing nucleons continue as a deuteron which will propagate and possibly emerge as such in the following cascade or pre-equilibrium evolution of the reaction. Figure 30 illustrates the improvement due to deuteron pre-formation when applied to the ¹²C(p,x)¹¹C and ¹⁶O(p,x)¹⁵O reactions.

4.7. Advanced Fermi break-up

Statistical evaporation of excited low-mass fragments is unsuitable due to the relatively few, widely spaced levels. Therefore, FLUKA makes use of an alternative de-excitation mechanism for these light (typically A ≤ 16) residual nuclei, the Fermi Break-up model [78, 79]. In the Fermi break-up approach the excited nucleus is supposed to disassemble in one single step into two or more fragments, possibly in excited states, with branching given by spin, Coulomb barrier, and plain phase space considerations. In particular, the probability for breaking-up a nucleus of N neutrons, Z protons, and U excitation energy (total mass M^* = U + M_A, Z) into n fragments (n ≥ 2) of the same total charge and baryon number, is given by:

$$\begin{array}{cc}\hfill W& =\frac{S_n}{G}{\left[\frac{V_{\mathrm{br}}}{{(2\pi ħ)}^3}\right]}^{n-1}{\left(\frac{1}{M^{\ast }}\prod _{i=1}^n{m}_i\right)}^{3/2}\hfill \\ \hfill & \times \frac{{(2\pi )}^{3(n-1)/2}}{\mathrm{\Gamma }\left(\frac{3}{2}(n-1)\right)}{E}_{\mathrm{kin}}^{3n/2-5/2}.\hfill \end{array}$$(22)

The spin multiplicity S_n, and the permutation factor G are given by (n_j is the number of identical particles of jth kind)

$$\begin{array}{c}\hfill S_n=\prod _{i=1}^n(2S_i+1),G=\prod _{j=1}^k{n}_j!\end{array}$$(23)

and E_kin is the total kinetic energy of all fragments at the moment of break-up given by

$$\begin{array}{c}\hfill E_{\mathrm{kin}}=M^{\ast }-\sum _{i=1}^n{m}_i-E_{\mathrm{Coul}}\end{array}$$(24)

where E_Coul is the Coulomb energy in case there is more than one charged particle. V_br is a volume of the order of the initial residual nucleus volume. Therefore, the final state can be conveniently selected by means of a MC procedure, by evaluating such an expression for all combinations of fragments energetically allowed and making a random selection. However, equation (22) implicitly assumes that the fragment emission occurs in L = 0 and neglects consideration of the initial spin and parity state, J^π, of the excited nucleus. If the initial J^π is known, now FLUKA applies suitable modifications to equation (22) in order to account for the spin and parity, in particular:

• The minimum orbital momentum L_min compatible with J^π and the spins and parities of the emitted particles must be computed.

• The spin factor S_n must be restricted to the spin projections compatible with an emission with L_min.

• In case L_min >  0, a suitable centrifugal barrier must be added to E_Coul.

Figure 31, right, shows an example of the improvement arising from the inclusion of spin/parity considerations. In this figure the calculated excitation curves for ¹²C(γ, n)¹¹C and ¹⁶O(γ, n)¹⁵O, reactions for which J^π = 0⁺ in the GDR energy range, are compared with experimental data before and after the application of the Fermi break-up enhancements. After their application, the excitation curve for ¹²C(γ, n)¹¹C is in very good agreement with the experimental data. The one for ¹⁶O(γ, n)¹⁵O is significantly improved, particularly in the energy range of the peak of the GDR, while there is still some deficit at higher energies. These reactions are of special relevance for underground experiments, particularly those using liquid scintillators, where virtual or real photons produced by high-energy muons penetrating through the underground rock can produce ¹¹C and neutron background.

4.8. Nuclear coherent elastic scattering

A simple but very powerful model for hadron–nucleus coherent elastic scattering was developed in the 90’s by Ferrari and Sala [178]. The model has been recently extended, and it is now used in FLUKA for nucleon-nucleus coherent elastic scattering below 200 MeV (for a few selected nuclei up to 1 GeV), and for whichever hadron–nucleus combination at energies in excess of 1 GeV.

The basic formula is presented in the following with a few examples of its application.

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\sigma }{\mathrm{d}t}& =\pi [(1-p_{J0})R_1^4{\left|\frac{J_1(q{R}_1)}{q{R}_1}\right|}^2{\mathrm{e}}^{-{\mathrm{\Lambda }}_1^2{q}^2}\hfill \\ \hfill & +\frac{p_{J0}}{2}{R}_0^4{J}_0^2(q{R}_0){\mathrm{e}}^{-{\mathrm{\Lambda }}_0^2{q}^2}]\hfill \end{array}$$(25)

where q is the momentum transfer. This formula is inspired by the grey disk model. Indeed for a black disk of radius R₁ the J₁ term alone would describe exactly nuclear coherent scattering. The presence of a diffuse edge is responsible for the introduction of a J₁ derivative term, the J₀ one. Expressing R_1/0 = r_1/0 A^1/3, where A is the mass number of the target nucleus, the model depends on 5 parameters, r₁, r₀, p_J0, Λ₁, Λ₀. r₀, r₁ are expected to be of the order of 1–1.2 fm and possibly very close one to the other. A careful parameterization based on available experimental data and DWBA calculations [179, 180] indeed matches these expectations and it allows for extrapolating well beyond the limited sets of data used for fixing the parameters. Examples of the model performances for neutrons and protons coherent elastic scattering at energies below 1 GeV can be found in Figures 32 and 33, left, compared with available DWBA calculations and experimental data. Furthermore, above 1 GeV, the parameters scale in a very simple and predictable way with the identity and energy of the hadron projectile, allowing the model to be used for any combination of incident hadron, energy and target. At high energies, experiments cannot distinguish between coherent elastic scattering and quasi-elastic (incoherent elastic) scattering with low momentum transfer. Therefore, comparisons with available experimental data must also take into account the quasi-elastic component within the experimental resolutions. The improved PEANUT Glauber model described in the previous paragraphs allows the calculation of the quasi-elastic cross-section (Eq. (4)) and the associated secondaries. In Figure 33, right, and in Figure 34 three examples of different projectile/energy/target combinations are shown and compared with available experimental data. Further details about this model and its performance will be given in a forthcoming paper.

4.9. Photo-nuclear interaction and photo-fission

The FLUKA photonuclear cross-sections and reaction models have been thoroughly revised in recent years.

At low energies, in the giant dipole resonance (GDR) region, updated cross-sections based on recent evaluations [179, 187] have been substituted in the FLUKA database [188] for several nuclei.

Up to and including Fluka2021.2, real and virtual photonuclear reactions in the Δ range and up to approximately 3 GeV were performed by turning a γ into a π⁰, still conserving energy and momentum, and using PEANUT. This approach had several disadvantages. In fact, photons behave somewhat differently, and in particular their couplings to Δ’s are not subject to isospin rules, unlike pions. Furthermore, the in-medium nuclear effects for γN → Δ, including the 2- and 3-nucleon absorption channels, are not the same as for pions. Last but not least, photons at these energies mostly probe the volume of the nucleus rather than the surface like pions.

Furthermore, γN and γA cross-sections in the Δ region and above were based on old parametrisations. More extensive and accurate data are now available. Above approximately 2–3 GeV, the photo-nuclear reactions assume that the photon behaves like a virtual vector meson, ρ⁰, ω, Φ, according to the Vector Meson Dominance model (VMD). Until Fluka2021.2, the interactions were performed with the old high-energy model of FLUKA, which lacks the sophisticated nuclear physics of PEANUT.

The cross-sections for γp/n and γA in the resonance energy region were re-evaluated as a first step. The current FLUKA cross-sections for γp/n are shown in Figure 35, together with the single pion channel cross-section. The data available at [189] were used for the latter.

The smoothing and eventual disappearance with increasing mass number of the bump in the cross-section of all γN resonances higher than the Δ(1232) is clearly seen in nuclear photoabsorption data on nuclei. The Fermi motion of the nucleon, Pauli blocking and the modification of the Δ properties and higher resonances in the nuclear medium all contribute to this effect. Indeed for all A >  4, the photon–nucleus cross-sections scale linearly with the atomic number of the target nucleus in the energy range below approximately 2 GeV. An example that includes the newly adopted FLUKA “universal curve”, is shown in Figure 37. For very light nuclei, where some residual traces of the higher resonances are still visible, the universal curve does not apply. Therefore, ad hoc cross-sections have been implemented for ²H, ³He and ⁴He and are shown in Figure 36 together with experimental data.

The linear scaling with A shows that the γA cross-sections do not exhibit any shadowing in the considered energy range, and therefore the photon probes the whole nuclear volume as a point-like particle. Now γA interactions for A >  1 and in this energy range are simulated with PEANUT as such (no longer as pseudo-π⁰). The initial target nucleon(s) is (are) chosen from the nuclear volume, and the new parameterizations for the γN cross-sections are used. The in-medium broadening of the intermediate Δ resonance and the resulting 2- and 3-nucleon absorption channels and their dependence on the local nuclear density are described according to the Oset formalism [45] for photons. The MAID07 data [189] for the γN single pion cross-sections and angular distributions (available up to approximately E_γ 1.25 GeV) are now used.

Resonant and non-resonant channels have been added in order to model multi-pion and strangeness production up to approximately 2–3 GeV. The γN cross-section is known to be partly due to resonant intermediate states (γN → Δ^( * )(x) or N^*(x)→ decay products) and to a non-resonant background which rapidly becomes dominant above the Δ(1232) region. In accordance with [190], the intermediate resonance states explicitly considered are: Δ(1232), N^*(1440), N^*(1520), N^*(1535), N^*(1680) and Δ(1950). For the non-resonant part, the following channels are included with an educated guess as to their cross-sections, based on the available experimental data:

$$\begin{array}{cc}\hfill \gamma N& \to \eta (548)N\hfill \\ \hfill \gamma N& \to {\rho }^{0/+/-}(770)N\hfill \\ \hfill \gamma N& \to \omega (782)N\hfill \\ \hfill \gamma N& \to \mathrm{\Phi }(1020)N\hfill \\ \hfill \gamma N& \to K^{+/0}\mathrm{\Lambda }\hfill \\ \hfill \gamma N& \to {\overline{K}}^{0/-}{K}^{0/+}N.\hfill \end{array}$$(26)

The previous channels do not exhaust the non-resonant cross-section. In the new FLUKA approach, the rest is split among the following channels, with weights set in such a way as to reproduce multi-pion production data, in absence of experimental data:

$$\begin{array}{cc}\hfill \gamma N& \to {\pi }^{+/0/-}{N}^{(+/0)\ast }(1440)\hfill \\ \hfill \gamma N& \to \sigma (500){\mathrm{\Delta }}^{+/0}(1232)\hfill \\ \hfill \gamma N& \to {\rho }^{+/0/-}(770){\mathrm{\Delta }}^{++/+/0/-}(1232)\hfill \\ \hfill \gamma N& \to \sigma (500)N^{(+/0)\ast }(1440)\hfill \\ \hfill \gamma N& \to {\rho }^{+/0/-}(770)N^{(+/0)\ast }(1440).\hfill \end{array}$$(27)

At energies progressively higher than 2–3 GeV, the photon is supposed to fluctuate into a vector meson, ρ⁰, ω(782), Φ(1020), with a probability corresponding to the experimental relative coupling. The interaction then proceeds through the following steps:

• A first choice is made between a (point-like) photon interaction (a volume one with no shadowing) and a vector-meson-nucleus one with a hadron-like shadowing;

• If the interaction is not point-like, it is chosen whether the interaction proceeds by diffractive-like coherent pseudo-elastic scattering γ + A → ρ⁰A → ρ⁰A (and equivalent for ω(782) or ϕ(1020)), or by an inelastic vector meson nucleus scattering, the latter including the pseudo-quasi-elastic channel γ + A → ρ⁰A → ρ⁰A^*;

• If the interaction is a point-like one, or a VMD non-elastic one, the vector meson is used as projectile in PEANUT, in the first case with a volume selection of a single target nucleon, in the second case with a hadron-like (Glauber) interaction.

PEANUT can now handle vector meson projectiles and their (re)interaction and decay in the nucleus to support the procedure outlined above.

In addition to the new models for the photon–nuclear interactions above the quasi-deuteron region, the photon–nucleus absorption cross-sections in the Giant Dipole Resonance (GDR) region have been revised and updated for several isotopes using different sources, in particular the evaluated data of references [179, 187].

An example is shown in Figure 39, where the PEANUT computed photo-fission cross-sections are compared with available experimental data. Figure 38 shows the newly adopted photonuclear cross-sections for a Carbon target over the whole energy range, from threshold up to several tens of GeV.

These procedures are also applied to virtual meson interactions, for example when describing Electro-Magnetic-Dissociation of ions [191], or electro- or muon-nuclear reactions.

4.10. ElectroMagnetic dissociation and electroNuclear interaction

ElectroMagnetic dissociation (EMD) in FLUKA has been tested and applied from a few GeV to LHC energies with remarkable results [191]. In the last decade, the model has been further improved, including the electric quadrupole (E2) multipolarity, and nuclear finite size, and higher-order effects. These are important for low-energy ions and at all energies for e^± and μ. At the same time, the FLUKA EMD treatment has been extended to include the photonuclear interactions of the muon below the threshold of γ + N → π + X and the electronuclear interactions. The deuteron Coulomb dissociation is also described using the same formalism. A description, including some preliminary results, of the revised FLUKA EMD approach can be found in [192]. In the following, the main physical assumptions are briefly recalled for completeness.

The EMD cross-section can be customarily expressed as a sum over all multi-polarities of the convolution of the number of emitted (virtual) photons n_i for a given multi-polarity, with the photon-nucleus cross-section σ_i γA for the same multi-polarity. n_i is the equivalent photon number for the ith multi-polarity already integrated over all impact parameters:

$$\begin{array}{c}\hfill {\displaystyle {\sigma }_{\mathrm{EMD}}\approx \sum _i{\int }_{E_{\min }}^{E_{\max }}{\sigma }_{i\gamma A}\left(E_{\gamma }\right)n_i\left(E_{\gamma }\right)\frac{\mathrm{d}{E}_{\gamma }}{E_{\gamma }}}\end{array}$$(28)

and where:

$$\begin{array}{c}\hfill {\displaystyle {\sigma }_{\gamma A}\left(E_{\gamma }\right)=\sum _i{\sigma }_{i\gamma A}\left(E_{\gamma }\right).}\end{array}$$(29)

It can be shown that the dominant components are E1 and E2, with E2 being (for ions) important mostly at low energies, while M1 is always negligible, so that equation (28) becomes:

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{EMD}}\approx {\int }_{E_{\min }}^{E_{\max }}& [{\sigma }_{E1\gamma A}\left(E_{\gamma }\right)n_{E1}\left(E_{\gamma }\right)\hfill \\ \hfill & +{\sigma }_{E2\gamma A}\left(E_{\gamma }\right)n_{E2}\left(E_{\gamma }\right)]\frac{\mathrm{d}{E}_{\gamma }}{E_{\gamma }}.\hfill \end{array}$$(30)

The equivalent (virtual) photon number is customarily expressed as a function of the adiabaticity parameter ξ, defined as (with the notation ω = E_γ^* and b = impact parameter):

$$\begin{array}{c}\hfill \xi =\frac{\omega b}{ħ\beta \gamma }\equiv \frac{\omega }{{\omega }_{\max }}\end{array}$$(31)

the equivalent (virtual) photon number can be expressed as [193]:

$$\begin{array}{cc}\hfill n_{E1}\left(E_{\gamma }\right)& =\frac{2\alpha Z^2}{\pi {\beta }^2}[\xi K_0(\xi )K_1(\xi )-\frac{1}{2}{\xi }^2{\beta }^2(K_1^2(\xi )-K_0^2(\xi ))]\hfill \\ \hfill n_{E2}\left(E_{\gamma }\right)& =\frac{2\alpha Z^2}{\pi {\beta }^4}[2(1-{\beta }^2)K_1^2(\xi )+\xi {(2-{\beta }^2)}^2{K}_0(\xi )K_1(\xi )\hfill \\ \hfill & -\frac{1}{2}{\xi }^2{\beta }^4(K_1^2(\xi )-K_0^2(\xi ))]\hfill \end{array}$$(32)

where K denotes the modified Bessel functions of the second kind, ${\omega }_{\max }=\frac{ħ\beta \gamma }{b_{\min }}$, and the minimum impact parameter depends on the emitting particle radius (δ_Ruth is a correction, see [191] for details):

$$\begin{array}{cc}\hfill b_{\min }^{\mathit{ion}}& =R_{\mathrm{AB}}+{\delta }_{\mathrm{Ruth}}\hfill \\ \hfill b_{\min }^{\mathrm{lepton}}& =\overline{\lambda }\mathrm{c}=\frac{ħ}{\beta \gamma m_{\mathrm{lepton}}}.\hfill \end{array}$$(33)

It is worthwhile to note in equation (32) the different β dependence of the two multipolarities, which enhances the importance of the E2 component at low energies. As mentioned before, the previous implementation of EMD in FLUKA did not include E2.

The FLUKA EMD model has been extended to the treatment of electronuclear interactions, including the E2 multipolarity and higher-order corrections, as will be explained below. Equations (32) and (33) are partially based on the first Born approximation. Many papers discuss the E1, M1 and E2 virtual photon spectra emitted by e^± including both higher Born approximations and the effect of nuclear/charge finite size (see for instance [194–197]). The numerical calculations of [196, 197] have been used to derive the corrections recently implemented in FLUKA. Accounting for higher Born terms and nuclear finite size effects results in a reduction of the E2 component, a reduction which is larger the higher the nuclear charge Z. For ω close to the emitting particle energy, and for high energies, the first order approach maintains its validity. Despite the E2 virtual photon spectrum being very large for low photon energies, the electric quadrupole contribution is relatively small due to the much lower cross-section for the Giant Quadrupole Resonance (GQR) cross-section when compared with the Giant Dipole Resonance (GDR) one.

The E2 contribution is also important in order to reproduce EMD data for ions at energies of a few hundred MeV/n. For instance, FLUKA predicts that E2 accounts for 25% of the total EMD cross-section for ²³⁸U on ²³⁸U at 120 MeV/n, a value in agreement with other theoretical predictions [198] and consistent with experimental total cross-sections.

FLUKA predictions for the electron and positron-induced fission of ²³⁸U are presented in Figure 40, together with available experimental data. It can be seen that the reproduction of electron-induced fission is very satisfactory, while the positron induced one is slightly over estimated.

4.11. Pre-equilibrium stage in rQMD and BME

The original Boltzmann Master Equation (BME) model used in FLUKA exclusively for ion-ion interaction below 100 MeV/n, and progressively phased out in between 100 and 150 MeV/n did contain a pre-equilibrium stage only for complete fusion reaction, and only for a sub-set of pre-tabulated projectile-target combinations. For the combinations not in the pre-computed database, and for all other channels (3-body, inelastic, nucleon transfer, and incomplete fusion), the excited fragments were passed directly to the FLUKA evaporation routines, with no pre-equilibrium stage. This was creating increasingly inaccurate results when BME was used towards its maximum energies. At the opposite end of the energy range, the improved rQMD model used in FLUKA was going from the QMD/intra-nuclear cascade stage directly to evaporation, without a pre-equilibrium stage. Symmetric to the BME situation, approaching the lower limit of the rQMD energy range, results were not completely satisfactory.

In order to improve both situations, an increasingly sophisticated approach has been developed over the last 12 years which allows the native PEANUT pre-equilibrium model to be used as an intermediate step between cascade and evaporation for rQMD and for all those configurations/reaction mechanisms of BME where no native pre-equilibrium step is available.

Three examples of the improvements due to this approach can be seen in Figure 41, for a pure BME case, in Figure 42, for a case where both BME and rQMD contribute, and in Figure 43, for a pure rQMD case.

4.12. In-flight de-excitation

Until Fluka2020 included, excited nuclear states underwent de-excitation “instantaneously” during a nuclear interaction, by emitting γ radiation only. This approach presents several drawbacks, particularly relevant at the lowest and highest energies treated by FLUKA:

• At high energy even ps/fs mean lives correspond to measurable decay distances, with the further observation that the emitted photon energy is boosted by a factor of the order of γ (γ being the relativistic factor of the emitting nucleus), for example:

  1. Fair (GSI): Au-Au (γ ≈ 10): 1 ns corresponds to a flight path of 3 m, 1 ps to 3 mm;

  2. LHC (CERN): Pb-Pb (γ ≈ 2500): 1 ns corresponds to 750 m, 1 ps still to 0.75 m;

  3. Ultra High Energy Cosmic Rays (UHECR), 0.001–10 EeV/n (γ ≈ 10⁶–10¹⁰): 1 ns corresponds to 0.3–3000 km, 1 fs to 0.3 m–3000 km, excited nuclei can even interact before de-exciting.

• At energies of interest for hadron therapy: the Doppler broadening of both target/projectile emitted photon lines is overestimated, and projectile-like de-excitation photons are emitted at the wrong position. Both artefacts can affect the simulation of therapy monitoring devices based on the detection of prompt gammas:

  1. Target like: ¹⁶O, E = 0.1MeV/n: already for 1ps the flight path is 4.4μm, while the residual range is approximately 3 μm. As a consequence, some excited states will decay at rest, and therefore will not exhibit any Doppler broadening;

  2. Projectile like: ¹²C, E = 150 MeV/n: 1 ps corresponds to 150 μm, 1 ns to 15 cm. Therefore the actual γ emission position for similarly lived nuclear levels can be significantly different from the original interaction point, taking into account the desired sub-millimeter accuracy of hadron-therapy treatments/monitoring.

• In the old implementation Internal Conversion (IC) transitions were neglected. This was in part due to the required availability of atomic electrons for those transitions to occur. Naturally, if fully stripped, an in-flight ion cannot decay by IC, and even if partially stripped its mean life for IC decay will be different. As a consequence:

  1. No IC transitions were performed, rather for those levels decaying both by γ emission and IC the former was always selected;

  2. For strictly forbidden γ transitions (e.g. 0⁺ →  0⁺), unphysical γ’s were emitted.

Now, by default, excited nuclei with measurable mean lives do not de-excite during the nuclear interaction that produced the excited state, but rather they fly until decay according to the mean life of the level. At the same time, IC competition is fully accounted for, and, if selected, IC transitions can occur only when the excited ion comes to rest fully neutralized.

4.13. Fission

The fission model has been improved in order to better reflect the competition between symmetric and asymmetric channels and to improve the shape of the fragment mass distribution. The improvements are particularly visible at energies above few tens of MeV. Figure 44 shows old and new results for 60 MeV protons on ²³⁸U in comparison with experimental data. At lower energies, the model continues to compare well to available data, as shown in Figure 45 for fission induced by low-energy neutrons. Moreover, a careful tune of the Coulomb energy back-shift achieved an almost perfect agreement with data on prompt neutron multiplicity (see Tab. 1). As far as low-energy neutrons are concerned, a better reproduction of inclusive spectra and multiplicities would be achieved by sampling directly from the wealth of ENDF available data. This is what FLUKA does when the group-wise neutron treatment is invoked. However, there are two reasons to develop and apply a native fission model for neutrons below 20 MeV. The former is that this algorithm allows the simulation of correlated exclusive events, as described in Section 5.2. The latter is that reactions at higher energies, or initiated by other projectiles, can result in a reaction product with A, Z, and excitation energy compatible with those of a low-energy neutron interaction: a reliable description of its fission is needed and is validated by the low-energy neutron data.

5. New developments in low-energy neutron transport

5.1. Updated group-wise transport

Starting with Fluka2020, neutron cross-sections for several isotopes have been updated with more recent evaluations, mainly ENDF/B–VIIIR0 [179], and more isotopes have been included in the FLUKA library (e.g. those of Os, Tl, Ho, Dy, Lu, Te, Rh, Se, Rb, and ¹³C, ¹⁷O, ¹⁸O, ²³⁷Np, ²³⁹Np, ²³³Pa, ²³¹Pa, ²²⁷Ac, ²²⁶Ra).

At the same time, fission fragment yields from low-energy neutrons have been updated with the latest evaluations, ENDF/B–VIIIR0 [179], JEFF–3.3 [204], and JENDL–4.0 [205]. Updates are ongoing and they will soon include the recently released JENDL–5 evaluation [180]. Recently, the neutron library has been supplemented with the information on isomer production contained in the European Activation File (EAF) [206]. Therefore, the production of isomers versus the ground state by low-energy neutrons is calculated according to the evaluated data when available.

5.2. Neutron point-wise, fully correlated cross-sections

The FLUKA point-wise library is built by processing the evaluated data with the PREPRO19 [207] and NJOY2016 [97] codes, complemented by a FLUKA-specific code that converts to the FLUKA format and takes care of correlations, as discussed in the next section. The FLUKA library includes all stable isotopes, some important unstable ones (e.g. ¹³⁵Xe) and the main transuranic isotopes. Data are pre-processed at different temperatures: 4, 87, 296, 430, and 686 K. Most isotopes are retrieved from ENDF/B-VIIIR0 [179], and a few from TENDL-19 [173]. The use of NJOY allows for the correct accounting of unresolved resonances (URR). These resonances, which are not resolved experimentally, induce fluctuations in the neutron cross-section and are typically situated between 1 keV and 1 MeV. NJOY prepares URR probability tables to be further employed by transport codes.

A careful validation of the FLUKA point-wise library has been carried out by comparing it with group-wise transport for each and every isotope. The test geometry consists of three concentric spheres, two made of the material under consideration, and the third, in the middle, made of water. 20 MeV neutrons are injected at the centre of the spheres, and neutron spectra are collected at all boundaries. An example of the excellent consistency between the two methods is shown in Figure 46, together with a schematics of the geometry. Figure 47 instead demonstrates the necessity to include a correct treatment of URR, in the absence of which the point-wise transport shows clear anomalies.

Data on the response of Bonner spheres to monochromatic neutron beams at different energies [208] provide another benchmark. As shown in Figure 48, the excellent agreement between data and FLUKA simulations already obtained with group-wise transport in [208] is replicated with point-wise transport.

Moreover, as expected, the point-wise treatment automatically accounts for self-shielding according to any given material composition.

The main advantage of group-wise transport is speed. Nevertheless, the new point-wise correlated algorithm only imposes a significant CPU penalty in worst-case scenarios involving low-energy neutron transport only, with full transport of the charged products down to low thresholds. CPU times for typical full shower simulations are marginally affected [211].

5.2.1. Event-by-event correlations

Evaluated data files contain uncorrelated and inclusive information. This means that, when using this information in Monte Carlo sampling, there is no guarantee that energy and momentum are conserved, neither in single interactions nor in complete “histories”. For example, the sum of photon energies after a neutron capture may exceed the capture Q-value, just because many photon energies are sampled independently from the same spectrum, each one without knowledge of the previous ones. Energy conservation is conserved on average, as well as inclusive particle spectra are reproduced on average. While most shielding applications are insensitive to this aspect, applications for particle detector studies can be severely affected.

In addition, nuclear recoil spectra may not always be available or reliable, which can have significant implications for the accuracy of radiation damage calculations.

Pioneering work to include correlations has been done in PHITS [212]. Starting with the Fluka2021 release, FLUKA provides fully correlated point-wise neutron interactions. For each interaction, the outgoing products, including the residual nucleus, are sampled to ensure that the cross-sections and spectra of the evaluations are reproduced as closely as possible while conserving energy and momentum. To this purpose, additional information or nuclear models are required to complement the inclusive data sets available in the evaluated data files.

Gamma de-excitation after interactions is always performed using the FLUKA de-excitation algorithm described in Section 3.2.9.

For “simple” reactions such as (n,el), (n,n′), (n,nα), (n,np), constrained sampling from available spectra and angular distributions are performed. Adjustments are made in the event of un-physical inconsistencies in the data, with priority always being given to the reproduction of the neutron spectrum.

The (n,2n) and (n,3n) reaction represents the most complex task. A sophisticated offline algorithm splits the evaluated inclusive neutron spectra into 1^st and 2^nd chance neutrons (3^rd for (n,3n)) such that the overall inclusive spectrum is conserved, and the reaction fulfils energy and momentum conservation. Thanks to this algorithm and associated checks several inconsistencies in some of the evaluated (n,2n), (n,3n) spectra have been found and corrected altering the evaluated spectra as little as possible so that they fulfil basic conservation laws.

For neutron capture, the branching ratios of the primary gamma transitions (those from the capture level to an underlying one) are taken from the Evaluated Nuclear Structure Data File (ENSDF) database [213], supplemented by the FLUKA statistical model [81] if the data are incomplete. The following steps of the γ cascade are computed in agreement with the data collected in [82]. Improvements in the treatment of incomplete primary gamma transition branchings are underway.

For more complex reactions such as (n,2p), (n,4n) and (n,2nα), the reaction channel is selected according to the ENDF data, and the reaction products are calculated with a Monte Carlo method based essentially on phase space considerations, similar to the Fermi break-up FLUKA model as in Section 4.7. The initial selection is much easier in cases where all reaction channels are clearly distinguished in the evaluated data. Cases where the reactions are grouped in the so-called (n,complex) (MT5) require support from models or correlation with residual nucleus production to disentangle the different channels.

Fission is simulated using the FLUKA Evaporation Model, forced to proceed through fission. This allows us to obtain correlated fission products and neutrons for each event while reproducing the fission product mass distributions as in the evaluated data.

The power of full correlations is evident in Figure 49, which shows the response, in terms of deposited energy, of a high-resolution diamond detector to a quasi-mono-energetic neutron beam. The device has the double role of target and detector. It was exposed [209] to neutrons in the range 18.9–20.7 MeV with an intrinsic energy spread of 0.2–0.25% FWHM. Due to its small dimensions (4.5 × 4.5 mm², thickness 500 μm) secondary neutrons and photons mostly escape without depositing energy. Several peaks are present in the experimental pulse height spectra, corresponding to reactions where only charged particles are produced. The same peaks are visible in the simulation. In Figure 49 the response to 20.49 MeV neutrons is shown and compared with experiment [209]. In the same figure, the contributions to energy deposition of the different reaction products are depicted separately. Peaks can be reproduced only if the energy of reaction products is generated in a correlated way in each interaction, adding up to the projectile energy minus the Q-value. The same is true for the shape of the continuum, corresponding to events where neutral particles escape the system.

Figure 50 illustrates another example. It shows the response of a Silicon Carbide (SiC) detector irradiated with 14.12 MeV neutrons [210]. Peaks corresponding to reactions on carbon and silicon are evident in data and simulations, particularly the one from the ¹²C(n,α₀)⁹Be reaction and those from ²⁸Si(n,α_i)²⁵Mg on various excited states of ²⁵Mg. The treatment of non-ionising energy losses affects the simulated response, as shown by the difference, small but visible between the red and green curves in Figure 50. The green curve is constructed by summing all charged particles’ energy losses, while in the red curve, only the ionizing fraction of the stopping power is considered, the non-ionizing one is discarded. The effect is larger for reactions on silicon, due to the higher mass of the recoiling nucleus, and improves the agreement with experimental data. This confirms the power of the new correlated treatment that allows faithfully reproducing the recoil spectra, thus investigating the detector response to the few percent levels. Residual discrepancies may point to residual experimental effects, or to the need to better model the partition between ionizing and non-ionizing losses, that here is performed only in an average way, without fluctuation, and with an approximate treatment of recoils below the transport threshold.

5.2.2. Point-wise thermal scattering laws

Starting with the 2024.1 release, FLUKA now supports point-wise neutron transport at thermal energies, also for materials with available Thermal Scattering Laws (TSL). In particular, the FLUKA implementation supports continuous energy and angular distributions for incoherent inelastic, aka S(α,β), coherent elastic, and incoherent elastic processes. Pre-processing of the ENDF files is performed using both ACEMAKER [214] and NJOY2016 [97]. The resulting ACE files are then post-processed by a FLUKA-specific tool in order to fix possible inconsistencies and to transform the discrete angular distributions for S(α, β) into continuous ones. For incoherent elastic scattering, only the total cross-section and the Debye–Waller integral are kept, and the angular sampling is done analytically without any discretization.

Figure 51 shows the impact of using the free gas treatment rather than the thermal scattering kernel for hydrogen in polyethylene when computing the response function of one of the Bonner spheres previously presented in Figure 48.

Figures 52 and 53 compare energy spectra and direction cosines of scattered neutrons in the broomstick benchmark tests [215, 216] for para-hydrogen at 20 K and 10% porosity reactor graphite at 296 K, computed with FLUKA and MCNP [217]. In this setup, neutrons are sent into an infinitely long and extremely thin cylinder of a given material along its axis, and the energy and angular distribution of the neutrons emerging from the cylinder are scored. The comparisons presented in Figures 52 and 53 are examples of the extensive verification performed on FLUKA’s newly implemented point-wise thermal scattering kernel algorithm at different energies with various materials, both versus results obtained with other codes and results obtained with the group-wise library.

6. Recent developments in the Electro-Magnetic FLUKA (EMF)

EMF is very stable and its physics models are mature and extensively benchmarked for several years. Only a few incremental additions/improvements have been implemented in the last decade which are described in the following.

6.1. New photon cross-section database

The photon atomic interaction cross-section database has been updated and it is now based on EPICS2017 [99]². Among the other advantages of the updated database, the individual shell binding energies are now consistent with currently accepted values hence resulting in fluorescence X-ray energies slightly different from those of the previous database, based on EPDL97.

6.2. Extended bremsstrahlung database

The bremsstrahlung cross-sections and spectra are generated offline by a dedicated, FLUKA specific, program which produces a database for all elements and primary and secondary energies. In order to allow the code to be used at the highest cosmic ray energies, the treatment has been recently extended to 100 EeV.

6.3. Anomalous form factors in coherent scattering

The coherent (Rayleigh) scattering cross-section previously listed in EPDL97 and now in EPICS2017 is evaluated taking into account also the anomalous real and imaginary form factors, besides the atomic elastic form factor. Up to Fluka2023.3, FLUKA did not consider the extra terms due to the real, f₁, and imaginary, f₂, anomalous form factors when computing the photon angular distribution following a coherent scattering event, and it was using only the atomic elastic form factor for that purpose. Starting with Fluka2024.1 the f₁ and f₂ related terms are now accounted for when sampling the photon angular distribution in coherent scattering. The availability of f₁ and f₂ is further exploited when computing the complex X-ray refraction indices involved in the X-ray reflectivity (see Sect. 6.5).

6.4. Ortho- para-hydrogen competition for positron annihilation at rest

Ortho-Positronium 3-photon annihilation is now included in FLUKA in competition with the para-positronium, 2-photon annihilation. A positronium e⁺e⁻ bound system is formed in about 40% of positron at-rest annihilations, in the residual 60% the positron annihilates directly into 2-photons without going through a positronium intermediate state. Out of the 40% positronium occurrences, roughly 10% are initially formed in the S = 0 para-positronium sate, which readily annihilates in to 2-photons with a mean life of approximately 0.12 ns. The other 30% are formed in the S = 1 ortho-positronium state which has a much longer mean life of about 140 ns. During this time the ortho-positronium can also decay into 2-photons by picking up a nearby electron with opposite spin. At the very end only a small fraction of the initially formed ortho-positronium states will eventually decay as such into 3-photons. Overall the probability of 3-photon versus 2-photon annihilation is 1:378. There is some evidence that the Ortho/Para ratio and Ortho decay constant could be material dependent. Such an effect has a direct application in the Positron Emission Tomography of ipossic tissue. The FLUKA implementation describes 3-photon annihilation with the theoretical matrix element, and taking into account atomic electron motion and binding. The ratio ortho/para annihilation is tunable per element and region/material, as well as the decay time constants, allowing for investigations about the impact of different ortho-positronium rates in different materials/tissues.

6.5. X-ray reflectivity

It is now possible to reproduce the reflections of photons at interfaces by activating a material-specific option, see Figure 54 for an example. The single-layer mirror reflections with a coating material laid upon a substrate are implemented in FLUKA by computing permutated expansions of recurrent relationships [218, 219]. For complex interfaces with multiple layers, a user routine needs to be prepared but, for reflections from thick mirrors, the code has the relevant information for all elements and it can automatically calculate the reflection probability by resorting to EPICS tabulated data. For setups other than the thick mirror case the refraction indices are computed internally by FLUKA and can be easily retrieved in the user routines.

7. Recent developments in particle transport

7.1. Transport in electric and magnetic fields

Starting with Fluka2021.2 issued in May 2021, FLUKA incorporates the possibility of accurately modelling the motion of charged particles in a gas where electric fields are present, such as inside ionization chambers or other detectors. A volume virtually “filled” with vacuum, such as the beamlines of a particle accelerator, would be a special case, also covered in the model. This is in addition to the long-existing capability to simulate transport under purely magnetic fields of arbitrarily geometrical complexity.

Since in the presence of electric fields charged particles see their kinetic energy altered, implementation of step-wise tracking

required a more computationally efficient and precise method than the one existing for magnetic fields, which was based on arcs steps sufficiently small to adapt to the geometric boundaries or the

magnetic field spatial variations. The Runge–Kutta–Gill [220] algorithm, which breaks steps into four sub-steps recursively computed as a function of

local derivatives was found to perform very satisfactorily when benchmarked against numerous analytical equations. For example, for

electrons of 0.1 MeV/c uniformly accelerated along a 1 m tube with a field ranging from 0–20 MV/m, the computed kinetic energy at the exit was less than 10⁻⁶% off from the analytical value. Similar tests carried out for complex, fully analytical 3D fields generated by a combination of charged planes, cylinders and spheres, yielded similar accuracies even when a magnetic field is also present. Yet Euler numerically resolved equations required much smaller steps to match FLUKA results, showing the efficacy of the Runge–Kutta–Gill approach as implemented in FLUKA. As for combined electric and magnetic fields, FLUKA was found to correctly reproduce the cycloidal motion of a particle subject to orthogonal electric and magnetic fields, matching well the expected period.

Charged particles in electric fields may decelerate, and even come at rest and reverse their direction. Consequently, FLUKA was adapted to perform EM tracking below the threshold energies set for transport in matter. Moreover, the existing TIME-CUT can be used to interrupt the tracking of charged particles trapped in an electric potential well, as well as low-energy particles spinning indefinitely under magnetic fields.

An important feature that was included in the functions for electric field tracking, and which was also added to those for magnetic fields, is the ability to define time variations. This has enabled the implementation of AC fields, including those in high-frequency applications, such as high-energy particle accelerators (see Fig. 55). With tools like the FLUKA LineBuilder [221] or MadFLUKA [222] that help automate the computer implementation of particle accelerators from their building blocks through FLUKA lattice and roto-translation directives, and with the new capability to include also accelerating cavities with their RF fields, etc., FLUKA can now practically simulate entire accelerators, possibly just excluding injectors with high space-charge effects. This might open the door to computing beam losses and the propagation of their induced showers all in one. For example, LCLS–II cryomodules have been modelled in detail, including 3D maps for their 1.9 GHz electromagnetic fields. This has allowed “test beams” to successfully synchronize with the RF fields, accelerating to the expected energies, and even has reproduced off-peak acceleration for phase shifts.

When FLUKA computes the electromagnetic field at a given position and time, the result is typically just used internally to steer charged particles, but there are other potential applications. For example, the emission of electrons from metal surfaces subject to strong electric fields has been recently modelled for the cavities of the aforementioned LCLS–II cryomodules, where the subsequent trajectory of the field-emitted electrons within the EM fields of the cryomodules was then computed making use of the new FLUKA capabilities.

Just like for magnetic fields, electric fields can be set to apply for all charged particles, just for prompt radiation, or only for residual radiation. This is an example of how FLUKA users who are familiar with the use of magnetic fields should be able to add electric fields without much difficulty, as the process is very analogous, in terms of input cards, region applicability, and user routines.

7.2. Synchrotron radiation

FLUKA can generate and transport synchrotron radiation photons produced by a charged particle travelling in a magnetic field via a special predefined source routine. The latter allows for a detailed description of the emitting particle (which itself is not transported), the respective circular arcs or helical paths, and the magnetic field properties to accurately reproduce the synchrotron radiation emission. (An example is in [Fig. 56]).

Polarization is implemented as a function of emitted photon energy and angle, and its effects are included for Compton, Photoelectric, and Coherent scattering while taking into account bound electron effects.

8. Examples of FLUKA applications

8.1. Medical applications

FLUKA is widely used in medical physics, in particular in hadron therapy. Below we discuss several applications of the FLUKA code related to medical physics.

8.1.1. Charged hadron therapy

FLUKA plays a pivotal role in supporting charged hadron therapy, particularly as a key component of the clinically-oriented workflow at the Heidelberg Ion-beam Therapy centre (HIT) and the National Center for Oncological Hadrontherapy (CNAO), in Europe. The demanding nature of these medical applications has been tackled by continuous refinement of nuclear models and extensive improvements in FLUKA modelling for cross-sections and interactions of light ions, to meet stringent reliability and predictability criteria [226]. Figure 57 shows the excellent reproduction of experimental depth-dose profiles for carbon ions.

At HIT and CNAO, for example, FLUKA has contributed to the generation of the synchrotron library and the derivation of experimentally validated datasets for the treatment planning systems (TPS), which are now essential components of clinical practice in proton and light ion therapy [224, 227].

Additionally, FLUKA has proven its efficacy in bolstering clinical applications, as MC dose calculations continue to serve as a valuable tool for patient treatment verification, supporting dosimetric measurements [228]. For computations based on computed tomography (CT) imaging, the patient model relies on stoichiometric calibration, with adjustments made to account for facility and CT-number dependencies of electromagnetic and nuclear processes, ensuring alignment with the CT-range calibration curves utilized by the TPS across all CT protocols [229, 230].

Both physical and relative biological effectiveness (RBE) weighted dose distributions can be computed for individual treatment beams and overall treatment fractions, using either constant or variable RBE models for protons and incorporating the major biological models adopted in the hadrontherapy community for other particle species: the local effect model (LEM) [231], the micro-dosimetric kinetic model (MKM) [232], and the repair-miss repair-fixation model (RFM) [233]. In physical and biological computations, dose-to-medium can be dynamically converted into dose-to-water, facilitating dose distribution analysis in both frameworks. Optionally, dose-averaged linear energy transfer (LET) distributions can be generated. Thanks to this feature, FLUKA has enabled the validation of treatment protocols at CNAO for carbon ion radiotherapy, based on the extensive clinical experience of the National Institute for Radiological Sciences (NIRS) [234]. FLUKA has also been extensively used in the evaluation and application of different RBE models for proton therapy [235].

As for proton and carbon ions, FLUKA has become the gold standard for helium ion therapy at HIT due to the high level of agreement found comparing FLUKA against dosimetric data in simple and complex conditions (e.g. with anthropomorphic scenarios, see [225, 236]). In fact, three decades after the end of the Lawrence Berkeley National Laboratory (LBNL) clinical program with helium ion beams, the first patient treatment with raster-scanned helium ions was performed in July 2021 at HIT. Over the past few years, several works have investigated the potential of helium ions using FLUKA. In Figure 58, the predictions of FLUKA are compared against experimental data for double differential neutron spectra emitted by the interaction of 135 Mev/n ⁴He ions with carbon [199]. The sound prediction of the mixed radiation field generated by ⁴He ion beams together with the well-benchmarked atomic models of FLUKA result in a solid prediction (within few %) of the depth-dose distributions in water (see for example Fig. 59). Along with the renewed interest in the therapeutic potential of ⁴He ions, also their promising role as a possible future imaging modality replacing X-ray CT for treatment (re)planning has been thoroughly explored with FLUKA simulations in comparison to protons and carbon ions. Investigations carried out using an ideal single particle tracking detector and the realistic HIT beam properties showed advantages of helium ion imaging in terms of the lowest error in retrieval of tissue stopping power ratio (relative to water) at the same imaging dose for a human head [233].

In a broader context, FLUKA offers the flexibility to interface with various models adhering to the linear quadratic formalism of the cell survival theory, even including the incorporation of the oxygen enhancement ratio (OER) [237].

The FLUKA coupling with BIANCA [238], a biophysical model of cell killing and chromosome damage, represents, in that sense, a practical example easily generalizable. Whenever according to FLUKA a certain amount of energy is deposited in a target voxel by a given particle type of given energy, FLUKA reads the corresponding linear and quadratic coefficients describing the cell survival probability as a function of energy/LET from tables produced by BIANCA. Afterwards, to take into account that a mixed radiation field is present in the voxel, dose-averaged values of the two coefficients are calculated as described in [239]. This allows for calculating the cell survival probability and the RBE-weighted dose in that voxel. As described in [240], an analogous approach allows calculating yields of non-lethal chromosome aberrations, which are correlated with cell conversion to malignancy and thus can be used to estimate the probability of cancer induction in healthy tissues. After benchmarking against in vivo experimental data on carbon ion and proton irradiation of the rat spinal cord [239], this approach was applied to re-calculate RBE-weighted dose distributions for carbon ion irradiation of cancer patients previously treated at CNAO, showing that BIANCA is suitable for treatment plan optimization in ion beam therapy [241], including irradiation with helium ions [242]. Furthermore, following the extension of BIANCA to heavier ions up to Fe [243], the same approach was applied in the framework of astronauts’ radiation protection, including the prediction of Galactic Cosmic Ray doses during a mission to Mars [244] and of organ equivalent doses in case of exposure to an intense Solar Particle Event [245].

8.1.2. Developments of PET simulation tools

Since its development in the second half of the 20^th century, positron emission tomography (PET) technology has become well-established in nuclear medicine. Modern PET scanners are usually coupled with other technologies, such as computerized tomography (CT) or magnetic resonance imaging (MRI), to complement the information about the location of β⁺ emitters with other anatomical or metabolic processes. The development of proton and heavy ion therapy treatments has forced PET technology to move beyond its initial diagnostic applications into monitoring and range verification in an extremely noisy environment. The simulation of such extreme backgrounds with Monte Carlo codes such as FLUKA is therefore highly recommended. Indeed, several studies already analyzed the feasibility and reliability of the FLUKA code for dose prediction and hadron therapy treatment monitoring by PET scanners [246–249] prior to the development of the FLUKA PET tools described in this section.

The FLUKA PET tools [250–252] are a dedicated package aimed at exploiting FLUKA full predictive potential to reproduce the PET/CT scanner response for diagnostics imaging and hadron therapy monitoring. These tools can provide images reconstructed from a PET signal generated by a predetermined radioactive source or beam irradiation as acquired in a modelled PET scanner, as well as user-customizable simulation data for further analysis.

Currently, thanks to a user-friendly graphical interface, a custom PET scanner can be assembled directly in the FLUKA input with few parameters or by selecting some of the templates available in the FLUKA advanced interface (FLAIR [223]). The detector’s geometry elements are interpreted by scoring routines to acquire the signals generated throughout the simulation. A general output is then recorded, listing various details of the interactions with the scanner, which can be further analyzed by the user. At the post-processing stage, another set of routines converts the output information into a rearranged list of coincidences and eventually into sinograms (including separated files for true, random, and scattering coincidences).

These tools are mainly intended to aid in developing new scanners or to help validate treatment scenarios in hadron therapy, particularly at the dosimetry and range monitoring level, where the full description of the particle background can be challenging.

8.1.3. Application to range monitoring of particle therapy treatments

In the last decades, FLUKA has been extensively used in the context of range monitoring of particle therapy treatments, using positron emitters, prompt gammas and secondary particles, as described in [253]. An early MC study, including several types of secondary particles (prompt gammas, neutrons, protons and ions, β⁺-emitters), was performed by Robert et al [254]. They compared FLUKA with GATE/GEANT4, including several distributions of various types of secondary particles.

Non-invasive range monitoring with β⁺-emitters was studied repeatedly with FLUKA. For proton therapy, offline PET treatment monitoring was performed in a pilot clinical study at the Massachusetts General Hospital to assess the feasibility of offline in-vivo treatment verification by Parodi et al [247]. FLUKA was combined with external cross-sections to obtain the predicted PET images, and a mathematical approach was used to model biological washout. The same approach was applied later on by Knopf et al. for a larger patient population [255], where the effect of the tumour location on the accuracy of the technique was assessed. Regarding carbon ion therapy, at HIT, the offline PET technique was implemented and the first clinical cases were presented [256]. Here predictions simulated using the FLUKA internal nuclear models were compared with offline PET measurements acquired with a commercial full-ring PET/CT scanner, installed in close vicinity to the treatment rooms. The simulation included the production of β⁺-emitters, a post-processing procedure considering the time course of the respective treatment fractions, and a simplified model for estimated biological washout of the induced activity. These studies demonstrated the feasibility of post-radiation PET/CT for in-vivo treatment verification (See [black)]Fig. 60).

The interest in in-beam PET has grown in the last decade, and several studies on in-beam PET have been performed with FLUKA relying entirely on its sound nuclear models. For example, in-beam PET monitoring of mono-energetic ¹⁶O and ¹²C beams was performed at GSI [248]. FLUKA was used to simulate the production of the β⁺-activity by all ions of interest. In particular, ion transport, β⁺-active isotope production and decay, positron annihilation and the transport of the annihilation photons were simulated with FLUKA. In-beam PET measurements of the activity created by ¹⁶O beams of various energies on targets of PMMA, water and graphite were thoroughly compared with FLUKA predictions. More recently, FLUKA was used in pre-clinical offline measurements by the DOPET [257] system, acquiring phantom data shortly after irradiation [258–260]. FLUKA was also used by the INSIDE [261] bi-modal imaging system discussed in the following.

Range monitoring with prompt gammas in proton and carbon ion therapy was also extensively studied with FLUKA. For example, reference [262] performed a feasibility study for a pixelated prompt gamma imaging detector which can not only measure the proton range but also merge the beam profile and 3D beam image inside the patient. Moreover, a simulation study was performed to show whether it was possible to measure the prompt-gamma fall-off inside the high neutron background typical of carbon ion irradiation, using a knife-edge slit camera [263]. FLUKA was also applied to evaluate the potential of secondary neutrons to provide proton range information, possibly supplementing the signal from prompt gammas [264].

In the context of treatment monitoring with secondary charged particles, a recent work [265] used FLUKA to evaluate the performance of secondary ion tracking to detect, visualize, and localize an internal air cavity used to mimic inter-fractional changes in the patient anatomy at different depths along the beam axis. The Timepix3 semiconductor pixel detector technology has been fully simulated [265]. This work also contained data-simulation comparisons.

Finally, a bi-modal imaging system (INSIDE) was developed and tested at the CNAO treatment facility. A clinical trial was started, in which in-beam PET and secondary proton data were acquired during and shortly after proton and carbon ion patient treatment. The experimental setup was simulated in detail with FLUKA, including the CNAO beam line, spot scanning system, production of β⁺ emitters and secondary particles, patient geometry, and so on. The FLUKA simulation framework was extensively used for PET-data analysis feasibility studies [266–270]. Moreover, it was used for studying secondary proton distributions (multiplicity, direction, and emission position along the beam) [271, 272].

In addition to stable beams, radioactive ion beams were also studied with FLUKA in an in-beam PET framework, with experimental validation of FLUKA dosimetry and β⁺-emitter predictions of ¹¹C, ¹²C, ¹⁵O and ¹⁶O ion beams against experimental data obtained with an openPET prototype in various online PET acquisitions at the Heavy Ion Medical Accelerator in Chiba (HIMAC), in collaboration with the National Institute of Radiological Sciences (NIRS) Imaging Physics Team [273]. A dedicated estimator for the annihilation events at rest was also implemented in FLUKA, allowing the user to directly quantify and obtain the spatial distribution of β⁺-emitters’ activity, integrated in time or at different steps during or after irradiation. Moreover, the code can track the parent isotope that originated the β⁺ particle annihilated. More recently, FLUKA was used by the BARB (Biomedical Applications of Radioactive ion Beams) experiment at GSI (Darmstadt) to perform feasibility studies to investigate PET imaging for high-intensity ¹⁰C and ¹¹C beams, showing that range monitoring with sub-millimeter resolution was possible [274].

In addition to all these works, exploiting secondary emissions produced in nuclear interactions, FLUKA has also been used in connection with explorations of range verification by means of acoustic waves induced by ionizing radiation. Here, FLUKA was used to provide the initial energy deposition required for input into an external engine for acoustic wave propagation. Drawing upon its novel ability to generate and transport optical photons, the putative contribution of an additional photoacoustic effect observed when using contrast agents with optical absorption properties was investigated [106].

8.1.4. Application to FLASH radiotherapy

Recent pre-clinical studies suggest that therapeutic doses delivered at ultra-high dose rates (UHDR, average dose rate > 40 Gy/s) can substantially reduce toxicity in healthy tissues while maintaining the same cancer-killing efficacy compared to doses delivered at conventional dose rates (approximately 0.05 Gy/s). This phenomenon, referred to as the “FLASH effect”, has the potential to enhance the therapeutic index of radiotherapy and has received substantial attention in the field of radiation oncology in recent years [275–278]. The ease of producing and delivering electrons at UHDR makes their use for FLASH radiotherapy particularly intriguing.

To explore the clinical potential of very high energy electron (VHEE) beams (50–300 MeV) [279, 280], it is essential to develop treatment planning systems (TPS) capable of planning and optimizing VHEE radiation delivery to patients. Such a VHEE TPS was developed based on FLUKA’s evaluation of dose matrices. The system was then used to determine the optimal accelerator parameters for treatment delivery, including the energy, intensity, and direction of therapeutic radiation fields. Detailed information on VHEE beam properties and the potential of VHEE-based radiotherapy at both conventional and ultra-high “FLASH” dose rates are reported in [279–283].

Additionally, dedicated VHEE devices are necessary to further characterize the FLASH effect for UHDR VHEE beams in pre-clinical settings and to assess the feasibility and added value of VHEE-based FLASH radiotherapy in clinical settings. In this context, a detailed FLUKA simulation of the prototype VHEE LINAC under construction at La Sapienza University was developed. This simulation was coupled to the output of the beam dynamics simulations, to evaluate the stray radiation induced by the interaction between the beam and the accelerator components for radioprotection purposes [284]. FLUKA was employed to evaluate scattered radiation, optimal shielding strategies, and the impact of continuous machine operation in a hospital environment. Figure 61 shows an example of a VHEE irradiation beam from the Sapienza VHEE LINAC prototype impinging on a water phantom.

8.2. Astrophysical applications

The FLUKA code has been applied to astroparticle physics since the second half of the 90’s. Probably the first published work in this respect was a numerical study of the electromagnetic component of Extensive Air Showers (EAS) in the atmosphere, up to 100 TeV [285]. There it was mentioned for the first time the use of a stratified earth atmosphere, made of several layers, each one having uniform density, scaling with height according to a fit of standard US atmosphere commonly used in cosmic ray physics at that time. In this work, the curvature of the earth was still neglected. Later, the FLUKA code was applied to problems related to the detection of high-energy muons underground. In [286] the issue of transport codes at high energy in the rock was addressed, mainly investigating the problem of possible incorrect approaches in the models for radiative energy loss of muons in use at that time. However, one of the main successes of FLUKA applications in the physics of secondary cosmic rays was the first 3-dimensional calculation of atmospheric neutrino fluxes [287]. In this work, the previously mentioned approach of a stratified atmosphere was improved and extended to a spherical earth model. The interaction of primary cosmic rays in the atmosphere was fully simulated, using in input the all-nucleon spectrum already adopted in previous 1-dimensional calculations by the group of Bartol Research Institute, considered as one of the main references at that time [288]. The treatment of solar modulation was addressed, together with the effect of geomagnetic cut-off for a few experimental sites in the world. This calculation was later improved by using an updated primary cosmic ray spectrum which took into account the most recent experimental measurements [289] and then was extended to consider neutrino energies below 100 MeV [290]. This calculation is still of particular importance today since these low-energy neutrinos contribute to the irreducible background for the underground experiments devoted to the search for dark matter particles and the search for the diffuse supernova neutrino background[291]. All these FLUKA-based calculations of atmospheric neutrino fluxes needed to be validated. For this reason, the fluxes of muons generated in the same interaction of neutrinos were produced with the same simulation setup and compared with the existing experimental measurements. In this respect, the first work was the successful comparison with the 1994 data by the CAPRICE balloon experiment [292]. Along this line, further validations can be found in [293, 294] using the muon data from the BESS experiment. As far as higher muon energies are concerned, validation against the measurement of charge ratio in the atmospheric muon data collected by the L3 experiment at CERN was also presented [295]. Based on the experience acquired thanks to these achievements, a complete FLUKA-based package for cosmic ray simulation has been developed. This is described in chapter 16 of the FLUKA user manual [296].

The recent improvements in astroparticle (gamma rays, cosmic rays and neutrinos) detectors have brought a precision era for the study of the propagation, interactions and production of these particles in galactic and extragalactic environments [299, 302]. The current accuracy on cosmic-ray (CR) spectral data allows us to precisely test our models of propagation of charged particles in the Galaxy as well as to improve our predictions for the emission of secondary particles, such as neutrinos and gamma rays, which are produced from the interactions of CRs with the interstellar medium [303]. However, these studies are severely limited by the uncertainties related to nuclear interactions and, in particular, to the cross-sections for CR interactions [304–307], which are currently measured with very poor accuracy (20%), and often no data are available at the energies of interest (from the GeV to the TeV).

In recent years, FLUKA has been then applied to a series of different astroparticle physics studies [148, 308–313] in the past, raising high expectations in the CR community. Recently, the authors of reference [151] presented the results of a campaign dedicated to the computation of the full cross-sections (including inelastic and inclusive cross-sections) of the nuclei relevant for Galactic CR studies and that are needed to be implemented in CR propagation codes using FLUKA. This involves the calculation of cross-sections of all nuclei from protons up to iron (including all isotopes in between and considering both stable and short-living species) interacting with ¹H and ⁴He from a few MeV to tens of TeV. The result is a set of tables with the production cross-sections of all isotopes involved in the CR network, as well as secondary positrons and electrons, gamma rays, neutrinos, and antiprotons that will be made publicly available soon.

The authors of reference [151] investigated the use of the FLUKA cross-sections for the prediction of the propagation parameters characterizing the transport of CRs using the secondary CRs B, Be and Li. We show in Figure 62 the good agreement of the evaluated local interstellar spectra of the main primary CRs (left panel) and the ratios of secondary-to-primary CRs (right panel) with available experimental data. This study demonstrated that the total inelastic cross-sections derived from FLUKA from the MeV to the TeV region are in good agreement with the current data and compatible with elaborated parametrizations to the level of a few percent for nuclei with Z ≤ 10 (up to Ne), where there is no experimental data available in the GeV energy range. For heavier nuclei, the energy trend of the FLUKA inelastic cross-sections was in agreement with parametrizations but approximately 15–30% lower in absolute value. This may indicate that the extrapolation of these parametrizations for nuclei for which there is no experimental data may be biased toward an overestimation of the inelastic cross-sections.

In addition, the inclusive cross-sections involved in the production of the isotopes of Be, Li and B were compared to the data and the most recent dedicated parameterizations, showing a good general agreement, with discrepancies lower than 30%. This represents a remarkable success of the FLUKA calculations since previous attempts to calculate cross-sections for CR interactions have resulted in very inaccurate predictions that usually can be only used to have an order of magnitude estimations. Finally, this work demonstrated that the FLUKA interaction model reproduces simultaneously the spectral data on B, Be and Li and ³He that has significant consequences on the physics behind the propagation and production of these nuclei. Therefore FLUKA can be safely used as an alternative solution to the use of cross-sections parameterizations, which totally depend on the very scarce, limited and uncertain experimental data available for these interactions.

FLUKA was used in reference [151] to predict the diffuse gamma-ray emission in the Galaxy. The derived gamma-ray maps were shown to be consistent with the data and, remarkably, the authors demonstrated that the local emissivity on atomic hydrogen (HI emissivity) can be well reproduced from the models derived with the FLUKA cross-sections, including a low-energy break in the injection spectrum of electrons.

Similarly, FLUKA was used to evaluate the production of electrons and positrons from CR interactions [314]. In this work, the authors investigated the main uncertainties affecting our predictions on the local positron spectrum and showed that the FLUKA estimation is compatible with the most recent estimations of e⁺, e⁻ production cross-sections within an approximately 30% relative difference. Also discussed the potential of using FLUKA for improving the current galactic magnetic field model and for the prediction of low-energy positrons (in particular, for the study of the 511 keV signal), as well as the importance of having a precise evaluation of the positron production for dark matter indirect searches.

Likewise, the study and optimization of the antiproton cross-sections from CR interactions is ongoing. These cross-sections are currently the main source of uncertainty in the modelling of CR antiprotons, as discussed in reference [315] and are very important to improve the accuracy of indirect dark matter searches from CRs. At the moment, there are no detailed analyses of the FLUKA antiproton cross-sections, but a recent preliminary study demonstrated that the flux of antiprotons on Earth predicted using the last evaluation of the FLUKA cross-sections for antiproton production was not compatible with the experimental data from AMS–02, and more than 50% above the most recent cross-sections parameterizations. A dedicated campaign to improve these cross-sections is foreseen for the near future. Following a similar motivation, the collaboration has agreed to implement the anti-nuclei production in the code, because of the potential anti-deuterons and anti-helium to reveal the existence of exotic phenomena in the Galaxy and anti-deuterons are expected to be observed in the next decade.

Finally, the new generation of experiments designed for the study of gamma rays at very high energies (above the TeV) and the lack of accurate evaluations of the gamma-ray production cross-sections at these energies motivates the use of FLUKA for this purpose. Recent studies demonstrate the importance of having refined cross-sections predictions for neutrinos and gamma rays at these energies, in order to characterize better the transport of CRs in different regions of the Galaxy [316, 317]. Therefore, the study of the cross-sections for gamma-ray and neutrino production at high energies is already running.

The Local Interstellar Spectra of primary and secondary CRs based on reference [151] are included in the FLUKA code starting from the release 2023 (they can be activated by a new option in SPECSOUR: “SDUM = GCR-AMS” [318] (see also chapter 16 of the FLUKA Ref. manual [296]).

8.3. Muon collider applications

The renewed interest of the scientific community towards studies for future multi-TeV muon colliders has prompted the application of the NUNDIS engine (Neutrino INteractions in FLUKA) to the simulation of the radiation environment due to neutrino interactions. The collider itself is designed to be operated underground, but neutrinos can travel very long distances, thus the radiation levels at the locations where the radiation resurfaces have to be assessed and compared with the limit value for the effective dose for members of the public, which in most legislations is set at 1 mSv/year.

Neutrinos from muon decays in a muon collider ring are strongly collimated, being emitted within a cone of approx. 1/γ radians around the muon direction (here γ is the muon Lorentz factor). At 1 TeV, 1/γ ≈ 10⁻⁴, resulting in a 10 m spot at a distance of 100 km from the production point. Neutrino interaction products are also forward-oriented. Despite the relative smallness of the neutrino interaction cross-section, the high intensity of the muon beams and the focused conditions may result in a radiation hazard for the population at the point where the neutrino interaction products reach the earth’s surface [319, 320]. The shape of the radiation field from the ensemble of the ring is that of a flat disc, and at large enough distances it does not depend on the circumference of the ring. Conversely, radiation originating from muon decays in the straight sections, for instance between bending magnets or at collision points, is contained in a narrow cone.

Simulations of the radiological hazards of neutrinos have been carried out with FLUKA [322]. For a preliminary assessment, simple descriptions of geometry and muon beams have been implemented. The Earth is assumed to be exactly spherical and its surface flat, without mountains or valleys. The density of the Earth is uniform, the muon beam is perfectly collimated, and the collider is placed at a fixed depth underground.

Neutrino interaction products are fully transported and the ambient dose equivalent (H^*(10)) is calculated online by convolution of particle fluence and conversion coefficients.

Exploiting the relation between the depth of the collider and the distance at which the neutrino beam reaches the surface, this simple setup allows us to calculate the radiation levels both as a function of distance from the collider and as a function of the depth of the collider.

Figure 63 shows a pictorial view of the propagation of radiation from the collider ring to the surface. Even after traveling almost 100 km, the vertical extent of the radiation field is of the order of a few tens of metres. Results are normalized to 1.2 × 10²¹ muon decays per year, which is a reasonable assumption [321].

FLUKA results have been compared with previous analytical approximations [320, 321] and with results [319] from the MARS [323] code. Results from the four approaches compare favorably, as in Figure 64. Also shown is the radiation limit for population.

Radiation levels from straight sections depend on the fraction of muon decays within the straight section, thus on the ratio of its length to that of the full machine. Figure 65, shows the calculated H*(10) from a straight section with length equal to 1/10 000 of the ring circumference, as a function of distance.

It is evident from Figures 64 and 65 that solutions will have to be devised in order to protect the population. Straight sections must be limited, the beam orientation could be varied with time in order to spread the radiation, the site can be chosen in order to optimize the resurfacing location, etc.

Besides radiation levels, FLUKA has been used to evaluate the background and the dose levels in the detectors surrounding the collision points, as described in [324]

8.4. Radiation protection applications

Since FLUKA early application to radiation protection in the high-energy accelerator domain, the code’s capabilities evolved to tackle complex radiation safety challenges in a broader context. Some features that greatly facilitate shielding studies include the various predefined geometry shapes, a library with all elements, and several compounds that can be easily customized by the user to design complex models of shielding structures. Its capacity to handle induced radioactivity, providing information about nuclide inventory, their decay, and the transport of residual radiation, has greatly expanded FLUKA range of applications in the framework of radiation protection.

Another characteristics of the code is its built-in variance reduction techniques that can be easily employed to many deep-shielding problems and, since recently, the inclusion of the already detailed neutron point-wise treatment.

Various radiological quantities are provided as built-in estimators in FLUKA. The default operational quantities for assessing potential effective doses in a given area now correspond to ambient dose H^* coefficients as defined by the joint ICRP116/ICRU95 publications with a few higher energy points extrapolated [325]. Moreover, the code allows the user to select alternative radiation weighting factors, for instance from ICRP74, or even to provide the energy-dependent coefficients for the conversion of fluence to effective dose and ambient dose equivalent [325, 326].

The dose conversion coefficients exist only for protons, neutrons, muons, photons, electrons/positrons and charged pions. The ⁴He conversion coefficients are also available for ICRP116 based quantities only and, for all other particles, the conversion factors were considered to be null until Fluka2023. Since then, an attempt is made to return a non-zero conversion coefficient for heavier ions based on the rescaling of known weighting factors. Although this approach can be deemed as a reasonable effort to avoid neglecting completely the contribution of heavy ions, it has inherent risks that should not be overlooked particularly in dosimetry calculations where heavy ions are relevant.

Dose averaged LET maps are now available as USRBIN estimators, they are computed by multiplying charged particles’ track segments with its restricted LET in water in order to obtain the corresponding water dose deposition, weighting it with LET_∞ in water, and dividing the result by water dose deposition and bin volume.

Recent examples where FLUKA was extensively used in the context of radiation protection include TRIUMF’s IAMI and ARIEL facilities. In the former case, the IAMI shielding design was optimized with several simulations to validate the occupancy of all accessible spaces in the vicinity of a TR–24 proton cyclotron, as well as its target stations’ routine operation for isotope production [327]. In the second case, the entire ARIEL facility was modelled with FLUKA with its shielding layers being constructed in stages [328]. The simulation results particularly guided the shielding design of two target stations which are to be simultaneously irradiated with proton and electron beams for radioactive ion beam production [329]. Additionally, these studies also helped to optimize the target design and to define optimal strategies to mitigate exposure to personnel during target exchange, due to stored targets, radioactive gas and activated coolant [328]. Among the various scenarios studied with FLUKA were also failure scenarios related to beam losses [330] and target failures [331], up to long-term activation of shielding in view of the decommissioning of the facility.

FLUKA is the standard Monte Carlo radiation transport tool of choice at synchrotron light source facilities (see, e.g. [333] and references therein), where it is used to design and validate radiation shielding for accelerator structures, beam lines, and upgrades. An example of a typical calculation to evaluate the radiological impact in the vicinity of a beam line’s first optical enclosure, due to gas bremsstrahlung generated in the NSLS–II synchrotron storage ring, is shown in Figure 66. Another example, shown in Figure 67, illustrates a residual dose rate map resulting from the activated electron beam dump at the MAX IV synchrotron light source.

The capability of FLUKA to work with voxel geometries from computed tomography (CT) scans has also been utilized to investigate new radioprotection perspectives in a radiotherapy setting, estimating organ doses and effective dose for caretakers or others present in the room during pediatric proton therapy [334].

8.4.1. Radiation effects

A wide range of radiation effects in materials can be estimated with FLUKA built-in estimators. These can be used to predict structural damage to materials under irradiation, and thus are very relevant in applications such as nuclear reactors, particle accelerators and space technology, where these effects are of critical importance.

Some of the effects and quantities that are relevant for studying a material behaviour under irradiation with FLUKA include the energy based estimators such as total energy deposition, which can be used to assess temperature rise in a target material, stresses and deformations; energy deposition due to EMF which can be applied to radiolysis calculations and; absorbed dose calculation which can be used to assess material deterioration. Material damage can be also evaluated by scoring in-target gas production, or the Non-Ionizing Energy Losses (NIEL) in a material. Moreover, FLUKA allows the users to assess damage to materials using the displacement per atom estimator. To do so, one only needs to provide a single material-dependent parameter, the damage threshold E_th, for each material of interest. The conversion to DPA is performed via a modified Kinchin-Pease damage model, which accounts for the effect on the displacement efficiency for higher recoil energies due to the recombination and migration of the Frenkel pairs [335]³.

In the context of electronics, FLUKA can be used to estimate both single event effects and cumulative effects. For instance, the rate of single event upset or latch-up can be estimated by scoring the fluence of hadrons with energy above approximately 20 MeV, coupled with high energy hadron or thermal neutron equivalent estimators. On the other hand, for cumulative effects, total ionizing dose can scored, as well as the Silicon 1 MeV Neutron Equivalent estimator to assess displacement damage.

8.5. Nuclear well logging applications

Monte Carlo particle transport simulation software plays a crucial role in the research and development of nuclear well logging tools [336]. We performed a study to investigate new applications of FLUKA in nuclear well logging industry. Instead of FLUKA, MCNP (Monte Carlo N-Particle Transport Code) [217] is widely used as a standard tool due to its accurate simulation results in the industry. However, some countries face accessibility issues due to licensing export restrictions. FLUKA comes into consideration as a potential alternative. Currently, there is a lack of published studies on FLUKA’s application in nuclear well logging [337–339]. Therefore, here we compare performance of FLUKA and MCNP using selected benchmark scenarios. The findings could offer valuable insights into effectively utilizing FLUKA to address scientific problems and challenges in nuclear well-logging.

Natural gamma-ray spectral logging

Model 1 (Fig. 68a) is built to validate the scenario of natural gamma-ray spectral logging. It simplifies the borehole conditions and logging tool structure. The model consists of three main components: formation, borehole, and detector. The formation is composed of 50% sandstone, 40% limestone, and 10% water, the borehole is filled with fresh water, and the detector is a virtual region to accept gamma rays. The source particles are uniformly distributed in the formation and their energies are sampled from the gamma-ray spectra of ²³⁸U decay series, ²³²Th decay series, and ⁴⁰K, with an intensity ratio of 2:3:15 [340].

Figure 69 illustrates the simulated responses of FLUKA and MCNP. It shows that FLUKA and MCNP demonstrate a good agreement. Specifically, when considering the characteristic energy peak corresponding to ⁴⁰K at 1.46 MeV, the relative error is less than 1%. Similarly, for ²³⁸U’s characteristic energy peak at 1.76 MeV, the relative error remains below 3%. In the case of ²³²Th’s characteristic energy peak at 2.62 MeV, however, the relative error increases to 14%. Despite this increase, when considering statistical errors, the results from FLUKA and MCNP are within the error margins.

Neutron porosity logging

Model 2 is established based on the previous Model 1 to verify the neutron porosity logging scenario. As shown in Figure 68b, the model has a shield and a neutron source positioned below the detector. The cylindrical shield is made of cadmium. The neutron point source is located 32 cm below the detector. It emits mono-energetic and isotropic neutrons with an energy of 2.4 MeV to mimic neutrons from a deuterium–deuterium (D–D) neutron generator. Thermal neutrons can be detected by the virtual detector in this arrangement. The geological formation is characterized as sandstone with sequential porosities of 0, 10%, 20%, and 30%, with each pore space filled with freshwater. Figure 70 presents the detected relative thermal neutron counts in strata with varying porosities. The average relative error between FLUKA and MCNP simulation data is within 2%, indicating good agreement between the results obtained from both platforms.

Compensated density logging.

To further assess the applicability of FLUKA in the field of nuclear well logging, a more complex Monte Carlo simulation model (Model 3) is developed. This model is associated with the compensated density logging scenario, illustrated in Figure 68c. This comprehensive model encompasses the complete logging instrument geometry, with the formation lithology set as sandstone and the borehole filled with fresh water. A typical compensated density logging tool utilizes ¹³⁷Cs as a gamma-ray source with an energy of 0.662 MeV. Additionally, it consists of both near and far detectors for gamma-ray detection. The tool is placed against the borehole wall.

Different simulations are conducted to obtain the response of the compensated density logging tool under different formation densities. We vary the porosity of the sandstone formation from 0 to 95% in 5% increments. Figures 71 and 72 plot the logarithm of relative counts against formation density for both far and near detectors. The curves displayed in the figures exhibit consistent trends, suggesting a robust correlation between the variables. For the near detector, the root mean squared error (RMSE) of the logarithm of relative counts from FLUKA and MCNP simulations is within 0.6%. Similarly, the far detector exhibits an RMSE within 1.4%. The data further emphasize the reliability of the results.

In summary, we applied FLUKA in nuclear well logging scenarios, addressing three different logging methods and we compared the obtained results to MCNP.

Models 1 and 2 are designed to simulate natural gamma-ray spectral logging and neutron porosity logging scenarios, respectively. A comparative analysis of FLUKA and MCNP shows a very good agreement between the two generated gamma spectra and thermal neutron counts. When considering statistical errors, no significant differences are observed between the two.

Model 3 focuses on the modeling of compensated density logging, incorporating a comprehensive description of the logging instrument’s geometry. The study analyzes the relationship between count rate and formation density for both far and near detectors, and identifies consistent trends in the curves. The results indicate that there are discrepancies between FLUKA and MCNP in terms of count rates that fall within the range of statistical errors.

In conclusion, FLUKA demonstrates excellent performance in modeling and simulating nuclear well logging problems, serving as a reliable candidate to replace MCNP. FLUKA accurately captures the responses of logging instruments to various logging signals and accurately simulates complex geometries and physical processes. However, further exploration of nuclear well logging scenarios and experimental verification using field-measured data are necessary to thoroughly assess the applicability and accuracy of FLUKA in this field.

9. Technical improvements

As well as improvements to the physics models in FLUKA, the new version also includes a number of technical improvements:

• Dynamic allocation of memory allows for “unlimited” number of regions, as well as an automatic increase of the memory size in case the scoring configuration or voxel geometry should require it;

• Up to two nested roto-translations are now allowed in geometry, both for bodies and lattices;

• The possibility of define a spatially distributed source from a previously created USRBIN file;

• Built-in source routines for neutron spectra from Am-Be, Am-B, ²⁵²Cf, deuterium-deuterium, and deuterium-tritium neutron sources are now available;

• The interface to the CORSIKA7 [341] air shower simulation program has been completely reworked and made modular. It now optionally supports, besides non-elastic and (new) quasi-elastic, also elastic hadron–nucleus interactions, photo-nuclear interactions, electro-nuclear interactions, and nucleus-nucleus Electro-Magnetic Dissociation events. The interface is supporting also the new CORSIKA8 [342] air shower simulation code;

• New geometrical bodies, pyramids along X-, Y- or Z-axis and generic tetrahedra have been added to the geometry;

• The rQMD event initialization has been improved and made significantly faster.

10. Conclusions

The modern FLUKA has a long history dating back 35 years. During these years, its models and their applications have evolved to cover a wide range of problems in nuclear physics, high energy physics, cosmic ray physics and medical physics, far beyond the traditional field of radiation protection. The fourth generation of FLUKA has more than 6000 registered users worldwide and the number is growing every day.

This paper illustrates roughly the last decade of FLUKA developments with examples of old and new applications made possible by FLUKA capabilities.

Along this period of time, the quality of physics models of FLUKA has been constantly improved, trying to keep up with the progress in the comprehension of the different processes of interest, the appearance of new experimental data and the general developments achieved by the scientific community. This line of development will of course continue also in the future, and, among the many different subjects of interest, we believe that the most important directions concern the following items:

1. High-energy hadronic interactions: following the now widely accepted paradigm that universal hadronization is no longer a viable concept and that new mechanisms are therefore required (see Sect. 4.3).

2. Further refinements and developments are also underway for photonuclear interactions and, very importantly for the next generation of experiments, for neutrino interactions.

3. The relatively recent development of fully correlated point-wise neutron interactions has opened up many new fields for FLUKA, and further refinements are expected in the coming years.

Thanks to a highly motivated group of collaborators from various fields who have contributed to past achievements and will contribute to future ones, the pace of development is now even faster than in the past.

---

Acknowledgments

The authors thank Dr. Mattias Lantz for the help in carefully proofreading the manuscript.

Funding

This work was partially supported by Department of Energy contract DE-AC02-76SF00515. We wish to thank INFN for the support received from 2002 to 2023 through the CSN5 projects FLUKA, FLUKA2 and MC-INFN. PDL was supported by the European Research Council under grant 742104 and the Swedish National Space Agency under contract 117/19, at the first stages of this project. PDL is currently supported by the Juan de la Cierva JDC2022-048916-I grant, funded by MCIU/AEI/10.13039/501100011033 European Union “NextGenerationEU”/PRTR. The work of PDL is also supported by the grants PID2021-125331NB-I00 and CEX2020-001007-S, both funded by MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”. PDL also acknowledges the MultiDark Network, ref. RED2022-134411-T.

Conflicts of interest

The authors declare that they have no competing interests to report.

Data availability statement

This article has no associated data generated and/or analyzed.

Author contribution statement

All authors contributed at least to one of the developments and applications presented in this paper. All authors participated to the validation of FLUKA. More in detail: F. Ballarini, G. Battistoni, M.G. Bisogni, T.T. Böhlen, M.P. Carante, A. De Gregorio, R. Dos Santos Augusto, A. Ferrari, G. Franciosini, A. Kraan, J. Lascaud, G. Magro, A. Mairani, I. Mattei, M.C. Morone, S. Muraro, P.G. Ortega, K. Parodi, V. Patera, R,L. Ramos, V. Rosso, T. Tessonnier, K. Ytre-Hauge contributed to the medical, radiobiological and PET developments and applications of FLUKA presented in this paper. K. Batkov, D. Chen, A. Fassò, A. Ferrari, W. Li, J. Liu, Z. Liu, M.C. Morone, P.R. Sala contributed to the low energy neutron models and related applications. G. Battistoni, P. De la Torre Luque, R. Engel, A. Fedynitch, A. Ferrari, M.N. Mazziotta, S. Muraro, L. Pinsky, P.R. Sala, J. Ranft contributed to the space and cosmic ray related models and applications. T.T. Böhlen, R. Dos Santos Augusto, A. Fassò, A. Ferrari, G. Magro, A. Mairani, K. Parodi, V. Patera, P.R. Sala contributed to the FLUKA atomic models and their applications. G. Battistoni, M. Campanella, A. Ferrari, S. Müller, P.R. Sala contributed to the FLUKA software architecture. D. Chen, W. Li, J. Liu, Z. Liu contributed to the well logging application. K. Batkov, P. Degtiarenko, R. Dos Santos Augusto, A. Fassò, A. Ferrari, M. Santana Leitner, L. Zana contributed to the development or validation of the FLUKA photonuclear models. K. Batkov, R. Dos Santos Augusto, A. Fassò, A. Ferrari, A. Ferrari, S. Müller, M. Santana Leitner contributed to the radioprotection and shielding applications and associated developments. R. Engel, A. Fedynitch, J. Ranft contributed to the development and validation of DPMJET-3. G. Battistoni, P. De la Torre Luque, R. Engel, A. Ferrari, M.N. Mazziotta, J. Ranft, P.R. Sala contributed to the development and validation of the FLUKA nuclear interaction models. R. Engel, A. Fedynitch, A. Ferrari, A. Mairani, L. Pinsky, J. Ranft, P.R. Sala contributed to the development and validation of the FLUKA nucleus-nucleus interaction models. G. Battistoni, A. Ferrari, A. Ferrari, P.R. Sala contributed to the development of the neutrino interaction models and their application to a possible muon collider. A. Ferrari, M. Santana Leitner contributed to the development and validation of the new FLUKA electric and magnetic field transport algorithm.

References

All Tables

All Figures