OPTIMIZATION PROBLEM ON MODELLING OF BORON NEUTRON CAPTURE THERAPY SYSTEM

The paper discusses dynamics of charged particles and neutrons in boron neutron capture therapy system (BNCT) as well as geometrical and physical optimization of BNCT system elements. Our choice is BNCT system with linear accelerator. BNCT track includes ion injector, RFQ, DTL, neutron-producing target and neutron moderator which provides an exit (last collimator) flux of epithermal neutrons satisfied to International Atomic Energy Agency (IAEA) requirements. The following software tools IBSimu, LIDOS, COMSOL Multiphysics and PHITS were used for modelling BNCT system


Introduction
The status of accelerator-based worldwide BNCT projects at the present time (end of 2019) is given in paper [Hashimoto, 2014]. The list includes fourteen BNCT installations. They are distinguished by composition of facilities and main output parameters of accelerators (energy and current) : linear RF accelerators (6 units), electrostatic accelerators of different constructions (5 units) and cyclotrons (3 units). Two types of neutron-producing targets are used in these installations: Beryllium and Lithium ones. All the types of accelerators which are used for BNCT may be separated in two of groups "high energy -small current" and "small energy -high current" in according to the main output parameters. It is hard to accelerate large current in cyclotrons, until now cyclotron current limit is 2 mA. So, contemporary BNCT system with cyclotron has on surface neutron-producing target 30 MeV energy and 1 mA current. It is impossible to obtain "high beam energy" in electrostatic accelerators (until now it is limited by 2-2,5 MeV) if we want to have reasonable accelerator dimensions. Linear RF accelerators doesn't have such limitations. Linacs can provide all range of energies and currents which need BNCT system. However linear RF accelerator works usually in pulse mode and needs more power consumption than cyclotron or electrostatic accelerator. That is way different accelerator laboratories in the world choose different accelerator types for BNCT. Thus, choice of optimal accelerator type depends on several criterions. But main criterion is a necessity to obtain neutron flux of epithermal neutrons at the moderator exit equal or more than 10 9 n/cm 2 /s under minimal level (≤ 5%) thermal and fast neutron fluxes. It is not simple problem because neutrons are not governed by electromagnetic fields (practically). This problem must be solved by optimization of physical and geometrical parameters of moderator.

Ion Beam Transport and Neutron Producing
We will consider the reaction 9 Be(p,n) with incident proton beam of energy10 MeV as source of neutron flux for BNCT system. According to the results of the paper [Kiyanagi, 2019], 10 MeV is an optimal energy because it makes possible to obtain necessary proton beam power (30 kW or more) with enough small current (about 5 mA) and satisfy to criterion of limitation of neutron induced radioactivity from neutron producing target. It is Figure 2. Parameters of injection system supposed that an accelerator will operate in CW mode with the value of current 5 mA. Small current makes it possible to use resonators with working frequencies of Pdiapason (400-450 MHz). Choice of power amplifier for RF system working in such diapason was discussed in paper [Svistunov et al, 2019]. Block diagram of BNCT system is given in the Figure 1 Injection system was modelled by the software tool IB-Simu [Kalvas, 2013]. IBSimu is a software for modelling ion optics, plasma extraction, ion beam transport taking into account space charge effects. In contrast to other codes IBSimu indirectly takes into account plasma state in ion source when ion beam take starts. Initial data including physical plasma parameters, electrode potentials and parameters of iterative process during dynamics simulation are given in the Table 1 and Figure 2. Table 1. Initial data for dynamic's simulation Current density particle flux J = 150 A/m 2 Plasma potential ϕ p = 13 V Magnetic field suppression threshold potential ϕ sup = 10 kV Transverse electron temperature T e = 2 eV Transverse ion temperature T t = 1 eV Electron to ion beam ratio eir = 10 Initial drift energy extracted particles E 0 = 2 eV Temperature of thermal compensating ions T p = 0.5 eV Ratio of fast compensating ions to total negative charges R f = 0.5 Induction of magnetic field taking off electrons from the beam B = 0.025 T Radius of aperture of plasma electrode r 0 = 0.008m Point of start on Z axes z s = −0.003 m Grid step size h = 0.0004 m Number of particles in mesh of net n perh = 50 Accuracy requirement for particle iterator ε part = 10 −4 -10 −5 Space charge under relaxation factor α = 0.5 The authors of [Kalvas, 2013;Kalvas et al, 2014] selected RADIS ion source and extraction system as prototype. Configurations of electrodes were preserved but potentials were changed. Considered injector includes plasma source and Low Energy Beam Transport (LEBT) system. Plasma flux is extracted from source by puller. Electrons are taking off from the beam by transverse permanent magnetic field on dump electrode. LEBT system forms and accelerates H − beam up to 50 keV before injecting into RFQ. Layout of LEBT system and beam trajectories are presented in the Figure 2. Beam phase volumes at planes xx and yy are presented in the Figure 3 and Figure 4 together with distribution of particle density. Output transverse rms emittances (grey ellipses in the pictures) which are used as characteristics of transverse dynamics in transport codes equal 6.59·10 −6 m·rad and 6.57·10 −6 m·rad at phase planes xx and yy accordingly. Twiss parameters α, β.γ are (0.3538, 0.059, 70.973) and (-2.033, 0.0359, 143.135) at planes xx and yy accordingly. . Emittance plot at plane xx , z=0.125 m Such a beam may be inscribed in acceptance of RFQ with working frequency 401 or 402.5 MHz. For example, simulation of SNS injector by IBSimu determined input norm.rms.emittance of SNS RFQ equal to 5 · 10 −8 m · rad with the current 20 mA and energy 65 keV [Kalvas, 2013;Kalvas et al, 2012]. Working frequency of SNS RFQ is 402.5 MHz. In our case working frequency of resonators is 401 MHz. Design and calculate parameters RFQ are given in the Table 2.     [Bondarev et al, 2001] was used. For dynamics optimization, additionally numerical optimization methods described in paper [Altsybeyev et al, 2018] were used. Results of dynamics simulation are presented in the Figure 5 and Figure 6. After RFQ ions are accelerated by DTL structure from 2.5 up to 10 MeV. DTL includes two of RF resonators placed in common volume (tank). Resonators are partitioned by thin dividing wall with a hole for beam passing. Beam focusing in DTL is accomplished by the permanent quadrupole magnets placed inside the drift tubes. Optimal magnetic lattice structure FFODDO was determined in [Skudnova, 2019]. Other design and calculated parameters of DTL are presented in Table 3. Beam phase volumes at the DTL exit are presented in Figure 7 (gray ellipses on the picture are RMS emittances).  Proton beam (10 MeV, 4.9 mA, ∆E/E=±1% for 90% particles) leaves DTL and impinges Beryllium target. Be-target has 1 mm thickness, copper can 1cm and cooling water channel 0.3 cm. These parameters are chosen to decrease the probability of blistering. As a result of reaction 9 Be(p,n) 1.25·10 12 n/cm 2 s neutron flux is appeared with energy distribution in diapason 0.1-8 MeV. Further motion, beam forming and neutron flux parameters at the exit of last collimator are determined by the Beam Shape Assembly (BSA).

Problem of BSA and Moderator Optimization
Neutron-producing target is a part of acceleration tract, but structurally it is a part of BSA. BSA includes neutron target, filter, moderator, thermal neutron shield, proton shield and collimators. BSA scheme is given in Figure 8. The main aim of BSA is to load out from neutron beam thermal and fast neutrons and to form flux of epithermal neutrons with the required density at the aperture of the last BSA collimator. Because more than ten different materials are proposed for using in moderator composition, optimal choice is quite difficult According to the present literary sources [Atomic Energy Agency, 2001;Koay et al, 2017] Fe, Al, AlF 3 , MgF 2 , CaF 2 are most effective materials for moderation. For simulation neutron beam dynamics after Beryllium target program product PHITS version 3.10 was used. This program was working out by leading Japan laboratories JAEA, RIST, KEK and others [Skudnova, 2019]. Program PHITS is intended to simulation of different tasks when beam dynamics modelling requires Monte Carlo methods. In case of modelling neutron interaction with matherials of moderator, PHITS useS library JENDL4.0. Moderation system consists of a few layers of different materials placed in sequence after neutron-producing target. Also BSA includes reflector and shield materials placed around THE moderator. Reflector reverses part of neutrons which are flying perpendicularly to longitudinal moderator axis to inside moderator. The shield protects the outer space from of neutron radiation to outer space. The main question in moderator design is a choice of materials for moderation system (it is discrete set of parameters). Base task of moderation system is to obtain at the moderator's exit intensive flux of epithermal neutrons with energy distribution from E min =0.5 eV up to E max =10 keV under minimal level of fast neutrons with energy greater than E max and thermal neutrons with the energy from 0 up to E min . Levels 0.5 eV and 10 keV were determined by IAEA-ECDOC-1223 recommended values of neutron beam characteristics in 2001. They are presented below in the Table 4. But these requirements were formulated by IAEA for BNCT, which worked with reactor neutrons and possible needs to be inspected. Let k be the number of different material layers of moderation system and x i is the thickness of i-layer. Then, vector x = (x 1 , x 2 , ..., x k ) describes moderation system as sequence material layers of determined thickness. Let Φ epi (x), Φ th (x) and Φ fast (x) be epithermal, thermal and fast neutrons fluxes on moderator exit correspondingly, that are where function n(E), E ∈ (0, +∞) is density of neutron flux distribution by energy.
On the bases of physical reasons it is possible to impose on x i restrictions in form α i ≤ x i ≤ β i , so set the admissible solutions One may consider three of quantitative criterions during moderator design: where Instead of multicriterion optimization problem (1) we may find minimum of the function where I is ideal meaning of f 1 (x) ≤ I. In our case, it is epithermal flux equal to 10 9 neutron/cm 2 /s. In spite of condition (1), strict minimization of functions f 2 (x) and f 3 (x) is not necessary. For example, if each out of functions f 2 (x), f 3 (x) is less then 0.05 simultaneously and the value of f (x) is near I for several solutions x ∈ X, one may to choose among these solutions an optimal one in a Pareto sense. Search of optimal solution is produced with help of inertial gradient method where δ n and γ n are chose on results of iteration process. Search is continued as long as obtained solution x ∈ X will satisfy the problem's conditions. Kit of layers x and hence functions f 1 (x), f 2 (x), f 3 (x) are determined by the results of dynamic simulation. Partial derivations of function f (x) are calculated by numerical derivation. Initial approach is determined as fixed choice moderator materials with equal thicknesses of layers. Fe, Al, ALF 3 , LiF, Cd materials were considered for moderator in different layer sequences. To withdraw thermal neutrons out of the beam, one considered LiF and Cd. Use of LiF gives better result. Fe was better then Al in terms of stopping the fast neutrons. In Figure  9 and Figure 10 intermediate and final results of investigations are shown. Best sequence of layers is given in Figure 10: Be -0.1 cm, Fe -15 cm, AlF 3 -30 cm, LiF -1 cm. In Figure 9, flux of epithermal neutrons is larger than in Figure 10, but the fluxes of fast and thermal neutrons exceed the admissible level too much. Sequence of layers for this case is Be -0.1 cm, Fe -22 cm, AlF 3 -34 cm.

Conclusion
Some comments may be added to the obtained results.
1. This result is close to the goal, but the required epithermal flux was not achieved. Epithermal neutrons flux may be increased if proton energy on Beryllium target is increased up to 12 MeV. But it contradicts with the condition 10 MeV as optimal energy from point of view radiation limitation. One may note that broadening of epithermal energy interval from 0.5 eV up to 30 keV gives in the BNCT system epithermal neutron flux 0.87·10 9 n/cm 2 s.
2. On the other hand, obtained result corresponds to the local extremum only. Possibilities of optimization of BSA and total BNCT tract are not exhausted. Data in Table 4 emphasize on the hardness of the optimization problem. 3. In these calculations we supposed CW current regime for resonators of BNCT system. It corresponds to light current loading but serious thermal loading for resonators. On the other hand, regime of patient's irradiation is not continuous and this must simplify watercooling system design.
4. So, in principle it was shown that it is possible to design BNCT system based on compact linear accelerator with P-frequency diapason.