Simulation of a beam rotation system for a spallation source

With a nominal beam power of nearly 1 MWon target, the Swiss Spallation Neutron Source (SINQ), ranks among the world ’ s most powerful spallation neutron sources. The proton beam transport to the SINQ target is carriedoutexclusivelybymeansoflinearmagneticelements.InthetransportlinetoSINQthebeamisscattered in two meson production targets and as a consequence, at the SINQ target entrance the beam shape can be described by Gaussian distributions in transverse x and y directions with tails cut short by collimators. This leads to a highly nonuniform power distribution inside the SINQ target, giving rise to thermal and mechanical stresses. In view of a future proton beam intensity upgrade, the possibility of homogenizing the beam distribution by means of a fast beam rotation system is currently under investigation. Important aspects which need to be studied are the impact of a rotating proton beam on the resulting neutron spectra, spatial flux distributionsandadditional — previouslynotpresent — protonlossescausingunwantedactivationofaccelerator components. Hence a new source description method was developed for the radiation transport code MCNPX. This new feature makes direct use of the results from the proton beam optics code TURTLE. Its advantage to existing MCNPX source options is that all phase space information and correlations of each primary beam particle computed with TURTLE are preserved and transferred to MCNPX. Simulations of the different beam distributions together with their consequences in terms of neutron production are presented in this publication. Additionally, a detailed description of the coupling method between TURTLE and MCNPX is provided.


I. INTRODUCTION
The PSI high intensity proton accelerator (HIPA) generates a continuous wave 1.3 MW beam. Protons are provided by an electron cyclotron resonance source and brought to 590 MeV energy by an accelerator chain composed of a Cockcroft-Walton generator followed by an injector and a ring cyclotron [1]. The 1.3 MW beam is transported through a 60 m long channel equipped with two graphite meson production targets, target M and E which are 5 and 40 mm thick, respectively. The rms beam parameters after Target E are (for the nomenclature see Table V): x ¼ 0.9 mm, θ ¼ 6.2 mrad, y ¼ 0.8 mm, ϕ ¼ 6.8 mrad, dp=p ¼ 0.1%. The 575 MeV beam leaving target E (around 10% is lost here) is reshaped by a collimator system (losing again approximately 20% of the beam) and the remaining about 70% of the 1.3 MW beam is transported to the SINQ target through a 55 m long beam line [2].
In Fig. 1 the layout of the SINQ beam line is illustrated starting with the first bending magnet located downstream of target E (the collimator dimensions are given in Table I). Figure 2 depicts the geometrical representation of the SINQ target region used for the particle transport simulation with MCNPX 2.7.0; the model consists of the last section of the beam line where three water cooled copper collimators (KHN31-33) shield the rim of the target beam entrance window and, at the same time, prevent activation of the beam line by backscattered neutrons. The shapes of the collimators resemble the proton beam envelope to inhibit a focused proton beam at the SINQ target.
The visual monitor and sieve (VIMOS) [3] is one of the three independent safety systems in place to prevent an insufficiently scattered beam reaching SINQ. All three systems are able to switch off the beam within 100 ms when 10% of the protons bypass Target E. VIMOS is the closest one to the target, sitting on top of the collimator KHN33. A mirror-camera system monitors optically the glowing of a nearly monocrystalline tungsten mesh.
The SINQ target is composed of a bundle of Zircaloy cladded Lead rods, colloquially called Cannelloni target. The target is actively cooled with heavy water and the neutron production is enhanced by a ring-shaped lead reflector [4]. The target is submerged in a cylindrical tank of 2 m diameter filled with heavy water which contains two moderator inserts, an ambient temperature light water volume of about 1 liter and a cold source of approximately 20 liter D 2 at about 25 K with a reentrant hole.
Since the proton beam transport is carried out only by means of linear magnetic elements, the beam shape at the SINQ target can be described by Gaussian distributions in both x and y coordinates (the z-coordinate is the beam direction), whereas tails are cut by collimators installed downstream of target E and upstream of the SINQ target. This highly inhomogeneous proton current distribution coupled with the large heat load causes thermal and mechanical stresses that potentially could lead, on the long term, to damage of the SINQ target, in particular in view of the 1.8 MW upgrade planned for HIPA. This issue triggered a systematic study aiming for the design of a magnetic proton beam rotation system that would produce a more uniform proton beam distribution at the SINQ target, thus reducing the maximal heat load and its large gradients. At the same time, the change in the shape of the proton beam profile has an impact on the neutron flux which has to be quantified. Furthermore, a suitable rotation frequency has to be estimated which gives acceptable variation in the time-dependent neutron flux from the neutron instruments point of view.

II. DESIGN ASPECTS
In a proton beam transport line, two different techniques can be employed in order to produce a more uniform beam distribution. One can either apply a static approach, where beam fringes are folded by means of nonlinear elements (octupoles and/or dodecapoles), or a dynamic one, where beam rotation or rastering are provided by high frequency AC or sawtooth wave dipoles. The static method has several drawbacks which make it very difficult to apply in the SINQ facility. In its latest version, the SINQ target possesses a 120 mm diameter circular cross section, thus the folded beam-which would have a rectangular footprint with sharply peaked edges-could damage the target rim or cause large activation. Moreover, the installation of octupoles and/or dodecapoles magnets in the already almost fully packed proton beam line would require an extensive modification of the present setup. On the other hand, the two dipole magnets needed for the beam rotation  could be integrated in one relatively short beam line segment that could be placed in the only free section located 24 m upstream of the SINQ target ( Fig. 1). Such a solution would also permit an easy and reasonably safe preliminary test. An issue related to this method is that a relatively large rotation frequency is needed in order to allow for a quasi time independent spatial neutron distribution and, at the same time, avoid local overheating of the target. For the SINQ spallation source, a rotation frequency of several 100 Hz has to be foreseen (see the timedependent calculations in Sec. IV); yet more investigations are necessary in order to fix this value.

III. BEAM OPTICS CALCULATIONS
A. Analytical approach The beam rotation concept foresees to bring the beam on a circular path at the SINQ target. This can be achieved using two magnetic deflectors in the beam line, a horizontal and a vertical one which are modulated according to a cosðωtÞ and b sinðωtÞ, respectively. To achieve a circle at the target location, the modulation amplitudes a and b must be chosen according to the beam optics at the deflectors location. In polar coordinates the resulting averaged beam profile can be computed by a convolution of the Gaussian beam distribution on a circular path: The coordinate system is given in Fig. 3. Here r 0 is the rotation radius, σ the rms beam width and I 0 is a Bessel function. The function is normalized to give the total current I after integration over the radius r. Ideally the outer boundary of the beam should not be larger after applying the rotation, thus the beam size must be decreased. The decrease can be achieved either by emittance reduction or by stronger focusing in the transport channel.
The ideal situation is illustrated for a few cases in Fig. 4: the two parameters r 0 and σ are varied. A reduction of the peak beam density by a factor 2 can be achieved. For the practical application at the HIPA the beam width can not be easily reduced, so we rely on collimation of the beam tails to avoid higher losses with the rotated beam.

B. Simulations
The simulation of the proton beam optics for the rotating beam was carried out using the standard tools Graphic TRANSPORT and TURTLE [5]. Using TRANSPORT, a modified beam optics setup for SINQ was simulated with the aim to reduce the beam footprint radius at the SINQ target by 20 to 30% to preserve the area hit by the proton beam if it is rotated. Figure 5 shows the x (bottom) and y (top) 2σ SINQ beam envelope for the current proton beam line setup (green, tagged as reference case in the following), as well as the modified one (black, tagged as rotating case  density normalized radius r 0 = 0.00, σ = 1.00, peak 100% r 0 = 1.00, σ = 0.77, peak 73% r 0 = 1.12, σ = 0.71, peak 60% in the remainder of the paper) for which beam rotation will be applied, i.e., in the picture the black curve is a snapshot of the rotated beam with kick angle zero. In Fig. 5, the label "ROT" represents the location of the rotation system while "TSNQ" is the position of the SINQ target. The beam rotation was then introduced in TURTLE, a ray-tracing code to calculate particle trajectories, the resulting beam distributions and losses. The beam rotation was simulated by adding a 4.5 mrad kick to the beam centroid at the location foreseen for the rotation dipoles ("ROT" in Fig. 5). The amplitude of the kick was chosen in order to match the SINQ target cross section. This procedure was repeated 36 times with an azimuthal step size of 10°. A crucial boundary condition for any new beam optics is to minimize additional losses in the collimators with respect to present beam line setup.
The TURTLE simulation results are presented in Fig. 6 where the reference (top) as well as the rotating (bottom) x and y beam distributions are displayed. The simulations start at the upstream end of target E and the initial number of protons is the same in both cases (3.6 million). The timeaveraged spatial distribution of the rotating beam is more uniform in both coordinates than the current beam line setup. As a consequence the peak current density is reduced by a factor 2 from 23.7 to 11.5 μA=ðcm 2 mAÞ.
Compared to the reference beam, the time-averaged footprint of the rotating beam is smaller in x and larger in y, becoming almost circular and matching therefore much better the geometry of the SINQ target. After Target E, the beam pipe is in the shadow of the collimators KHE0-3, thus the losses are negligible for both the reference and the rotated case (less than 0.005%).
These results are partially obtained at the cost of increasing the losses at the collimator system KHN31-33. While in the reference case losses in this region are very low (0.009% corresponding to 0.11 kW in KHN33 and negligible in KNH31-32), in the rotating case beam losses as high as 0.318% (4.13 kW), 0.080% (1.05 kW) and 0.147% (1.91 kW) would take place in KHN31, KHN32 and KHN33, respectively. The aforementioned losses were calculated with TURTLE. Slightly different values are observed using MCNPX (see Table II). The loss calculations in the collimators KHN31-33 were also confirmed with OPAL [6].
Using TURTLE, the phase space information of each proton is written to a text file at the entrance of collimator KHN31. These data are then directly used in MCNPX to determine the effect of the rotating beam on the neutron production compared to the reference case, and for loss calculations.

A. MCNPX model
A detailed MCNPX ( [7]) geometry model of the SINQ including the target, the thermal scatterer and the cold source, the beam ports, all submerged into the D 2 O moderator tank is used for the neutronics calculations. The model was extended to include the collimators KHN31, KHN32, KHN33, their surrounding cooling system and the beam position monitor VIMOS (Fig. 2). A text file written by TURTLE containing all phase space information for each beam proton at the entrance of collimator KHN31 is used for the source description in MCNPX. Since heavy shielding surrounds the collimators, tracking of the protons which pass through any of the collimators is omitted; these protons are counted as losses.
A detailed description of the new source definition option of MCNPX can be found in the appendix of this paper.

B. Comparison with experimental results
In 2006, a liquid metal target-MEGAPIE-was successfully operated at SINQ, PSI [8]. From the simulation point of view, the only difference to Fig. 2 is that the target contains liquid lead-bismuth eutectic instead of Zircalloy rods. The target received an integrated current of 2.8 Ah of 575 MeV protons. The PIE included the γ-mapping of the cusp-shaped AlMg3 beam entrance window (BEW). Any part of the BEW could be precisely positioned in front of a 2 × 2 mm pin-hole=collimator at which exit-at a distance of 1.41 m-a Germanium detector was located. This allowed the measurement of the 22 Na activity, which is proportional to the proton fluence since the 22 Na production is dominated by the reaction 27 Alðp; XÞ 22 Na. For more details see [9]. Based on the measurements a detailed source description could be constructed which is then used on the SDEF card in MCNPX.
Additionally, the beam optics calculations are simulated with TURTLE based on the magnet settings at the time of MEGAPIE operation. The output file from TURTLE is fed to a subsequent MCNPX calculation. It has to be noted that the experimental determination of the proton beam profile revealed a displacement of the proton beams centroid by several mm from the beam line axis of symmetry. To account for this shift the TURTLE beam optics simulation was transformed accordingly. The proton beam profile upstream of the BEW is shown in Fig. 7.
The difference between the measured, integrated proton beam profile onto the MEGAPIE target and the simulated profile, using the new coupling method of MCNPX and TURTLE, is depicted in Fig. 8. The central part of the beam distribution is well described by the simulation. For the major axis of the elliptical beam footprint deviations above 25% appear beyond a radius of 4 cm, while the minor axis first shows an overprediction of the simulated values at radii larger than 2 cm. The most outer region of the beam footprint is poorly described by the calculation; there is a clear underprediction with respect to the measurement.

C. Comparison of the reference and the rotating beam
As already pointed out, the use of two dipole magnets for the beam rotation is a relatively simple way to produce a more uniform beam distribution. A major criterion for a realization of the beam rotation is that the performance of the neutron source (basically the neutron flux) does not deteriorate significantly. Furthermore, the losses along the proton beam line have to be small.

Beam losses
The losses calculated with MCNPX and with TURTLE are given separately for the three collimators in Table II. The differences come mainly from the different physics treatment in the collimators. This has been verified with calculations where physics was turned off, i.e., particles hitting the collimator surface were instantly absorbed. Comparing these results the agreement was better than 6%.

Beam profile at the target
The time-averaged proton beam profiles upstream of the BEW of the target are plotted in Fig. 9. The peak current density is about 50% smaller in the rotating case compared to the reference case. Furthermore, the distribution is more uniform and its envelope changed from an elliptical to an almost circular shape.

Power density distribution
Due to the more uniform proton beam, the peak power density-in the central cannelloni rods along the targetsignificantly decreases (see Fig. 10); in the rotating case the maximum power density is about 3 cm from the center (see Fig. 11). Furthermore, the power density profile, perpendicular to the beam direction, flattens as depicted in Fig. 12. Due to the more uniform power distribution the maximum temperature will be smaller in the rotating case. The resulting thermomechanical stresses have to be calculated in the future and compared to the reference case. However, as temperature induced stresses are depending on temperature gradients mainly, a reduction of the stresses in the central part of the rods can be expected.

Time-averaged cold and thermal neutron flux
It has been shown that the more uniform distribution of the rotating beam profile is advantageous with respect to the power density distribution in the target. Now a significant decline of the neutron production in the target and the decrease of the neutron flux at the instruments has to be ruled out. Since most of the neutron scattering instruments are positioned to view the cold source and the thermal scatterer via beam ports reaching into the D 2 O tank of SINQ, the time-averaged neutron flux is calculated for different wavelength bands in two control volumes as shown in Fig. 13. The fluxes from these two control volumes are representative for the cold source (control volume 3017) and the thermal beam lines (control volume 21). The neutron flux results are summarized in Table III. The differences of the neutron fluxes for the two proton beam options are a few percents. The rotating case gives 4-5% higher values for the thermal flux, whereas it gives 2-3% less for the cold flux compared to the reference case. Further optimization of the magnet settings and collimator alignments will potentially result in even lower differences.
It should be noted that a reduction in the peak power density does not translate into the reduction of the timeaveraged neutron flux, because the envelope of the beam profile for the rotating case is still inside of the target.

Time-dependent neutron flux for the rotating case
The results from the previous subsection showed that the time-averaged neutron flux remains the same when a beam rotation is applied. However, due to the rotation, the centroid of the proton beam will periodically be closer or further away from the moderator inserts and the neutron beam ports. Therefore, a variation of the neutron flux in time at these positions is expected. The beam rotation in TURTLE is simulated as a series of 36 cases; each case is created by adding a 4.5 mrad directional kick to the beam for every 10°segment. Each case of the 36 TURTLE output files are run separately with MCNPX, recording the timedependent neutron flux in the two control volumes, see Fig. 13. The results are then superimposed by taking into account the rotation frequency of the proton beam. Figure 14 shows two characteristic proton beam profile cases with a phase shift of 180°. The distance between the  With an infinitely high rotation frequency, the neutron flux would be constant, somewhere between the values from the two cases in Fig. 14. This is the time-averaged neutron flux, calculated in Sec. IV C 4, and depicted as a black line in Fig. 16. On the opposite, with a very small (close to zero) rotation frequency the neutron flux would oscillate between the two cases. Thus, with a finite frequency, the neutron flux varies between these limits. Figure 16 shows the oscillation of the cold neutron flux (between 2.5 and 5.0 Å) in the liquid D 2 source as a function of rotation frequencies. Table IV summarizes the results for the cold source and the control volume in front of the thermal neutron beam port.

V. CONCLUSION
A preliminary design study of a proton beam rotation system for the SINQ spallation source at PSI was presented. A new method was developed which couples the beam optics code TURTLE with the Monte Carlo particle transport code MCNPX, thus allowing the simulation of more realistic beam profiles. The MPI-version of the tool was successfully tested and used for the calculations presented in this publication.
Results show that the proton beam rotation would bring substantial benefits in terms of the power distribution in the SINQ target without affecting the neutron spectrum and yield. The beam peak current density on the target can be reduced by a factor 2. The magnitude of the proton beam rotation frequency is still a matter of investigation. This issue is of particular importance because of its implication on the technology needed for the rotation magnets. Furthermore, the final value will certainly depend on the requirements of the neutron scattering instruments: assuming that a 10% change in the neutron flux is acceptable a rotation frequency of several 100 Hz will be necessary.
Future studies should cover abnormal operation conditions, e.g., an unexpected failure of the magnets responsible for the rotation. In this case the beam would hit the SINQ target with a reduced beam footprint. Although this accident scenario is similar to an overfocusing event for which safety systems are already in place, it has to be checked in the future. Moreover, studies for an improvement of the proton beam losses in the collimators upstream of the SINQ target are under way.
Based on the above findings it has to be decided in the future if such a system is feasible to be built and implemented at SINQ.

ACKNOWLEDGMENTS
We would like to thank A. Adelmann and V. Rizzoglio for the calculations with OPAL.

APPENDIX: NEW SOURCE ROUTINE FOR MCNPX
TURTLE offers the possibility to create ASCII files containing the phase space information of the transported, and not lost, particles at any location of interest. Although MCNPX provides a great variety of source definition options, prior analysis with an averaging of the phase space information evaluated with beam transport simulation packages is necessary. As a consequence of the averaging some phase space information is always lost. Therefore, a new source definition feature has been implemented in MCNPX which uses the output file(s) of TURTLE without any averaging processes and information loss.
TURTLE only provides phase space information about spatial coordinates and their derivatives from the coordinates perpendicular to the main beam direction (see Table V). Hence, in addition the user has to provide the starting coordinate in beam direction (the z-coordinate is considered to be the direction of beam propagation).

Name
Dimension Description x mm x coordinate of the proton θ mrad Angle between the z-axis and the projection of the momentum vector on the xz-plane y mm y coordinate of the proton ϕ mrad Angle between the z-axis and the projection of the momentum vector on the yz-plane p MeV=c Momentum of the proton w Á Á Á Weight of the proton n Á Á Á Counter