Extending deterministic transport capabilities for very-high and ultra-high energy electron beams

Focused Very-High Energy Electron (VHEE, 50–300 MeV) and Ultra-High Energy Electron (UHEE, > 300 MeV) beams can accurately target both large and deeply seated human tumors with high sparing properties, while avoiding the spatial requirements and cost of proton and heavy ion facilities. Advanced testing phases are underway at the CLEAR facilities at CERN (Switzerland), NLCTA at Stanford (USA), and SPARC at INFN (Italy), aiming to accelerate the transition to clinical application. Currently, Monte Carlo (MC) transport is the sole paradigm supporting preclinical trials and imminent clinical deployment. In this paper, we propose an alternative: the first extension of the nuclear-reactor deterministic chain Njoy-Dragon for VHEE and UHEE applications. We have extended the Boltzmann-Fokker-Planck (BFP) multigroup formalism and validated it using standard radio-oncology benchmarks, complex assemblies with a wide range of atomic numbers, and comprehensive irradiation of the entire periodic table. We report that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$99\%$$\end{document}99% of water voxels exhibit a BFP-MC deviation below \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\%$$\end{document}2% for electron energies under \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text{1.5 GeV}$$\end{document}1.5 GeV. Additionally, we demonstrate that at least \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$97\%$$\end{document}97% of voxels of bone, lung, adipose tissue, muscle, soft tissue, tumor, steel, and aluminum meet the same criterion between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text{50 MeV}$$\end{document}50 MeV and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text{1.5 GeV}$$\end{document}1.5 GeV. For water, the thorax, and the breast intra-operative benchmark, typical average BFP-MC deviations of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.3\%$$\end{document}0.3% and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.4\%$$\end{document}0.4% were observed at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text{300 MeV}$$\end{document}300 MeV and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text{1 GeV}$$\end{document}1 GeV, respectively. By irradiating the entire periodic table, we observed similar performance between lithium (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z=3$$\end{document}Z=3) and cerium (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z=58$$\end{document}Z=58). Deficiencies observed between praseodymium (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z=59$$\end{document}Z=59) and einsteinium (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z=99$$\end{document}Z=99) have been reported, analyzed, and quantified, offering critical insights for the ongoing development of the Evaluated Nuclear Data File mode in Njoy.

1 GeV , respectively.By irradiating the entire periodic table, we observed similar performance between lithium ( Z = 3 ) and cerium ( Z = 58 ).Deficiencies observed between praseodymium ( Z = 59 ) and einsteinium ( Z = 99 ) have been reported, analyzed, and quantified, offering critical insights for the ongoing development of the Evaluated Nuclear Data File mode in Njoy.
More than half of the 19 million cancer patients, diagnosed worldwide each year, receive radiotherapy (RT) treatment during the course of their disease 1,2 .Survival rates differ starkly between cancer types and conventional RT (CONV-RT) techniques 3,4 .Schulz and Kagan 5 outlined that CONV-RT 5-year survival rate for endometrium, breast, bladder and colorectal cancers is 87.9% , 74.6% , 72.6% and 49.8% , respectively, while it is only 2.5% , 3.8% , 4.6% and 10% for pancreas, liver, esophagus and lung cancers.
Cure and mortality patterns depend on Holthusen's therapeutic window (TW) width 6 ; i.e., the difference between the tumor control probability (TCP) and the normal tissue complication probability (NTCP) 7,8 .Oncologists acknowledge that RT survival rate can be improved if (1) current tumor control is guaranteed 9 ; (2) the risk of developing fatal secondary malignancies-draining regional lymphatics and lymph nodes-is minimized or preserved 10 ; (3) with a marked reduction in the collateral damage inflicted on healthy tissues 11 .

Introduction to NJOY workflow
Proposed in 1973 as a natural successor to the Multigroup Interpretation of Nuclear X-sections (MINX) 107 , the development of Njoy was entirely funded by the U.S. Fast Breeder Reactor 108 and Weapons Programs 106 .Njoy benefited from the merging of several algorithms and codes, namely ETOPL kernels for Pointwise Evaluated Nuclear Data Files (PENDF) libraries 109 , RESEND for union grid with resonance reconstruction 110 , SIGMA for Doppler-broadening 111 , ETOX for unresolved resonance self-shielding 112 , LAPHANO0 for photonic production 113 , CAMLEG for photonic interaction 114 , FLANGE-II 115 and HEXSCAT 116 for thermal neutron scattering, TRANSX for library formatting 105 , SAMMY for Reich-Moore-Limited resonance representation 117 and MACK for heat production and radiation damage energy production 118 .In 1977 and 1979, Njoy garnered support from the U.S. Electric Power Research Institute (EPRI) and the U.S. Magnetic Fusion Energy Program for the development of Groupwise Evaluated Nuclear Data Files (GENDF) libraries, enabling data formatting compatible with EPRI neutronic solvers and covariance production, respectively.Between 1981 and 2016, Njoy underwent numerous improvements, building upon previous work and leveraging the accumulated experience of numerous international contributors.This collaborative effort resulted in enhanced stability, code maturity, an open-source release in 2016, and a versatile capability for processing various types of data and formats.These include the US.ENDF/B 119,120 , the JEFF libraries in Europe 121 , JENDL in Japan 122 , TENDL 123 , CENDL in China 124 , BROND 125 and RUSFOND 126 in Russia, and the specialized libraries of the International Atomic Energy Agency (IAEA) Nuclear Data Section (NDS) [127][128][129] .Njoy has earned its distinctiveness and near-domination over the years through its enduring capacity to concurrently handle all these evaluations and their updates.A single problematic isotope can obstruct the system's operation, underscoring the necessity of these qualities in the processing of certified nuclear data for nuclear reactor physics.For over 51 years, the Njoy nuclear data processing system (Fig. 1) has Figure 1.Njoy nuclear data processing system simplified workflow.RECONR resconstructs PENDF cross sections from ENDF; BROADR computes Doppler-broadened cross sections; UNRESR computes self-shielded cross sections in the unresolved range; HEATR produces heat KERMA and radiation damage cross sections; THERMR computes free and bound scatters thermal cross sections; GROUPR applies multigroup theory for neutron cross sections; GAMINR applies multigroup theory for photon cross sections; ELECTR applies multigroup theory for electron cross sections; PURR prepares unresolved-region probability tables for MC neutron transport; GASPR generates gas-production cross sections; ERROR computes multigroup covariance matrices; COVR performs covariance plotting; MATXSR, WIMSR, ACER, DTFR, POWR, CCCR and RESXSR format multigroup and pointwise data for specialized codes and applications.played a pivotal role in supporting certified and qualified Boltzmann and MC codes' missions, including neutral particle transport, regulatory licensing processes, advanced fission and fusion reactor design, safety assessment, criticality safety benchmarking, stockpile stewardship modeling, radiation shielding, and nuclear waste management.Throughout these years, Njoy has remained limited to neutron and photon transport.Our introduction of ELECTR (Fig. 1) represents the first expansion of the system's functionalities to accommodate light charged particles 130,131 .ELECTR operates under two distinct modes: ENDF 120 and CEPXS 132 .This paper concerns the extension of CEPXS-mode capabilities for VHEE and UHEE beams.The multigroup CEPXS formalism from the U.S. Sandia National Laboratory has never been validated beyond 20 MeV 104 .The objective of this paper is to demonstrate that an extension beyond its design limit of 100 MeV , announced at its release in October 1989 132 , is possible.It is important to place this effort into context.The CEPXS-mode in ELECTR serves as a foundation for the development of the ENDF-mode.The Lorence-Morel-Valdez postulates 132 , that we will discuss in what follows, provide a basis for a first application of multigroup theory to ENDF data.For this reason, we propose this final step: to extend the state-of-the-art to VHEE and UHEE beams prior to any release of a fully operational open-source ENDF-mode.Unlike the last work proposed for the CONV-RT 104 , we are now proposing a detailed presentation of the analytical cross sections implemented in the CEPXS mode in VHEE-and UHEE-RT.

VHEE multigroup theory
Let g, n and r be the VHEE/UHEE energy group, its discrete ordinate direction and position, respectively.If ψ e (e − /cm 2 /s) designates the electron flux, t (cm −1 ) the macroscopic total cross section, x l (cm −1 ) the lth Legendre macroscopic scattering cross section coefficient associated to interaction x and Q e,ext.(e − /cm 3 /s) the external source, the S N multigroup form of the first-order BFP equation is given by: where R m l denote the real spherical harmonics.The electron direction ( ˆ n ) is uniquely defined by the polar angle cosine ( µ n ∈ [−1, +1] ) and the azimuthal angle ( φ n ∈ [0, 2π] ).Total and scattering cross sections in Eq. ( 1) are restricted to catastrophic collisions while stopping powers ( β g [MeV/cm]) and momentum transfers ( α g [cm −1 ]) are restricted to soft ones.The mth moment of the lth Legendre order group-to-group microscopic transfer cross section ( barns • MeV m ) is given by: Two quantities are needed for each x, the DDSC ( σ x in Eq. 2) and the VHEE/UHEE transfer groups ( I x g in Eq. 2).The ELECTR [CEPXS-mode] considers the ionization ( x = i ) as a two-body mechanism, made on a free electron.No binding energy is involved in the process.Two particles emerge from a typical CEPXS i-event; the scattered primary electron and the delta ( δ ) ray.The i-DDSC is that of Møller.Unlike Class-I and-II MC codes, no straggling model is needed for the Njoy-Dragon chain.The energy distributions are derived directly from Møller's DDSC.The scattering angles for the primary and δ electrons are determined from the i-collision kinemat- ics, with µ i e/δ = e e/δ (e ′ e + 2)/e ′ e (e e/δ + 2) .e e/δ energies are in reduced units.One can verify that both µ i e and µ i δ are forward peaked.For this reason, some MC codes assume that µ i e = 1 and account for the primary angular deflection by a modification to the nuclear elastic kernel.This is not the case for the Njoy-Dragon chain.The same kinematics provide e δ ∈ [0, e ′ e /2] .However, because of (1) the singularity of Møller's DDSC (Eqs. 3-5)at e δ = 0 ; and ( 2 ] .Integral boundaries for multigroup transfer matrices can be deduced easily by crossing the Njoy[ELECTR] user group structure with the kinematics restrictions. (1) g,n ( r).
(2)   www.nature.com/scientificreports/where β ′ e refers to the incident electron's velocity in the speed of light units and r 0 to the classical electron radius.H is the Heaviside function.From Eqs. ( 4)- (5), it is clear that the degree of correlation between the δ and primary emission spectra is dictated by the user multigroup structure in Njoy.Moreover, the minimum δ-ray energy cannot fall below the Njoy user-selected lower energy limit.Several authors, including Morel 133 , Lorence 132 , Olbrant 134 , Cullen 135-137 and Salvat 138 , recommend a 1 keV transport cutoff.This is also the default value used by Njoy-Dragon for CONV-RT 104 , VHEE-RT and UHEE-RT.A 100 eV limit is possible only in the ELECTR [ENDF-mode].Otherwise, in the CEPXS-mode, if e δ < 1 keV , only the primary electron undergoes deflection, and the δ-ray is not produced.In this case, the associated soft energy loss is addressed by the Fokker-Planck continuous slowing down (CSD) operator (Eq.1).
The catastrophic production of VHEE/UHEE bremsstrahlung electrons ( x = b ) is possible in an electric atomic field or a nuclear one.Electrons are emitted in the incident particle's direction ( µ b e = 1 ), while photons follow a Sommerfield angular distribution 139 .The b-DDSC is that of Berger-Seltzer 140 based on Koch-Motz developments 141,142 and Born's assemblies.Like Møller, Berger-Seltzer's DDSC shows a singularity at e γ = 0 , so the same cutoff e c must be applied for the catastrophic radiative emission.Moreover, a high-frequency limit ( e m γ ) must be defined to avoid the divergence of the atomic-dependent Elwert screening factor f e when e γ = e ′ e .We have: σ 1 and σ 2 refer to Sauter-Gluckstern-Hull 143,144 unscreened and screened scattering kernels, respectively.σ 3 refers to Schiff 145 unrelativistic kernel, while σ 4 to Olsen-Maximon's 146 Coulomb correction for small angles.The first three kernels are derived in the Born approximation, while σ 4 goes beyond that to include a Sommerfeld-Mauve wave function.The DSCs assembly (Eq.6) is made in such a way to avoid discontinuities arising from abrupt switching from one kernel to another.ξ r represents the Berger-Seltzer atomic-dependent correction factor.ω(e ′ e ) is a weighting function introduced by Berger and Seltzer 140 to switch off Coulomb correction at low energies where it becomes unreliable.e t refers to the total incident electron energy.e s denotes the total energy lost during the b-event.p is the parti- cle's momentum in reduced units, L = 2 log[(e t e s + p e p ′ e − 1)/e γ ] and ξ ′ e = log[(e s + p ′ e )/(e s − p ′ e )] . 1 and 2 are the screening factors.ζ is a b-spectrum dependent factor.f(Z) is a Z-dependent empirical correction factor.The nuclear contribution to bremsstrahlung is taken into account in the CEPXS-mode by modifying Z(Z + 1) by Z 2 .One should note that all scattering kernels (Eqs.7-10) vanish in the high-frequency limit.For the latter, the CEPXS-mode uses Berger-Seltzer extrapolations for consistency between the theoretical work of Fano-Koch-Motz 147 , McVoy-Fano 148 , Jarbur-Pratt 149 and the experimental data.
Elastic scattering ( x = e ) is predominantly treated in a nuclear field.The electron changes direction without any loss of energy.In the non-relativistic domain ( E ′ e ≤ 256 keV ), the CEPXS-mode implements Riley's DDSC (5)  154 for medium and high Z and tabulated by Birkhoff and Sherman 155 for discrete scattering angles and incidence energies.ν e is an empirical energy-dependent correction factor.Unlike the ENDF-mode, the CEPXS-mode in ELECTR avoids using angular quadrature for Eq.11 integration.A Goudsmit-Saunderson (GS) 156,157 moment-based semi-analytical approach is used for Eq. ( 11).The Mott-based GS moments are given by: A high-order expansion of G l is required for forward-peaked multiple scattering (e.g., L = 60 for µ e = 0.986286 (9.5 • ) 158 ).To avoid the slowly converging Legendre series in the deflection angle, Spencer's numerical substitution is introduced 159 : This representation has been dem- onstrated to be accurate to within 1% with J = 5 .The latter is referred to as the five Bartlett-Watson tabulated cosine angles 154 and is utilized in the CEPXS-mode.Combining Eqs.11-13 and �(µ e ) , G where recursion forms can be derived from the Legendre polynomials as proposed by Spencer 159 : It is worth noting that the accuracy of Eqs. ( 11)-( 16) remains less than that achieved by Brown's (1961) method 158 , which involves a direct numerical resolution of the Dirac equation with a screened Coulomb potential.Moreover, the Elsepa kernel 161,162 , designed for relativistic Dirac partial-wave elastic scattering and implemented within the Penelope framework 138 , is expected to yield further improved accuracy.Conversely, the Riley DDSC can be derived from a partial-wave expanded Dirac equation resolution in a Poisson spherical central atomic static potential field 163 .This exact solution can be fitted using 12g-dependent parameters 164 .
where the first sum is developed to fit the small-angle deflections and the second one to fit moderate and largeangle tail.B is a screening parameter and A m and C n are in Å/Sr .The idea of the fit is in response to the need in both MC and BFP for an accurate, rapid, and analytic scattering kernel at low energies.Equation ( 17) is validated against Fink and Kessler 165 absolute small-angle measurements and Ibres-Vainshtein 166 Born approximation results.The Riley-based GS moments and Spencer associated functions are given by 159 : Finally, both Riley and Mott's DDSCs follow a transport correction to bring back the highly-forward peaked scattering kernels reducible to a lower-order Legendre expansion 167 .( 11) dµ e P l (µ e ) e e + 1 e e (e e + 2) � e e + 1 e e (e e + 2) www.nature.com/scientificreports/Unlike the ENDF-mode, the Auger ( x = a ) and fluorescence ( x = f ) relaxation cascades in the CEPXS-mode are uncorrelated with i-events.Consequently, an impact ionization cross-section is separately introduced in the CEPXS-mode.All a-and f-events involving secondary emission below the Njoy-Dragon chain cutoff ( e c ) are excluded.As a result, only K, L 1 , L 2 , L 3 , M, and N shells and subshells, involving Z > 10 , Z > 27 , Z > 29 , Z > 29 , Z > 51 , and Z > 84 elements, are considered as cascade candidates.The energy of the emitted particle is determined by the difference in binding energies of the involved subshells.The Auger electron emission is isotropic ( l = 0 ), and thus, the Auger transfer matrix is given by: where j refers to the ionized shell and k to the allowed line radiation.The maximum number of subshells ( N s ) and transition lines ( N t ) are restricted to 5 and 28, respectively.υ j e denotes the number of electrons in the jth subshell.η a jk indicates the relaxation efficiency for a kth Auger emission following the jth subshell ionization event.The g k transfer group arises from loop convergence, assigning the emitted Auger electron to one of the Njoy user-selected groups.Given that e j b represents the jth subshell binding energy in reduced units, the Gryzinski 168 impact ionization cross section is expressed as: where β j = e All N shells intervene with e N b = 0 .The relaxation process, triggered by an inner vacancy, ceases if the latter is transferred to the outermost shell.Each relaxation cascade efficiency (radiative or non-radiative) is computed from a multiplication of the specific conditional branch probabilities.
For a given incidence group, the total catastrophic cross-section (Eq. 1) is obtained by first integrating all scattering kernels (Eqs.4, 6-10, 12-16, and 17-19) over all Legendre orders and all possible emission energies, and then summing the resultant values.This second integral must adhere to the definition of the catastrophic jump, the Møller non-divergence condition, the associated kinematic restrictions, and the bremsstrahlung highfrequency limit.
In ELECTR [CEPXS-mode], soft VHEE/UHEE stopping powers ( β g ) for both radiative and non-radiative events are derived neither from single-event cross sections nor from down-scattering in adjacent energy groups.Instead, Berger's collisional ( β c t ) and radiative ( β r t ) total stopping power tabulations are extrapolated at the incident VHEE/UHEE midpoint group ( e m g ′ ) and corrected for the plasmon effect using a Lorence correction.Soft stopping powers are obtained by subtracting catastrophic stopping powers, derived by integrating the first moment of group-to-group transfer cross sections across all Legendre orders and energy losses, from β c t and β r t .Soft collisions do not involve energy loss straggling or angular deflection.
Soft collisions as well as catastrophic i-, b-and a-events define the energy deposition cross section.We can infer catastrophic energy deposits from local absorptions.Therefore, the total energy deposition cross section can be expressed as follows: Once the MATXS library is produced, the numerical solution of Eq. 1 can be found using the Dragon-5 solver.First, we assume a continuous variation of the flux over each group and spatial domain and a linear variation of the angular flux and scattering source over each group.Next, by (1) applying the Galerkin method of weighted residuals; (2) substituting the flux and scattering source with a second-order energy expansion; (3) multiplying Eq. (1) (in 1D) by normalized Legendre polynomials; (4) integrating over the appropriate support E ∈ [−1/2, +1/2] ; and (5) canceling the first moments of the flux and scattering sources, it can be shown that: with the energy propagation relation given by: where β ∓ g ( r) and ψ e,∓ g,n ( r) refer to the soft stopping powers and angular fluxes at the gth group upper and lower boundaries, respectively.The base cosine angles ( µ n ) and associated weights ( w n ) are given by an N-point Gauss-Legendre quadrature.Prior to inverting the system, spatial discretization is performed.Denoting ψ e,i ± g,n and ψ e,i g,n as the electron flux's mesh-edge and mesh-centered values for a specific sub-mesh i, the same procedure applied in energy can be employed.By (1) expanding the flux and sources' spatial expansions of order M using normalized Legendre polynomials; (2) weighting Eqs. ( 27) and (28) with normalized (in space) Legendre moments; (3) integrating Eqs. ( 27) and ( 28) over each submesh region to obtain the spatial moments of Eq. ( 27) and spatial Legendre moments of the slowing down angular flux (Eq.28); and (4) applying a diamond-difference discretization of the electron flux: With the variable change Bienvenue and Hébert 169 showed that the discre- tized 1D BFP equation can be resolved recursively by initiating with known entering fluxes (generally a vacuum boundary condition) and computing the other mesh-edge and flux moments: ( (29) g,n,i + 2 √ 5ψ [2] g,n,i − ψ e,i ± g,n ∀M = 2 ...  (36).Utilizing these values, the subsequent iteration's spatial and angular Legendre moments, L g,n,i , are calculated 169 .

Results
The 1D, 2D and 3D high-order diamond difference schemes as well as the classical, linear and quadratic discontinuous Galerkin schemes are already implemented in the Dragon-5 solver [169][170][171] .This study aims to extend and detect any anomalies in the VHEE/UHEE multigroup formalism in the CEPXS-mode.Since VHEE and UHEE interact over much larger ranges than CONV-RT, larger spatial dimensions are necessary for complete beam attenuation.To accurately evaluate the Njoy-Dragon chain's multigroup performance for VHEE/UHEE-RT, it is crucial to minimize interference with the BFP equation's numerical resolution and avoid error compensation effects related to ray effects or the S N discretization of the 3D Boltzmann catastrophic kernel.Consequently, deterministic transport is limited to one spatial dimension.Pure numerical investigations for higher dimensions can be found elsewhere 171 .Njoy-Dragon computational schemes for VHEE/UHEE consist of N g = 300 groups, an S 16 Gauss-Legendre quadrature, a P 15 Legendre anisotropy, and D r ≥ 100 voxels.A convergence criterion of 10 −5 was imposed on the internal iterations of the electron flux.VHEE medium self-polarization is corrected by the classical Sternheimer-Peierls density effect correction 172,173 .Geant-4, which has been qualified against experiments for VHEE and UHEE beams 66,[68][69][70][71][72]74 , is the reference MC code used for validation purposes. A tatistical uncertainty of 0.2% or better is achieved in each voxel with the G4EmLivermore phys- ics list constructor.Other legacy constructors, such as G4EmPenelope and GEmStandard 174 , were also used to verify the stability of Geant-4 response at VHEE and UHEE.The classical Lewis condensed-history (CH) theory 175 , implemented in Urbán's algorithm ( G4UrbanMscModel) 176 , is used below 100 MeV while the Goudsmit-Saunderson condensed-history theory 156,157 , implemented in the advanced Bagulya algorithm ( G4GoudsmitSaundersonMscModel) 177 , is used for higher energies.The default G4MultipleScat- tering step size is used for MC soft collisions.Fluorescence and bremsstrahlung photons are produced in Geant-4 and Dragon-5 and immediately eliminated at the point of birth.Our aim is to retain purely electron transport, eliminating any coupling effect or dose deposition of a photonic nature.The resulting dose is therefore strictly electronic.This strategy ensures the BFP-MC comparison remains unaffected by this operation, owing to the equivalence of photonic elimination processes in both Geant-4 and Njoy-Dragon-5.Both codes use a 1 keV secondary production threshold and transport cutoff in compliance with the CEPXS-mode requirements in the ELECTR module.

Radiation oncology benchmarks
Figure 3a-d shows the dose profiles for four ascending transport complexity benchmarks 104 .The dimensions of each benchmark are determined iteratively based on the maximum range of the incident electron beam.The investigated beams range from 1 MeV to 6 GeV , with a primary focus on the VHEE and UHEE domains.Insets display the relative deviations ( ǫ r ), expressed as a percentage, between the Njoy-Dragon chain and Geant-4.
The error is computed using a fine piecewise cubic Hermite grid unification applied to both BFP and MC dose detectors.We have optimized the irradiation process for full beam attenuation to enhance electron transport challenges in this study, deliberately not conforming to typical clinical settings with smaller dimensions and, thus, less demanding buildup, backscattering, and attenuation scenarios.In Fig. 4a-d, the percentage of voxels within the Njoy-Dragon chain that satisfy a deviation criteria of less than 1% ( ǫ 1 ) and 2% ( ǫ 2 ) compared to Geant-4 is quantified for the same benchmarks.The ǫ 2 criterion is deemed, within this context, the optimal instrument for assessing the performance of the BFP response.This premise is rooted in the standards set by the American Association of Physicists in Medicine (AAPM) for satisfactory dose prediction in CONV-RT 178,179 .The ǫ 1 criterion, on the other hand, is seldom mentioned in existing literature.Its introduction in this study serves to underscore specific vulnerabilities.However, it should be perceived as an overly cautious exaggeration as the stated uncertainty on the cross-sections exceeds 1%.
All benchmark samples undergo irradiation with incremental energy levels as follows: increments of 1 MeV in the range of 1 to 20 MeV , 5 MeV between 20 and 100 MeV , and 20 MeV in the 100 MeV to 1 GeV interval.Beyond 1 GeV , the increment increases to 500 MeV .The selection of beams, as presented in Fig. 3a-d, has been made deliberately to enhance visual clarity and minimize overlap.However, the performance of the entire set of considered beams is comprehensively presented in Fig. 4a-d.We truncated the graphical performance display beyond 2 GeV as the extrapolation routines for the CEPXS-mode fail past this energy threshold.
Figures 3a and 4a depict the fundamental case of a spatially homogeneous water (W) slab.The high BFP-MC agreement seen in Fig. 3a leads to 99% of the W-voxels satisfying the AAPM ǫ 2 criterion (Fig. 4a) below 1.5 GeV .Fig. 4a can also be compared to our previous result in CONV-RT 104 , wherein 100% of voxels satisfied the ǫ 2 criterion.The decrease in ǫ 2 by 1 voxel beyond CONV-RT is attributed to our newly introduced deterministic computation scheme, which required a Legendre order of P 15 (in VHEE and UHEE) compared to P 8 and P 12 (in CONV).This necessitated a commensurate reduction in the S N order ( S 16 in VHEE and UHEE vs. S 64 in CONV) and a five-fold decrease in spatial discretization.A Gauss-Legendre quadrature with N = l + 1 points is used here to ensure that the highly forward-peaked scattering is correctly integrated by the Boltzmann kernel.This forms what Morel defined as a Galerkin quadrature 180 , a discrete ordinates method that exactly integrates scattering represented by a delta function.When applying a more conservative ǫ 1 criterion, the percentage of voxels decreases to 70.8-84.4% between 1 and 20 MeV , 91.0-97.9% between 25 and 70 MeV , 98.0-98.7% between 70 and 600 MeV , and 97.9-97.7% between 601 MeV and 1.5 GeV .Figure 4a displays a monotonic decrease in the conformity of BFP-MC above 1.5 GeV .The adherence to ǫ 2 ( ǫ 1 ) reaches 91.8% ( 38.3% ) at 2 GeV , subsequently decreases to 78.0% ( 18.7% ) at 3 GeV , and further diminishes to 53.3% ( 12.6% ) at 5 GeV .Optimization attempts beyond 2 GeV didn't improve the downward trend in ǫ 2 .Thus, we deduce that ELECTR's extrapolation routines are valid up to 1.5 GeV , aligning with the interest limit in medical physics, negating the need for further develop- ment.The ǫ 1 criterion emphasizes (Fig. 4a) the superior performance in VHEE and UHEE domains compared to the CONV range and beyond CONV.Figure 4a highlights that, for water, the ǫ 1 criterion is affected during the 1.0 − 1.5 GeV transition without affecting the compliance with ǫ 2 .Taking all the W-voxels into account, the average absolute difference ( ǭ ) between BFP and MC is 0.19% , 0.23% , 0.30% , and 0.38% at 100 MeV , 300 MeV , 500 MeV , and 1 GeV , respectively.It is noteworthy that this study remains fundamentally a proof-of-concept demonstrating the possibility of transcending the 100 MeV design limit of the CEPXS-mode.The qualification of Dragon-5 for a particular clinical routine will initiate iteration calculations on N g , P l , S N and D r deterministic parameters.Consequently, greater accuracy than those presented and tailored computation times could potentially be reported.Such parameter studies, as in nuclear reactor physics, are typically conducted and constructed on a case-by-case basis.
Figures 3b and 4b refer to the first level of heterogeneity with the thorax (T) benchmark.The thicknesses of the 4 slabs (muscle [ 13% ], bone [ 7% ], lung [ 22% ] and soft tissue [ 58% ]) are deduced from the maximum range of the incident beam's energy.The BFP-MC T-dose profile agreement (Fig. 3b) is maintained regardless of the biological tissue's nature, density, or thickness.The slight decrease in BFP-MC W-compliance concerning the AAPM ǫ 2 criterion (Figs.4b vs. 4a) can be attributed to the Geant-4 boundary crossing effects at the bone inter- face.This effect pertains solely to one voxel and appears only on one interface, the bone-lung interface.Its slight propagation in Fig. 4b is attributed to its spread to artificial voxels following a grid unification operation based on Hermite polynomials for BFP-MC comparison.As observed in Fig. 4b, the percentage of T-voxels satisfying the ǫ 2 criterion reaches 99.2% for energies between 1 and 20 MeV .It then decreases and plateaus at around 97.0% for energies up to 1.5 GeV .Like the W-benchmark, the ǫ 1 criterion in Fig. 4b emphasizes the BFP-MC compliance improvement at the 20 MeV threshold, i.e., the CONV-VHEE transition.The BFP-MC ǫ 1 conformity increases from 65.8 to 93.6% in CONV range, stabilizes around 97.0% for energies between 25 and 200 MeV , and then undergoes a monotonic decline from 97.0 to 80.1% within the range of 200 MeV to 1.5 GeV .We report, however, that ǭ is 0.34% , 0.37% , 0.39% and 0.44% , respectively, at 100 MeV , 300 MeV , 500 MeV and 1 GeV.
Figures 3c and 4c introduce a second level of heterogeneity with the IORT Mobetron benchmark.The thicknesses of the tumor [ 40% ], aluminium [ 40% ], steel [ 15% ] and tissue [ 5% ] slabs are deduced in the same iterative methodology as for W-and T-slabs.Fig. 3c confirms that the IORT objective of a quasi-homogeneous dose within the tumor, followed by an immediate attenuation in the high-Z slabs inserted by the surgeon and finally a complete protection of healthy tissue is achieved for all beams.Fig. 4c illustrates that, for energies between 25 and 100 MeV , 97 to 99% of IORT-voxels meet the ǫ 2 criterion.This is followed by a quasi-plateau around 99%   of the VHEEs insensitivity to the boundary crossing effects and heterogeneity level is better translated by the ǫ 1 criterion.Figure 4c shows that the percentage of IORT-voxels meeting this criterion rises from 80.2 to 97.7% as energy increases from 1 MeV to 100 MeV .This percentage stabilizes around this value up to 1.0 GeV , then slightly decreases to 95.2% at 1.5 GeV .Localized losses of accuracy, amounting to 1.2 and 1.8% compared to the discussed performance, are observed at 300 MeV and 450 MeV , respectively, as highlighted by both ǫ 1 and ǫ 2 in Fig. 4c.The superior BFP-MC compliance on the IORT-benchmark, when compared to the T-one (Fig. 4b vs. 4c), can be attributed to the absence of the thin osseous slab.Averaging over all IORT-voxels, ǭ is 0.28% , 0.37% , 0.32% and 0.27% , respectively, at 100 MeV , 300 MeV , 500 MeV and 1 GeV. .The agreement between BFP and MC doses is maintained across all slabs, depths, buildups, densities, and beams.What was reported for the previous discussion regarding the CONV-VHEE transition remains true for HH-benchmark.Figure 4d indicates that for VHEE below 100 MeV , the ǫ 2 criterion is met by 93.7 to 96.8% of the HH-voxels.The fulfill- ment rate slightly rises to 96.8-97.7% for energies from 100 MeV to 1 GeV , then further increases and stabilizes around 98.1% up to 1.5 GeV .When the more stringent ǫ 1 criterion is applied, these percentages drop to 70.4% , 70.4-91.1%,and 92.9% , respectively.Restricting consideration to the ǫ 2 criterion, a slightly improved BFP-MC compliance can be readily observed in Fig. 4 during the transition from the T to HH-benchmark.This is due to the compensatory effects of backscattering from adjacent slabs, given the tripling of the total number of slabs during the T-HH transition.However, the ǫ 1 criterion indicates the opposite.This can be attributed to the sub- stantial difficulty in achieving compliance below 1% within tissues when transitioning from 4 slabs to 11 slabs.When averaged across all HH-voxels, ǭ is 0.55% , 0.61% , 0.68% and 0.70% , respectively, at 100 MeV , 300 MeV , 500 MeV and 1 GeV .All violations of the ǫ 2 criterion observed in the insets of Figs.3a-d are attributed to specific interface issues on the Geant-4 side.
The hypothesis of a BFP-MC dose equivalence fails to be rejected by the Kolmogrov-Smirnov goodness-of-fit test at a significance level of ∼ 0.02 below 1.5 GeV for all benchmarks, indicating a high level of BFP-MC compli- ance.The Dragon-5 CPU time, which depends solely on deterministic parameters ( N g , P l , S N , and D r ), remains consistent at 5 seconds across all beams.In contrast, the Geant-4 CPU time increases with the beam energy, with estimated times of 10.4 days, 24.9 days, and 36.36 days at 100 MeV , 500 MeV , and 1 GeV , respectively.These comparisons were conducted on an Intel Xenon E5-2683 v4 CPU.www.nature.com/scientificreports/

High-Z benchmarks
We propose to amplify the complexity of transport by incorporating slabs of high heterogeneity, high density and high atomic number.Simultaneously, we enhance unidirectional beam incidences by initiating irradiation from both sides of the benchmarks.Figure 5a 5d illustrates that the percentage of NM1-voxels satisfying the ǫ 2 criterion fluctuates around 98.8% up to 1.5 GeV .Beyond this, a continuous decline is observed, reaching 95.6% , 79.7% , 73.1% , 69.3% , and 68.7% at 2 GeV , 3 GeV , 4 GeV , 5 GeV , and 6 GeV , respectively.The ǫ 1 criterion demonstrates a significant loss of accuracy for certain well-distinguished energies, while for other energies, the criterion remains satisfied at 95.1-97.5%. Figure 5a reveals that significant BFP-MC discrepancies are largely localized within the carbon slab for specific beams, most notably at 300 MeV , 400 MeV , and 850 MeV , but this pattern is not universal.Our impending meta-analysis will underscore that carbon does not pose any significant challenge for the BFP solution within VHEE and UHEE domains.However, the origins of these discrepancies within specific beams remain unclear.When averaged across all NM1-voxels, the mean BFP-MC discrepancy ǭ stands at 0.34% , 0.44% , 0.51% , and 0.65% for 100 MeV , 300 MeV , 500 MeV , and 1 GeV , respectively.Figure 5b presents dose profiles for the NM2 benchmark: As observed with the NM1 case, NM2 displays similar trends with subtle differences in conformity.In Fig. 5d, the BFP-MC conformity for the ǫ 2 criterion oscillates between 98.3 and 99.4% within the 100 MeV-1.5 GeV energy range.Above this range, we observe a consistent decline.Below it, the BFP-MC conformity to the mentioned criterion climbs from 86.1 to 98.3% .With the ǫ 1 criterion, we can distinguish the clear pattern of BFP-MC convergence that escalates with increasing beam energy: BFP-MC conformity rises from 77.2% at 61 MeV to 95.0% at 1 GeV , beyond which the conformity starts to diminish.Transitioning from NM1 to NM2 brings increased complexity in transport, seen in the expanded number of slabs and the more intricate irradiated materials.This increase in complexity correlates with a decrease in conformity for the stringent ǫ 1 criterion.The diverging NM2-voxels are located in the fourth slab of francium, which possesses the highest atomic number among all NM2 slabs.The upcoming meta-analysis will show that francium categorically exhibits a BFP-MC deviation above 2% at the point of maximum energy deposition for all VHEE and UHEE beams.Taking an average over all the NM2 voxels, ǭ measures 0.62% , 0.52% , 0.60% , and 0.43% at 100 MeV , 300 MeV , 500 MeV , and 1 GeV , respectively.
Figure 5c relates to the benchmark with the highest level of complexity, NM3, composed of 15 slabs, from left to right: The complexity of transport is illustrated by compliance with both ǫ 1 and ǫ 2 criteria, as depicted in Fig. 5d.Between 100 MeV and 1 GeV , compliance with the ǫ 2 criterion fluctuates between 87.2 and 87.9% , but it drops significantly to 20.5% at 1.5 GeV .Conversely, it increases from 80.4 to 87.2% as energy transitions from 45 to 99 MeV .It was previously demon- strated that 1.5 GeV was the maximum limit to which the CEPXS-mode could be extended from its design limit of 100 MeV .This was the case for all previous benchmarks, except for NM3.For the latter, the compliance with the ǫ 1 criterion is the lowest observed, with a maximum of 68.0% at 150 MeV (Fig. 5d).Notably, voxels show- ing significant deviations are predominantly located in the seventh slab of phosphorus, the largest slab of the benchmark.Intriguingly, just as observed in NM1, phosphorus does not pose a transport challenge.Hence, the observed deviation potentially implicates the type of geometric encapsulation of the slab.This scenario underscores the NM3-benchmark's need for meticulous optimization of deterministic computation schemes, emphasizing spatial discretization, quadratures, and anisotropy order, to enhance BFP-MC compliance.Averaged over all NM3 voxels, ǭ are 1.35% , 1.38% , 1.35% , and 1.36% at 100 MeV , 300 MeV , 500 MeV , and 1 GeV , respectively.

VHEE and UHEE meta-analysis
Here, we probe the boundaries of the extended CEPXS-mode in ELECTR at VHEE and UHEE, scrutinizing the Class concept's validity from CONV-RT 104 .We irradiate the entire periodic table, from hydrogen to einsteinium, with unidirectional beams under VHEE and UHEE conditions.For each irradiation, an iterative procedure autonomously determines the slab dimensions in Geant-4 and Dragon-5, taking into account the beam energy, irradiated material, and corresponding stopping power.Figure 6a shows the behaviour of criteria ǫ 1 and ǫ 2 as a function of Z, while Fig. 6b highlights the same behaviour for ǭ .Figures 7-8 illustrate the spatial distribution of ǫ r for most elements.
Tracking the behavior of the maximum relative error across Figs.7-8, we observe the compliance of the BFP-MC deviation with respect to the ǫ 2 criterion.We note that the maximum deviation is often located at the benchmark's hottest point, i.e., the point of maximum energy deposit ( d max ).Except for a few exceptions, Figs.7-8 show that the ǫ 2 criterion is respected at VHEE and UHEE from Z = 3 (lithium) to Z = 58 (cerium).This con- sistency constitutes what we define as the first response Class, C 1 , which aligns with the criteria established for CONV-RT 104 .Illustratively, for C 1 -elements, and within an energy range of 30 MeV to 1 GeV , the combined maximum deviation for all voxels does not exceed 0.6% for 5 B, 0.8% for 14 Si, 1.6% for 21 Sc, 1.2% for 29 Cu, 1.2% for 41 Nb, 1.9% for 48 Cd and 1.98% for 57 La.
From Z = 59 (praseodymium) to Z = 92 (uranium), we observe the formation of a second response Class, C 2 .Three distinct characteristics define this Class: (1) the systematic appearance of a deviation above 2% at d max ; (2) the persistence of this behavior across all C 2 -elements and all beams studied; and (3) a gradual yet consistent increase and lateral expansion in this deviation with increasing Z.
A shift towards a third response Class ( C 3 ) is observed beginning at Z = 93 (neptunium) and continuing up to Z = 99 (einsteinium).In C 3 , the BFP-MC deviations initially quadruple, then increase tenfold as Z increases.in the violation of the ǫ 2 criterion (Fig. 8).We therefore conclude that, under VHEE and UHEE conditions, the extended CEPXS-mode is operational for C 1 , rejected for C 2 and entirely ineffective for C 3 .The deviation from the anticipated behavior of C 1 is manifested across four distinct categories of elements: three noble gases-helium ( 2 He), neon ( 10 Ne), and argon ( 18 Ar); three halogens-fluorine ( 9 F), bromine ( 35 Br), and iodine ( 53 I); two gases, namely hydrogen ( 1 H) and oxygen ( 8 O); and two metals-lithium ( 3 Li) and boron ( 5 B).The transitions between the response Classes-from C 1 to C 2 , and subsequently the fall into C 3 -are dis- tinctly depicted across the full VHEE and UHEE range, as showcased by the ǫ 2 criterion in Fig. 6a.Additionally, the ǫ 1 criterion (Fig. 6a) elucidates how the convergence of BFP-MC strengthens with increasing beam energy within the same Class.Remarkably, a significant decrease in both ǫ 1 and ǫ 2 values within the same Class (Fig. 6a) facilitates the identification of the previously mentioned exceptions to our classification methodology.In terms of ǭ , Fig. 6b shows that the C 1 -elements demonstrate a consistent ǭ below 1.0% , clearly distinguishing them.The transition to C 2 is marked by a noticeable shift to ǭ = 1.0%, which then gradually settles around ǭ = 1.98% .The move to C 3 is even more dramatic, as depicted in Fig. 6b.This transition initiates an abrupt increase from ǭ = 2.00% for 92 U to ǭ = 10.73% for 93 Np, yet it doesn't inhibit the recording of an extraordinary average devia- tion of 1719.88% for 98 Cf.Ultimately, the irregularities observed in Fig. 6a of an atom relative to its neighbors are intrinsic to the CEPXS mode.Our ongoing study indicates that the ENDF mode in ELECTR fully resolves these irregularities and failures beyond Z=50.

Study limits
Transverse transport: A deliberate decision was taken to validate the proposed extension through deterministic transport, rather than contrasting multigroup cross sections.Therefore, there's a skipped layer which involves extensive comparison of multigroup cross-sections-namely catastrophic transfer matrices (Eqs.3-22), total catastrophic cross-sections (Eq.23), soft stopping powers (Eqs.24-25), and energy deposition cross-sections (Eq.26)-with their credited, certified, or evaluated counterparts.It was then essential to pinpoint the pure influence of multigroup atomic data on electron transport.One approach to accomplish this is by restricting transport to one dimension, thereby eliminating any external interference or additional factors that might affect the direct observation of the multigroup cross sections' impact on BFP transport.BFP 1D resolution is known to be exact.In contrast, 2D and 3D transverse transport intervene with ray effects and SN catastrophic kernel discretization effects.Hence, it is imperative that 2D and 3D deterministic numerical challenges don't skrew the relative BFP-MC discrepancy, whether by amplifying it through additive factors or lessening it with compensatory ones.Transverse electron transport interferes our main objective: to clearly identify and separate the effects of multigroup cross-sections on electron transport.Dragon-5 transverse transport capacity can be found elsewhere (e.g., Fig. 11, p.12 of Ref. 169 ).Moreover, since the cross sections are spatially anisotropic, the transverse transport capabilities are not affected by the cross sections.
Photon transport: The BFP dose calculation is predetermined by the energy deposition cross section (Eq.26).This quantity is (1) complex; (2) deterministic; and (3) has no equivalent in MC paradigm.The successes and failures reported in this study are more a reflection of the quality of this specific cross section.Therefore, to validate "pure electron" transport, it is important to maintain a pure electron energy deposition cross section.The production of fluorescence and bremsstrahlung spectra is already implemented in ELECTR.A coupled photon-electron calculation in Dragon-5 is straightforward.However, from the Njoy perspective, it requires the coupling of two multigroup nuclear data modules; GAMINR for photons and ELECTR for electrons.Such coupling is essential to address the photons that emerge from electron transport-including fluorescence and bremsstrahlung-as well as the electrons produced by photon interactions, namely the photoelectric effect, Compton scattering, pair production, and triplet production.This coupling goes beyond the scope of this study and has not yet been performed.
Crossing boundary effect: A portion of the BFP-MC discrepancy is due to Geant-4's inability to accurately predict the dose at the interface point.In contrast, Dragon-5's predictions at the interface are highly realistic.This topic has been thoroughly addressed in our previous paper (see Figs. 3 and 4 and p. 12-13 of Ref. 104 ).Two methods can be applied to correct Geant-4's predictions: (1) by forcing the electron step size not to exceed 25% of the voxel's thickness; or (2) by optimizing the control constants of the G4UrbanMscModel condensed history (CH) algorithm.Such modifications to the MC computational scheme are recognized as extremely timeconsuming.Therefore, they are not practical for VHEE and UHEE beams, which are already very demanding in terms of CPU time.

Conclusion
Very-High Energy Electron (VHEE) and Ultra-High Energy Electron (UHEE) treatments are currently in their advanced preclinical testing stage.The Monte Carlo (MC) solution currently serves as the only operational transport support for experimental trials and forthcoming clinical deployments.However, its computational intensity and time consumption make it less suitable for real-time, or on-the-fly, clinical applications.In this paper, we propose a potential alternative by extending the Boltzmann-Fokker-Planck (BFP) chain, Njoy-Dragon, for VHEE and UHEE beams.The energy range examined covers 1 MeV to 6 GeV .We validated our multigroup solution in comparison with Geant-4, which has been previously qualified against experiments for VHEE and UHEE conditions.The validation process was conducted along three fronts: the first involves typical radio-oncology benchmarks for increasingly complex scenarios, the second focuses on benchmarks characterized by high levels of heterogeneity and intricate atomic structure, and the third involves irradiation of the entire periodic table.
By comparing the dose in each voxel, we found that 99% of water voxels exhibited a BFP-MC deviation below 2% under 1.5 GeV .Similarly, applying the same criterion, 97% and 99% of thorax and breast intra-operative voxels, respectively, met this criterion above 50 MeV .In the thorax, we observed a loss of 1 voxel of accuracy, attributable to the bone-lung interface.However, this issue was not further investigated on the Geant-4 side.For a heterogeneity level consisting of 11 slabs-including muscle, adipose, lung, and bone, 97-98% of the voxels displayed a BFP-MC deviation below 2% between 60 MeV and 1.5 GeV .These findings suggest that the accuracy gain following the mentioned energy thresholds is likely due to the insensitivity of VHEE/UHEE to the level of heterogeneity.While BFP-MC conformity at 1% was analyzed, it is crucial to note that this criterion is lower than the uncertainties on cross-sections.Applying this criterion, we noted a decrease in voxel complicity: a reduction of 1% for water voxels, 1-8% for thorax, 2-12% for Mobetron, and 7% for the highly heterogeneous benchmark.The average BFP-MC discrepancy remained around 0.3% and 0.4% for water, thorax, and Mobetron at 300 MeV and 1 GeV .However, it slightly increased to 0.6% and 0.7% for the highly heterogeneous benchmark.
For the second study, we involve, permute, and assemble high and medium Z materials while gradually increasing the level of heterogeneity and multiplying the incident beams.Below 1.5 GeV , we achieved a minimum of 98.2% BFP-MC compliance to the 2% criterion in the Fe-As-C-Zr assembly.Meanwhile, for the Si-Mo-Cr-Fr- Mg-Cu assembly, compliance to the same criterion ranged between 97.8% and 99.4% within the [75 MeV, 1.5 GeV] interval.However, in the more complex Au-S-Zn-Sn-Na-Se-K-Sm-V-Pd-B-Y-In-Yb-Ti assembly, BFP-MC compliance failed to exceed 87.9% of all involved voxels and beams.This suggests that increased heterogeneity leads to a loss of accuracy at lower energy, particularly for the highly stringent 1% criterion.The average BFP-MC discrepancy for assemblies 1 to 3 was found to be 0.44% ( 0.65% ), 0.52% ( 0.43% ), and 1.38% ( 1.36% ) at 300 MeV ( 1 GeV ), respectively.VHEE and UHEE irradiation of the entire periodic table has identified three classes of response.For elements with atomic numbers lower than praseodymium ( Z = 59 ), a BFP-MC discrepancy below 2% was consistently achieved for all voxels within the 1 MeV to 1.5 GeV range.For elements from praseodymium to uranium, a systematic BFP-MC deviation slightly above 2% was observed at the point of maximum energy deposit.This deviation not only persisted but also increased and widened laterally as Z increased.A third category was identified between neptunium and einsteinium, where the discrepancies were at least four times as large.However, we noted some exceptions to our classification, such as hydrogen, the noble gases helium and argon, and the halogens bromine and iodine.
The present study avoided any interference or counterbalancing effect between errors related to the multigroup formalism, ray effect, or the S N discretization of the Boltzmann kernel.For this reason, deterministic transport was limited to one dimension.The optimization of the deterministic computational scheme must be adapted to clinical routine.Even better results can be reported in the future with a refinement of the computational scheme (spatial discretization, Legendre order, S N quadrature and number of groups) and with the Evaluated Nuclear Data File (ENDF) mode in Njoy.
) the quantum wave-particle duality at low energy, we restrict e δ ∈ [e c , e ′ e /2] .This immediately implies e e ∈ [e ′ e /2, e ′ e − e c ] .The cutoff energy ( e c ) is thus identified as the lowest energy a δ-ray could have as a result of a catastrophic i-collision.In multigroup theory, the ELECTR [CEPXS-mode] follows Lorence's 132 and Morel's133 proposal to consider a catastrophic collision, an i-event in which the scattered VHEE does not appear in two adjacent groups, i.e., e c = e ′ e g l [M] moments can be redefined as a linear combination of Spencer's integrals for Mott scattering 160 : https://doi.org/10.1038/s41598-023-51143-8

Figure 2 .
Figure 2. Simplified partial relaxation cascades from K, L 1 , L 2 , L 3 , M and N shells and subshells in ELECTR [CEPXS-mode].η a and η f refer, respectively, to Auger and fluorescence relaxation efficiency.Relaxation stops when the vacancy is transferred to the outermost subshell.The line radiation for a given relaxation is obtained by considering a conditioned multiplication of the probabilities of all the involved branches.The complete cascade is depicted by replicating the de-excitations at each call of the subshell.For illustration purposes, the duplication in this figure is demonstrated only for a few branches.