Feasibility of a multigroup Boltzmann–Fokker–Planck solution for electron beam dose calculations

Legacy nuclear-reactor Boltzmann solvers start clinical deployment as an alternative to Monte Carlo (MC) codes and Fermi–Eyges semiemprical models in radiation oncology treatment planning. Today’s certified clinical solvers are limited to photon beams. In this paper, ELECTR, a state-of-the-art multigroup electron cross sections generation module in NJOY is presented and validated against Lockwood’s calorimetric measurements, EGS-nrc and GEANT-4 for 1–20 MeV unidirectional electron beams. The nuclear-reactor DRAGON-5 solver is upgraded to access the library and solve the Boltzmann–Fokker–Planck (BFP) equation. A variety of heterogeneous radiotherapy and radiosurgery phantom configurations were used for validation purpose. Case studies include a thorax benchmark, that of a typical breast Intra-Operative Radiotherapy and a high-heterogeneity patient-like benchmark. For all beams, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$100\%$$\end{document}100% of the water voxels satisfied the American Association of Physicists in Medicine accuracy criterion for a BFP-MC dose error 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%. 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.0\%$$\end{document}97.0% of adipose, muscle, bone, lung, tumor and breast voxels satisfied the \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% criterion. The average BFP-MC relative error was about \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.56\%$$\end{document}0.56% for all voxels, beams and materials combined. By irradiating homogeneous slabs from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z=1$$\end{document}Z=1 (hydrogen) to \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 (einsteinium), we reported performance and defects of the CEPXS mode [US. Sandia National Lab., SAND-89-1685] in ELECTR for the entire periodic table. For all Lockwood’s benchmarks, NJOY-DRAGON dose predictions are within the experimental data precision for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$98\%$$\end{document}98% of voxels.

www.nature.com/scientificreports/ recovered by the upgraded Dragon-5 BFP discrete ordinate solver 67 , where a high-order diamond differencing (HODD) scheme is used for spatial discretization, with a P ≥8 Legendre expansion for scattering anisotropy and S 64 angular quadrature. Finally, we proceed with validations against Lockwood's experimental calorimetric measurements 68 , Egs-nrc, Geant-4 on water slabs, thorax 69 , intraoperative Mobetron electron radiotherapy 70 and high-heterogeneity benchmarks 41 as well as typical slabs covering the entire periodic table. This paper is limited to the case of CEPXS feed functions in ELECTR. For validation reasons, only electrons are transported. Bremsstrahlung and fluorescence photons are produced and eliminated at the place of birth.

Methods
Boltzmann-Fokker-Planck equation. Let n be the electron discrete ordinate direction ( ˆ n ), g its energy group, r its position, the multigroup form of the S N BFP equation is then given by: ψ e designates the electron flux, t the total cross section, L B the Boltzmann operator, L FP the Fokker-Planck operator and Q e,ext. the external electron source. All these quantities are discretized on the phase space D . According to Fano's [ 1 keV , 100GeV ] classification 71 , the considered interactions are the bremsstrahlung (b), inelastic collision/ionization (c), elastic collision (e) Auger, Coster-Kronig and super Coster-Kroning cascades (a).
Although fluorescence cascades are implemented in ELECTR , here we limit our study to pure electron transport. The electron flux satisfies the periodicity condition. Its development in spherical harmonics ( R m l ) is written as: where L refers to the Legendre order. Flux moments ψ m l,g ( r) can be obtained by integrating Eq. 2 over all incidence directions and normalizing with the complex conjugate R m * l . L B describes the rate of catastrophic collisions. L B has a rotation invariant kernel, i.e., its eigenfunctions are the R m l spherical harmonics and the associated eigenvalues are the Legendre scattering cross section coefficients, x s,l . x designates the electroatomic interaction, i.e., x ∈ {a, b, c, e} . Using Eq. 2, decomposing the scattering cross sections into Legendre polynomials, then using Legendre Addition Theorem and finally exploiting the orthogonality of the R m l harmonics, the Boltzmann operator is written as: g,n operator is obtained by a Taylor polynomial expansion of the Boltzmann L B g,n operator around small energyloss and small-angle scattering 72 . Considering only the first-order Taylor expansion, L FP g,n can be written as the sum of the continuous slowing-down ( L FP 1 g ) and the continuous scattering operator ( L FP 2 g,n ).
β g and α g designate, respectively, the restricted stopping power and the restricted momentum transfer. Restriction refers to soft collisions. Here, the considered source is monoenergetic and unidirectional. If Q e ( r) is the source intensity and D Q its spatial domain: The system Eqs. 1 and 5 is considered complete after applying the boundary conditions. [73][74][75][76] is used to define electroatomic multigroup Legendre coefficient transfer matrices ( x l ) and multigroup total cross sections ( t g ). A Gauss-Lobatto quadrature is used to evaluate all integrals in ELECTR. The nodes and weights are provided up to the 10th order. For a given interaction, x, microscopic catastrophic cross sections are given by:

Microscopic electroatomic catastrophic cross sections. Neutronic formalism
where g and g ′ refer to the incident and transfer energy groups, respectively. φ ℓ is a Legendre component of the electron flux's guess. For both collisional and radiative interactions, a catastrophic event is assumed to be one in which down-scattering happens into non-adjacent energy groups. F (E) is a unifying concept of the so-called feed function. This concept was introduced by MacFarlane et al. 77,78 in the NJOY GROUPR and GAMINR modules for neutron-and photon-production cross sections and is adapted here for electron-production cross sections. It alone changes for different data types. For matrices, F x,0 lg ′ (E g ) is the ℓ th Legendre moment of the normalized probability of scattering into secondary energy group g ′ from initial energy E g . For vectors, F x,0 (E g ) is the total (1) � n · � ∇ψ e g,n (� r) + � t,g (� r)ψ e g,n (� r) = L B g,n ψ e g,n (� r) + L FP g,n ψ e g,n (� r) + Q e,ext. g,n (� r).
(2) ψ e g,n (� r) L FP g,n = L FP1 g ψ e g,n (� r) + L FP2 g,n ψ e g,n (� r), where L FP1 www.nature.com/scientificreports/ cross section interaction from the same initial group E g . The general form of the mth moment of the ℓ th Legendre order group-to-group feed function (MeV m · barns) is given by: where ˆ and ˆ ′ stand for the incident and scattered particle directions, respectively. P l is the Legendre polynomial. I x g ′ is the electron transfer domain for interaction x. The differential scattering cross section (DSC) is given by: P x is the microscopic differential energy distribution kernel, while D x refers to the microscopic differential angular distribution kernel (both for the x interaction).
Ionization feed function. Two electrons emerge from each kth atomic subshell ionization by an incident electron of energy E g : the scattered electron and the delta ray (also called recoil electron). By convention, the particle with the lower energy is identified as the delta ray, and the one with the higher energy is the principal scattered electron. The total ionization feed function for group g ′ is then the sum over all subshells of the mth moment of the ℓ th Legendre order for primary electron transferring to group g ′ ( F k,ee,m ℓg ′ (E g ) ) and the mth moment of the ℓ th Legendre order for delta ray production and transferring to group g ′ ( F k,eδ,m ℓg ′ (E g ) ), all starting with incident energy E g impacting on the kth subshell: By substituting Eq. 7 in 9, the questions then become (i) how ELECTR proceeds for P k,ee/eδ and D k,ee/eδ distributions (Eq. 8); and (ii) what are the transfer domains I k,ee g ′ and I k,eδ g ′ (Eq. 7) for both scattered and delta electrons. It depends on the operation mode, i.e., if it is an ENDF or a CEPXS feed function mode. Here, the operation "mode" should not be confused with the data "format". In both cases, the energy and the scattering angle of the primary and secondary electrons are deduced from the conservation laws of the collision kinematics. In either case, it is assumed that no kinetic energy is transferred to the residual atom. For the ENDF mode, three-body collision kinematics are used, namely: the incident particle, the secondary emitted one and the subshell involved. P k,eδ distributions are all provided in the ENDF evaluation as MF=26, MT=534-572 for all materials. The MF and MT values are part of the ENDF-6 format for evaluated nuclear data. The MF value, running from 1 to 99, encodes the class of information (e.g., MF=23: interaction cross section, MF=26: angular distribution). The MT value, running from 1 to 999, encodes the interaction type (e.g., MT=102: neutron capture, MT=527: electron bremsstrahlung). D k,ee and D k,eδ are both Dirac delta distributions around the cosine of the scattering angle of the emitted particle. The subshell's binding energy and ionization total cross section ( σ k (E) in Eq. 8) are, respectively, provided in MF=28, MT=534-572 and MF=23, MT=534-572. The transfer domains for the primary and secondary electron are given by: where E cut is a threshold integration limit chosen to avoid feed functions' divergence for catastrophic collisions. E/2 is the lowest kinetic energy that the principal scattered electron can have. For the CEPXS mode, the ionization is done on a free electron and therefore there is no threshold condition nor any atomic shell or subshell involved. For both the primary and secondary electrons, the DSC (Eq. 8) is that of Møller, for which the angular distributions remain Dirac distributions around the cosine of the scattering angles obtained, in this case, by a two-body kinematics 79 . Unlike in condensed-history MC codes, energy-loss straggling models are not necessary in ELECTR since differential P k,eδ distributions already account for this effect.
Bremsstrahlung feed functions. One electron and one photon emanate from this interaction. Whether for ENDF or CEPXS mode, no angular deflection is allowed for the bremsstrahlung electron. D b in Eq. 8 is therefore a Dirac delta distribution around ˆ ·ˆ ′ = 1 . The photon, for both ENDF and CEPXS modes, is emitted according to a Sommerfeld angular distribution 80 . The energies of the two particles are obtained by the conservation balance. Eqs. 7 and 8 can be simplified, in this case, as follows: E min is the high frequency limit of bremsstrahlung divergence. It is fixed at 1 keV for both modes. E cut has the same definition as in Eq. 10. For the ENDF mode, P b and σ b in Eq. 11 are interpolated, respectively, from MF=26, MT=527 and MF=23, MT=527, taking into account bremsstrahlung production in nucleus and Coulomb electric fields. For the CEPXS mode, the distribution is implemented analytically according to Koch and Motz DSC assembly based on relativistic Born approximation with Coulomb correction 81 . Empirical parameters, e.g., Elwert corrections and screening factors, are retrieved from CEPXS Tape9 database 66 . The Koch and Motz DSC also diverges when the energy of the photon approaches that of the incident electron 82 Scientific Reports | (2023) 13:1310 | https://doi.org/10.1038/s41598-023-27376-y www.nature.com/scientificreports/ introduced in Eq. 11 is also used in this mode. Bremsstrahlung production from interaction with atomic electrons' field is accounted by modifying the Koch and Motz Z 2 factor by Z(Z + 1). Elastic feed function The nucleus-scattered electron is deviated without change in energy ( P e = δ gg ′ ). Equations 7 and 8 are then simplified as follows: For the ENDF mode, σ e (E g ) is interpolated from MF=23, MT=525. For large angular scattering ( µ ∈ [−1.0, 0.999999] ), the D e distribution is interpolated from MF=26, MT=525 83 . For forward scattering ( µ ∈ [0.999999, 1.0] ), Seltzer's analytical screened Coulomb DSC 84,85 is implemented. From Cullen's EEDL classification 83 , elastic scattering is considered either large or forward-peaked. For all atoms, ENDF angular distributions are provided below and above 256 keV . In the CEPXS mode, the Mott DSC with Molière screening 86 is used for relativistic energies ( E g > 256 keV ), while a Riley DSC 87 is used for lower energies. For both Mott and Riley distributions, and instead of evaluating the Legendre moments in Eq. 12 with quadratures, the CEPXS mode uses Berger's semi-analytical approach 88 based on Goudsmit-Saunderson distribution and Spencer functions 89 . All empirical parameters are retrieved from CEPXS database 66 . The final step in both modes is to modify the lth Legendre order of the elastic feed function by an extended transport correction, similar to the one proposed by Bell for neutrons 90 in nuclear reactor physics. The transport-corrected within-group elastic feed function is given by: where L is the maximum Legendre order. As shown by Morel for electrons 91 , such an approximation makes the highly forward-peaked elastic scattering kernel compatible with low-order Legendre expansion and, potentially, a replacement for L FP 2 in Eq. 4.
Auger relaxation cascade. In both modes, it is assumed that Auger electrons are emitted isotropically. Let j be the transition line, η e kj its relaxation efficiency (i.e., that the j th transition produces an Auger electron following the kth shell or subshell impact ionization), the m th moment of the l th Legendre order of Auger production feed function is given by: Higher-order Auger production feed functions are identically zero. g j is the electron group that contains the j th emitted Auger energy, while e e kj refers to the actual Auger, Coster-Kronig or super Coster-Kronig electron energy. N tr is the number of possible atomic transitions. For the ENDF mode, transition probabilities and the emitted particle spectra are both interpolated from MF=28, MT=534-572, which describe subshell ionization for the K1 to O5 shells. For the CEPXS mode, only the K, L1, L2, L3 and M shells are involved in relaxation cascades, which results in N tr = 28 different transitions after an impact event. Also, for this mode, e e kj is expressed in terms of the midpoint group energy rather than the line radiation one. Binding energies and probability parameters are stored in Tape9 relaxation quantities 66 . For both modes, the relaxation algorithm does not contribute to the total reaction rate.
Total macroscopic cross sections. For all interactions, the total catastrophic cross section can be obtained by summing the group-to-group feed function (Eq. 7) over all incident angles: The upper limit of I x g ′ (Eqs. 15 and 8) is set in such a way that the threshold energy is arbitrarily selected at the interface E g−1 between soft and catastrophic collisions. Its lower limit is the lowest energy that the scattered electron could have (e.g., E/2 for Møller scattering).
Microscopic restricted stopping powers. Total stopping powers (MeV · barns) are converted in ENDF-6 format from Berger's evaluation 92 . ELECTR interpolates collisional ( S c g ) and radiative ( S b g ) total stopping powers from the converted MF=23, MT=507 and MF=23, MT=508 sections, respectively. Since Bethe S c g shows imperfections below 10 keV , due to negligence of deviations from plasmon, inner subshells' electrons and conduction electrons 93 , ELECTR corrects S c g below 10 keV using Lorence's power-law extrapolation 66 . Catastrophic stopping powers, M c g and M b g , are taken as the first moment of the total catastrophic cross section. Before being averaged, restricted stopping powers are evaluated at all group boundaries except the last one. www.nature.com/scientificreports/ where, for M g , we followed the state-of-the-art recommendations by integrating Møller (in Eq. 16) and Koch and Motz (in Eq. 17) DSCs. Integrals' limits are identical to the previous ones, i.e., a non-divergence lower limit and a soft-catastrophic interface upper limit. The total restricted stopping power (Eq. 4) is the sum of Eqs. 16 and 17. In this study, lack of stopping powers below 1 keV impose a transport cutoff at 1 keV . This cutoff is in line with recommendations made by Salvat to Cullen for EPICS-2017 data 83,94,95 . It is also in compliance with the default cutoff limits of Penelope 96 , Egsnrc 97 and Fluka 98 .

Microscopic energy and charge deposition cross sections. Energy deposition cross section (MeV ·
barns) account for the energy deposited in each electron-matter interaction. It can take either a positive (local energy deposition) or a negative (local energy removal) value. The energy transported with the scattered electron or the emitted Auger, Coster-Kronig, or super Coster-Kronig electron are typical examples of negative energy deposition cross sections. This quantity is central for a deterministic dose assessment in when solving the BFP equation. Three phenomena are contributing: inelastic collision, bremsstrahlung and Auger cascades. In the first two cases, one should sum contributions from soft and catastrophic collisions. Inelastic collisional and bremsstrahlung energy deposition cross sections are given by Relaxation energy deposition is simply the negative first moment of the Auger production feed function (Eq. 14).
Charge deposition cross sections account for single electron reckoning, i.e., its deposition and removal, from the interaction site. They are therefore obtained in the same manner as to Eqs. 18-20 depositions, except with zero moments of the feed functions. Charge deposition feed functions are, by definition, associated with non-zero absorption interactions, and are given by: Both E x g ′ and C x g ′ concepts were introduced by Lorence et al. 66 , used by Acuros for within patient dose distribution and damage quantification 53 and are generalized here for the ELECTR module. Dose deposition. The dose deposited is calculated in the Dragon-5 solver. Two quantities are needed: (1) the multigroup flux distribution ( g ) obtained in Dragon-5 after solving Eq. 1 and integrating over all directions; and (2) the energy deposition cross sections provided by ELECTR (Eqs. [18][19][20]. Let T be the irradiation time and ρ the voxel density, the total dose distribution (in Gy) is given by: where E g refer to the midpoint energy for group g and E g to the flux-weighted average energy for group g. Eq. 24 shows three components: (1) catastrophic energy deposition; (2) continuous slowing down ( β g � g ) and 3) energy deposition below the cutoff ( β N � N ). If an electron reaches an energy below 1 keV , it is stopped and its energy is deposited locally. As explained by Morel et al. 99 , by manipulating Eq. 24 and using Eqs. 18-20, the dose deposited is simplified as follows: Algorithm coding and computational scheme. A pointwise-ENDF (PENDF) file, which contains electron cross sections on a unionized and linearized energy grid, is first generated using the NJOY RECONR www.nature.com/scientificreports/ For each interaction process, the NJOY GROUPR module panel logic is adapted to average microscopic catastrophic cross sections (Eq. 6). Since each integrand in Eq. 6 has its own characteristics, the first calls and operations aim to unify the integration grid and quadrature schemes and detect discontinuities. The etsig subroutine is called to interpolate PENDF σ x at E g incidence energy and returns the next grid energy. The etflx and eetff subroutines perform similar operations for the electron flux and feed functions. The subroutine eetsed is called in both eetff and epanel cores. In the first call, scratch storage is allocated and all the ENDF subsections are read in. In the second call, angular and energy distribution kernels ( P x and D x in Eq. 8) are interpolated using the unit-base interpolation technique. A union grid is then generated for the three integrands σ x , φ x and F x . The integration algorithm is equivalent to the Minx legacy adaptative algorithm 100 .
Group-to-group catastrophic transfer matrices are computed for all secondary energy groups and all Legendre orders simultaneously. After all reactions are processed for this atom, a special pass is called to compute and store total catastrophic cross sections, total energy and charge deposition cross sections, as well as total restricted stopping powers. The NJOY upgraded MATXSR post-processing module formats the generated GENDF library in a single stand-alone multigroup dataset that can be accessed by a variety of legacy discrete ordinates and lattice codes. The MATXS library is used in conjunction with the Dragon-5 BFP solver 67 . Our deterministic computational scheme goes as follows. Depending on the mixture density and composition, microscopic cross sections are extracted, prepared and converted to macroscopic ones, respectively, by the LIB: and MAC: modules. Medium self-polarization reduces collisional stopping power. The latter appears in both restricted stopping power (Eq. 16) and energy deposition cross section (Eq. 18). A Sternheimer density effect correction 101,102 is implemented and applied in Dragon-5 on both quantities (STERN 1). Region mixtures with vacuum boundary conditions are defined by the GEO: module. Discrete ordinate integration lines, region identification pointers and all tracking information are generated by the SNT: module. A High-Order Diamond Differencing (HODD) is used for the spatial discretization with a linear Legendre spatial order 103 . A convergence criterion of 1 × 10 −5 was imposed on the flux inner iterations. The PSOUR: module is called to compute Eq. 1 right-hand-side source term. Since we limited our study to pure electron transport, PSOUR: will not compute corresponding bremsstrahlung and fluorescence photon sources. The ASM: assembly module recovers tracking lengths and material numbers from previous sequential tracking and computes BFP group-dependent system matrices. The linear BFP multigroup equation is finally solved by the FLU: module. Multigroup electron flux and macroscopic energy deposition cross section are used to compute dose spatial distributions by the HEAT: module.

Results
Zero to medium heterogeneity benchmarks. We propose a validation of the NJOY-Dragon chain with a set of an increasing complexity benchmarks. Furthermore, a gradual increase in the beam energy is considered to challenge Legendre anisotropy, S N order, deterministic transport correction and number of energy  is solved in 1D. This strategy makes it possible to study the net performance of the multigroup library on the depth dose profile without there being an accumulation of numerical errors related to the resolution of Eq. 1 in 2D and 3D. We thus avoid ray effects, an exorbitant angular discretization of the transport kernel and negative consequences from the high dimensionality of the Legendre scattering matrices. The use of unidirectional beams increases the computational complexity as well as the CPU time for two reasons. The first is that more focused interactions per unit length are imposed. The second is that the counterbalancing effects of dose errors at the interfaces and at the buildup for isotropic incidences are avoided. The W-benchmark is a rudimentary case study, which remains of clinical interest for both daily calibration verification and treatment planning. The T-benchmark introduces a first level of complexity with high-and low-density thin heterogeneities. The latter brings out the interface effects. The IORT-benchmark introduces a second level of within patient complexity combining internal (tissues) and external (instrumental) heterogeneities. High Z heterogeneities highlight shielding effect before a typical OAR. Incident beam energies vary from 1 to 20 MeV , covering the entire energy spectrum of clinical and medical research accelerators 104 . This also covers the entire electron radiotherapy and radiosurgery clinical beams' spectra. Although electron beam treatments represent only 10-15% of daily workload in clinical practices 14 and are constantly declining with the development of IMRT 105 , developing pure electron transport capabilities remains the preliminary step before fully-coupled [ γ , e − , e + ] transport for photon beam treatments.
Here, our goals are to: (1) demonstrate the accuracy of the ELECTR module, and (2) highlight, if any, the limitations of the CEPXS feed functions. As an immediate consequence, reducing the CPU time is not of interest in this study. For this reason, further optimizations of deterministic parameters ( N g , P l , S N , D r ) are not discussed here, where D r refers to the discretization of the spatial domain. Conservative parameters are selected to eliminate deterministic bias when evaluating the accuracy of the ELECTR module. All computational schemes therefore relied on N g = 300 groups, S 48 order and D r ≥ 500 elements. Legendre anisotropy is fixed to P 9 for 1 to 6 MeV , P 10 for 7 to 9 MeV , P 11 for 10 to 14 MeV and P 12 for 15 to 20 MeV electron beams. Convergence calculations were performed to support such choices. Only electrons are transported. Incident beams are assumed clean without any contamination.
For the Geant-4 reference scheme,15 × 10 6 primary events are considered for a 0.2% statistical uncertainty (or better) in every material. Unless otherwise specified, we consider the physics constructor G4EmLivermore based on EPICS data (EADL 94 , EEDL 83 and EPDL 95 ). To be consistent with Dragon-5, transport cutoff and secondary production threshold are both fixed to 1 keV . The identity of each secondary particle is checked at each step of the track. The fluorescence and bremsstrahlung photons are thus eliminated at the point of birth. Figure 2a shows Dragon-5 versus Geant-4 depth-dose curves for the W-benchmark (Fig. 1a). As the beam's energy increases, build-up spreads and the Legendre requirements are therefore higher. We observe that the buildup inflection dictates the P l order requirement. This can be explained by Mott and Møller's increasing anisotropy with the beam energy. The insert shows the relative BFP deviation with respect to the MC scheme. This was obtained, in post-processing, following a piecewise cubic Hermite grid unification of Geant-4 and Dragon-5 dose detectors. The average BFP-MC error (in absolute value) is of 0.43% , 0.58% , 0.52% and 0.53% , respectively, at 1, 9, 15 and 20 MeV. This average is on the entire spatial domain. Table 1 shows voxels' percentage satisfying the 2% criterion. For all beams examined, 100% of the water voxels satisfy this criterion. This agreement is reduced to 85% for a 1 % BFP-MC deviation criterion. Figure 2b shows Dragon-5 vs. Geant-4 depth-dose curves for the T-benchmark (Fig. 1b). The Legendre anisotropy increase with energy is true for all benchmarks. At 9 MeV , the mean BFP-MC error is 0.76% , 0.53% , 0.68% and 0.50% , respectively, in the first water slab, bone, lung and the last slab of water. For the same materials, these deviations are of the order of 0.62% (0.62% ), 0.93% (0.80% ), 0.48% (0.40% ) and 0.57% (0.62% ) for the 15 MeV ( 20 MeV ) beam. From Table 1, 98.2% , 99.2% , 98.4% and 99.4% of the thorax voxels satisfy the 2% criterion, respectively, at 1, 9, 15 and 20 MeV . At 1 MeV , there is a loss of precision in lung voxels compared to higher energy beams. The slight loss of BFP-MC agreement, in some voxels, is caused by the bone interface. If one disregard the bone interface discontinuity, 99.8% of the thorax voxels satisfy the 2% criterion (all beams combined). Figure 3 points out that the way Geant-4 handles this discontinuity is ambiguous. We tested other physics list constructors, namely G4EmPenelope, G4EmLow-Physics, G4EmStandard and G4EmStandard_opt4. The same failures, as those of G4EmLivermore (Fig. 3), persist. This Geant-4 failure can be corrected either (1) by forcing the electron step size not to exceed 25% of the voxel's thickness; or (2) by optimizing control constants of the G4UrbanMscModel condensedhistory (CH) algorithm 106 . The Geant-4 team advises option 1. Although more physical, the option 2 is very time-consuming 107 . Finally, the post-processing piecewise cubic Hermite interpolation of the Geant-4 dose for a unique voxel grid MC-BFP comparison causes this single anomaly to propagate at intermediate voxels, thus further reducing the percentage at Table 1.
The same depth-dose curves for the IORT-benchmark (Fig. 1c) are depicted in Fig. 2c. Doses at the aluminum interfaces as well as within the tumor buildup are in close agreement between Dragon-5 and Geant-4. For the 9 MeV beam, the average BFP-MC deviation is 0.60% , 0.98% , 0.00% and 0.00% , respectively, in tumor tissue, aluminum, steel and breast tissue. These deviations are of 0.47% (0.46% ), 0.99% (0.92% ), 0.00% (0.00% ) and 0.00% (0.00% ) for the 15 MeV ( 20 MeV ) beam. Both Dragon-5 and Geant-4 electron doses are zero in the steel slab and breast tissue, which means that a complete attenuation of the beam was achieved in the aluminum slab for  Figure 4 confirms that the BFP-MC discrepancies are observed within 13 Al buildup and not at the tumor-13 Al interface.  www.nature.com/scientificreports/ Using other physics list constructors will not fix this error. We will show that there is an acceptable match (less than 2% ) between the dose predicted by Geant-4 and Dragon-5 in homogeneous 13 Al slabs for all depths and beams. The 13 Al buildup error can be fixed, as for the T-bone (Fig. 3), either by forcing tracking or refining the condensed history method. As soon as the particle enters a new volume or starts a new track, its step s is re-initialized according to s = f r max(R, 1 ) , where R is the electron range and f r ( ∈ [0, 1] ) is a control constant. Here, all cross sections are assumed to be constant along s. It is advisable either (1) to restrict the value of s so that the stopping power does not vary by more than 20% during the step 108 or (2) to reduce s to less than 25% of the detector's width. Therefore, the public SetStepFunction of the base class G4VMultipleSacttering should be called for further electron step control. Two internal parameters need to be finely optimized: (1) dRoverRange which imposes a maximum step size given by the ratio step/range and (2) the finalRange. The MC electron transport then becomes extremely time-consuming; resolving this is beyond the scope of the current work.
Moreover, there is no clinical interest in being precise in the prediction of the dose deposition in the 13 Al slab. However, the partial loss of precision in tumor tissue ( 0.5% of voxels, Table 1) is caused by the 13 Al backscattering effect. The maximum error remains below 2.48% (all voxels, beams and materials combined). One can question the validity of the 1 keV transport cut-off imposed by the CEPXS data. The range of a 1 keV electron in tissue is 5 µm . Taking the worst-case scenario studied here, a 1 MeV beam, typical for epithelial nonmelanoma skin cancer, will have a range of 0.5 cm . With the density of 500 voxels used in Dragon-5, this systematically results in a voxel having a 10 µm thickness; that is to say a voxel twice larger than the range of the electron.
High-heterogeneity benchmark. Dragon-5 depth-dose curves compared to Geant-4 ones on the high-heterogeneity (HH) patient-like benchmark (Fig. 1d) are depicted in Fig. 2d. A good BFP-MC agreement is observed for any beam, material, interface and/or buildup. This benchmark adds a new level of complexity greater than the previous thorax and IORT case studies. Electron backscattering effects are more frequent, more intense and have a direct consequence on the dose. Unlike the T-benchmark (Fig. 3), the insert in Fig. 2d shows that interfaces' discrepancies are much lower and we do not observe failures at the bone interfaces. This can be explained by errors' counterbalancing from the neighboring slabs' backscattering, given that we go from 4 (T) to 11 (HH) slabs. The average BFP-MC relative error is of 0.55% , 0.56% , 0.57% and 0.57% , respectively, at 1, 9, 15 and 20 MeV . Table 2 shows the percentage of voxels satisfying the 2% criterion for each of the 11 slabs. For the entire HH-benchmark, a minimum of 99.18% of voxels satisfy the 2% criterion, all beams combined. The percentage of the HH voxels that satisfy a limit of 1% BFP-MC discrepancy is about 77.82% , 78.55% , 84.36% , 82.55% , respectively, for 1, 9, 15 and 20 MeV electron beams.
The limits of the CEPXS mode in ELECTR: a meta-analysis. Previous benchmarks are limited to the radiation oncology context. We now answer the question of what would be the limit of the CEPXS mode in ELECTR for beams from 1 to 20 MeV . We then propose an electron irradiation of successive homogeneous slabs covering the entire periodic table. The slabs' dimensions are always fixed to the range of the incident beam. Transport remains restricted to electrons without considering the contamination effects. Figures 5, 6 and 7 show BFP-MC relative error in energy deposition profiles as a function of material depth from Z = 3 (lithium) up to Z = 98 (californium). Energy deposition is followed until complete beam attenuation is achieved. As a criterion for acceptable accuracy, we take a MC-BFP discrepancy below the threshold of 2% , i.e., in compliance with radiation oncology standards 14,50 . Apart from some exceptions, which we elaborate on in what follows, the criterion of 2% is met for almost all voxels and beam energies, from Z = 3 (lithium) to Z = 58 (cerium) (Figs. 5, 6).
The range from 3 Li to 58 Ce forms a first class of CEPXS interest, which we denote Class-I. From 59 Pr, we observe the emergence of a systematic MC-BFP deviation triggered around D max , the depth of the maximum dose. Overall, the latter error both increases and expands with Z, independent of the beam energy. For the 15 MeV beam, at D max , it goes from 2.268% to 2.296% , 3.080% , 3.381% , 3.392% respectively, for 60 Nd, 61 Pm, 62 Sm, 63 Eu Table 1. Percentage of voxels satisfying the AAPM 2% criterion for a relative BFP-MC error below 2% . Listed values include buildup, interfaces and straggling regions. From top to bottom, the three showed sections are: the W-benchmark, the T-benchmark (and its materials) and the IORT-benchmark (and its materials).  (Fig. 7). The increase of the MC-BFP gap in spatial expansion around D max is ∼ 4% (of the total voxels) per increase of Z from 59 Pr to 92 U. For this second class of interest from 59 Pr to 92 U, which we denote Class-II, the failure concerns only the region around D max . The buildup and the tail region are spared. A third CEPXS class of interest from 93 Np to 99 Es, denoted Class-III, arises for which the BFP-MC divergence is total and considerable (Fig. 7). There are exceptions to this classification. First, BFP transport is ambiguous for gaseous slabs, e.g., 1 H, 2 He, 9 F and 54 Xe. Understanding the origin of this deterministic transport anomaly in gaseous media is beyond the scope of this study. The literature does not address this topic either. Second, 32 Ge is an exception to the performance of Class-I materials (Fig. 5), i.e., the failure is total for this slab ( ∼ 15% at 15 MeV around D max ). Third, the excellent performance for both 85 At and 87 Fr slabs is unexpected for Class-II atoms (Fig. 7). We therefore retain that the CEPXS feed functions in ELECTR are acceptable from Z = 1 to Z = 58 , 85 At and 87 Fr, except for 32 Ge, 3 Li,4 Be and purely gaseous media. Moreover, we observe that decreasing the densities of Class-II materials from those of NIST to 1.0g/cm 3 systematically reduces the BFP-MC error. This is because of an ensuing decrease in the interaction rate. For example, decreasing the 92 U slab's density from 18.94g/cm 3 to 1.0g/cm 3 reduces the maximum error from 4.62% to 2.50% at D max . The error around D max is systematically corrected too. Therefore, Class-II materials can meet the 2% criterion at the limit of a density reduction. However, the same test will not work for Class-III materials. For example, decreasing the 98 Cf slab's density (from 15.1 to 1.0g/cm 3 ) reduces the error from 22% to just 12% at D max . While for Class-I and Class-II materials, the BFP-MC error is quite independent of the beam's energy (Figs. 5, 6 and 7), Class-III materials' errors are sensitive to the incident beam's energy. Figure 8a depicts BFP-MC error averaged over all voxels for Class-I and Class-II atoms as a function of the beam energy. However, the average error remains unrepresentative, given that in 95% of cases, it remains below the 2% criterion. On the other hand, Fig. 8a confirms that the MC-BFP deviation is (i) not very sensitive to the beam energy; and (ii) highly sensitive to Z. The distinction of Class effect is better identified in Fig. 8b, which shows the percentage of voxels with a BFP-MC relative error below 2.0% as a function of Z for different beams. The transition from Class-I ( 100% of voxels) to Class-II ( ∼30-80% of voxels) is obvious around 59 Pr. The large increase in the error for Class-III and the decrease in the number of efficient voxels to the limit of ∼ 20% is also evident, in Fig. 8b, around 93 Np. The exceptions, in particular gas slabs and 32 Ge from Class-I, 85 At and 87 Fr from Class-II, are noticeable in Fig. 8b. Figure 9 shows energy deposition profiles of the BFP Dragon-5 solver vs. Lockwood's experimental data 68 and Egs-nrc. Two Dragon-5 computational schemes are examined; the first with the legacy CEPXS-BFP multigroup library 109 and the second with the CEPXS mode in ELECTR. Lockwood's calorimetric measurements are limited to 1 MeV unidirectional beams. The analyzed benchmarks are characterized by heterogeneous, high and medium Z materials. Four observations can be made. First, the CEPXS feed functions in ELECTR are equivalent to the CEPXS-BFP code. This provides strong evidence that ELECTR-calculated cross sections are reliable. Second, the Dragon-5-Lockwood difference is lower than the experimental data precision ( ∼ 2% ) for 98% of the voxels (all benchmarks included). Third, the Dragon-Egs-nrc agreement is lower than the precision on the cross sections ( ∼ 2% ) for 100% of the voxels. Fourth, there is an unexpected agreement for the uranium slab (Fig. 9d) given that this atom is of Class-II (Fig. 7). The uranium agreement may be explained by (1) the beam's energy (in Fig. 9) which is very low for a clean classification; or (2) the G4EmLivermore library (in Figs. 5, 6, 7) could be problematic for this atom. Addressing this question would require an investigation beyond the scope of the current work.

Conclusion
In this paper, we presented and validated ELECTR, a state-of-the-art multigroup electron cross section generation module in NJOY. ELECTR produces a GENDF-formatted library containing: multigroup total catastrophic electroatomic cross sections, relaxation cascade production, anisotropic Legendre components of within-group elastic and group-to-group catastrophic inelastic scattering matrices, radiative and collisional soft stopping powers, multigroup Fokker-Planck momentum transfer coefficients, energy and charge deposition cross sections. This validation was limited to the CEPXS mode in ELECTR. Therefore, the interactions handled include: Møller inelastic scattering, Mott elastic scattering, bremsstrahlung, but also fluorescence, Auger, Coster-Kronig and super Coster-Kronig electron productions. The NJOY MATXSR post-processing module was upgraded to format the produced GENDF library in a MATXS format. The Dragon-5 Boltzmann-Fokker-Planck (BFP) solver was also upgraded to access the MATXS library and solve the BFP equation. The validation of the Njoy-Dragon chain was proposed against experimental calorimetric data and Egs-nrc on Lockwood's benchmarks, but also against Geant-4 on typical radiation oncology benchmarks. The limits of the CEPXS mode have been reported by irradiating slabs from Z = 1 (hydrogen) up to Z = 99 (einsteinium), from 1 to 20 MeV, covering the entire radiotherapy and radiosurgery clinical beams' spectra. For typical high and medium Z Lockwood's slabs, Njoy-Dragon dose differences, with respect to calorimetric measurements, were below experimental precision for 99% of voxels. Whereas for 100% of voxels, Njoy-Dragon vs. Egs-nrc dose discrepancies were within the order of the precision on the cross sections. For homogeneous water slabs, for beams from 1 to 20 MeV, all voxels combined, the average BFP-MC error was below 0.52% . Of particular note, 100% of water voxels satisfied the 2% criterion. For the thorax benchmark, all beams combined, the average BFP-MC deviation was 0.65% , 0.75% , 0.52% and 0.56% , respectively, in the first water slab, bone, lung and the last slab of water. Precisely, 99.2% and 98.4% of the thorax voxels satisfied the 2% criterion, respectively, at 9 and 15 MeV. A systematic BFP-MC deviation was observed at the bone interface discontinuity point. This error is due to a crossing boundary problem on the Geant-4 side which can be resolved by forcing the step size of the electron's track. For a typical Intra-Operative Radiotherapy (IORT) benchmark, all beams combined, the average BFP-MC relative error was 0.51% , 0.96% , 0.00% and 0.00% , respectively, in tumor tissue, aluminum, steel and breast tissue. For all beams, 99.7% of IORT tumor and breast voxels satisfied the AAPM 2% criterion. Finally, for the high-heterogeneity patient-like benchmark, all beams combined, the average BFP-MC relative error was about 0.56% , while, in the worst-case, 97.0% of adipose, muscle, bone and lung voxels satisfied the 2% criterion. By irradiating homogeneous slabs from Z = 1 (hydrogen) to Z = 99 (einsteinium), we have determined the limits of the CEPXS mode in ELECTR. Apart from gas exception, there is an excellent BFP-MC agreement (below 2% ) from 1 H to 58 Ce. Then, from 59 Pr, a systematic deviation appears around D max (the maximum dose deposition depth). This error increases and expands with Z, independently of the beam's energy. The error's gain in spatial expansion, around D max , is of ∼ 4% (of total voxels) per 1-unit of Z increase from 59 Pr to 92 U. The error goes from 2.27% for 60 Nd to 4.62% for 92 U. From 93 Np, a last category is observed for which the CEPXS data are no longer reliable. All the benchmarks studied are in 1D to exclude any deterministic bias in the validation of the multigroup library.
We believe that additional precision and quality assurance can be achieved with ENDF feed functions, i.e., EEDL-2017, EPDL-2017 and EADL-2017 data. One natural next step is therefore to study the potential performance of the ELECTR module in the ENDF mode. Such a library is expected to correct reported anomalies, cover conventional radiotherapy and radiosurgery treatment planning but also the very high-energy radiotherapy modality for deep-seated tumors. The ENDF mode will make it possible to use the same data on both sides (BFP and MC), which opens the door to sensitivity analyzes and the explanation of the origin of the differences. This library is expected to cover a wide range of energies from 100eV to 100GeV.

Data availibility
The datasets used and analyzed during the current study are available from the corresponding author on reasonable request. The DRAGON-5 code and the MATXS-formatted libraries are available under the GNU Lesser General Public License (LGPL). The ENDF mode in the ELECTR code will also be made available under open source license by 2024.