Focused VHEE (very high energy electron) beams and dose delivery for radiotherapy applications

This paper presents the first demonstration of deeply penetrating dose delivery using focused very high energy electron (VHEE) beams using quadrupole magnets in Monte Carlo simulations. We show that the focal point is readily modified by linearly changing the quadrupole magnet strength only. We also present a weighted sum of focused electron beams to form a spread-out electron peak (SOEP) over a target region. This has a significantly reduced entrance dose compared to a proton-based spread-out Bragg peak (SOBP). Very high energy electron (VHEE) beams are an exciting prospect in external beam radiotherapy. VHEEs are less sensitive to inhomogeneities than proton and photon beams, have a deep dose reach and could potentially be used to deliver FLASH radiotherapy. The dose distributions of unfocused VHEE produce high entrance and exit doses compared to other radiotherapy modalities unless focusing is employed, and in this case the entrance dose is considerably improved over existing radiations. We have investigated both symmetric and asymmetric focusing as well as focusing with a range of beam energies.

Radiotherapy is used in the treatment of around half of all cancer cases worldwide 1 , and 27% of cases in the UK 2 . The aim of radiotherapy is to deliver a lethal dose to the tumour whilst sparing healthy tissue as much as possible; since it was first invented the field has advanced greatly and highly conformal photon radiotherapy techniques such as intensity modulated radiotherapy (IMRT) and volumetric modulated arc therapy (VMAT) are now routinely used in cancer treatment 3 . In addition to conventional radiotherapy, FLASH radiotherapy in which ultra-high dose rates are used to deliver radiotherapy in a manner that reduces normal tissue toxicity is currently being investigated 4 with the first patient treated in Lausanne in 2019 5 . Proton radiotherapy is becoming more widely used due to the highly favourable dose distribution in the form of the Bragg peak, with minimal dose delivered downstream from the target in the beam direction. Proton radiotherapy has the disadvantage of requiring large halls and expensive accelerators and beam transport systems 6 ; as such there is current interest in other forms of radiotherapy that have similarly beneficial dose distributions but are less space-intensive.
For deep-seated tumours, the radiation must penetrate 15-30 cm into the body. Different types of radiotherapy have different integrated dose distributions, as shown in Fig. 1. Conventional electron radiotherapy peaks at or before 5 cm into the body before sharply falling off, making it unsuitable for the treatment of deepseated tumours. Clinical photon beams also sharply peak close to the surface of the patient before exponentially falling off, resulting in a high entrance dose, though this can be decreased by combining multiple beams from different angles. Proton beams can penetrate far into the patient and deposit most of their dose in the Bragg peak before falling off sharply, and have a relatively low entrance dose, although this does increase up to ≥40% of the maximum dose for higher energy ( ≥180 MeV) proton beams. The entire tumour region can be covered using multiple proton beams of different energies weighted by intensity using for example the method shown in the 1984 paper by Bortfeld and Schlegel 7 in order to produce a spread-out Bragg Peak, although this results in a significantly higher entrance dose (60-80%). Proton beams are highly affected by inhomogeneities 8 , with the result that a small change in composition in the body pre-treatment such as the presence of an air gap can result in the proton Bragg peak being deposited outside of the tumour. Motion due to for example breathing can also cause the Bragg peak to be deposited in the wrong region. Proton radiotherapy therefore requires robustness measures which increase the dose close to the target. Carbon-ion beams have similar dose distributions and characteristics to proton beams, but with reduced penumbral spreading 9  www.nature.com/scientificreports/ peak, and also a smaller percentage energy spread 10 resulting in a sharper Bragg peak. Carbon-ion beams also have a higher relative biological effectiveness (RBE) thereby causing more damage to the tumour, but with the disadvantage of requiring higher energy machines, which are more expensive and require significantly more space than proton radiotherapy machines 11 . Very high energy electron (VHEE) beams (with energies in the region 50-250 MeV) can penetrate far into the body. The benefit of this is that VHEE beams could be used to treat deep-seated tumours and have been shown to be much less sensitive to inhomogeneities than proton beams 12 . While no suitable machines explicitly developed for VHEE radiotherapy currently exist, advancements in accelerator technology have enabled linear accelerators (linacs) to robustly produce 100 MeV/m gradients [13][14][15] . This provides some confidence that similar robust structures could be used for medical linacs 16 . This combined with the fact that VHEE beams could be delivered rapidly makes VHEE an attractive candidate for FLASH radiotherapy [16][17][18] . There are currently plans for two VHEE FLASH radiotherapy treatment machines, the CERN/CHUV collaboration 16 and PHASER at Stanford University 19 . Studies by Bazalova Carter et al. 20 have shown that VHEE plans are consistently similar or better at treating certain tumours than the photon VMAT plans. However, VHEE dose distributions have a high entrance and exit dose. In order to circumvent this, the possibility of using quadrupoles to focus the VHEE beams before they enter the phantom has been investigated.
Studies by McAuley et al. 21 using proton beams and by Kokurewicz et al. 22 and Lagzda et al. 23 using VHEE beams have shown experimentally that beam focusing at shallow depths of approximately 5 cm (comparable to normal electron beam radiotherapy) in a water phantom is possible using quadrupole magnets. The aim of this study is to investigate the feasibility of using realistic quadrupole magnets to focus deeply penetrating VHEE beams by performing Monte-Carlo simulations to depths ≥ 15 cm. The resultant dose distributions inside a water phantom will be investigated for 250 MeV symmetrically and asymmetrically focused beams, with the asymmetrically focused beams also being investigated for 100 MeV, 150 MeV and 200 MeV. The effect of changing the final magnet strength on the on-axis dose distributions and beam σ will also be investigated, both using Monte Carlo and analytical calculations. The dose distribution results will be used to produce a spread-out electron peak over a target region in the water phantom by summing the dose distributions with weighting factors. This will be compared to a proton spread-out Bragg peak (SOBP) over the same region.

Results
The Elegant particle physics code 24 is a rapid, well-tested code used routinely in accelerator physics for machine design of particle orbits in a vacuum. The Elegant code is however unable to predict the effects of scattering in air and water; as such Monte-Carlo simulations in TOPAS 25 are also required to model focused VHEE beams in a water phantom. The Elegant code was used to rapidly find the initial optimal number, position and strengths of quadrupoles to focus 250 MeV beams symmetrically and asymmetrically in a vacuum. This optimisation required multiple iterations that would take many hours using the TOPAS code, compared to seconds in the Elegant code. Thus the initial simulations in Elegant formed a starting basis in the optimised focusing positions in TOPAS simulations. The method and quadrupole strengths used are described in the "Methods" section.
The normalised on-axis dose distributions for the symmetrically and asymmetrically focused 250 MeV VHEE are shown in Fig. 2, with unfocused 250 MeV VHEE for comparison. In both focusing cases the quadrupole strengths found by Elegant were meant to focus at 20 cm into the phantom; the resultant TOPAS simulations show that the actual position of the focal point in the water simulation varies due to scattering and the quadrupole set up used. The dose distributions in the transverse and longitudinal planes are shown in Fig. 3. These show that for approximately symmetrically focused beams that are less sharply focused in the on-axis dose distributions, the beam enters the phantom at a similar size in both transverse planes, and is not sharply focused at the focal point. For the asymmetrically focused beams, a stronger focusing effect occurs at the focal point but the effect is different in the x and y planes. This is due to the nature of the quadrupoles-each quadrupole produces focusing in one transverse plane whilst defocusing in the other transverse plane. The benefit to this type of focusing is www.nature.com/scientificreports/ that it spreads out the entrance dose over a wider area, resulting in a lower entrance dose at any point over the surface of the phantom, as shown in Fig. 3 and Table 1, whilst the exit doses remain similar for both focused and unfocused VHEE. The cross-sections of the dose at the focal point are also shown in Fig. 3. The corresponding Gaussian beam σ x and σ y at the focal point are shown in Table 1. For proton radiotherapy, the beam σ at the Bragg peak is typically of the order of 10 mm at high energies ( ≥100 MeV) 26 , meaning that the 250 MeV VHEE beam is smaller in both planes than a typical proton beam at the same depth. The Gaussian beam σ x and σ y at different depths through the water phantom is shown in Supplementary Fig. 1 in Supplementary Materials for symmetrically, asymmetrically and unfocused VHEE. This shows that the symmetric focusing has a minor effect on the shape of the dose distribution in the water phantom, whilst asymmetrically focused VHEE markedly changes this, with noticeable differences between the shape of the beam in the x-and y-planes. While pencil-beam spot scanning using FLASH radiotherapy is under investigation, it is currently clear the FLASH effects are apparent for ultra-high dose rates and also very short treatment times (in many < 0.1 s) 27 . For this reason, having a dose distribution that spreads out the entrance dose over a wide region and targets a volume in the patient could be a successful modality for FLASH VHEE radiotherapy. Therefore, due to the potentially beneficial nature of the dose distributions for asymmetrically focused VHEE, other beam energies were also investigated for the asymmetric case, with the aim of investigating the beam energies required for effective VHEE focusing to treat deep-seated tumours. The quadrupole strengths found for asymmetrically focused 250 MeV were adjusted for 100 MeV, 150 MeV and 200 MeV VHEE beams using the method shown in the "Methods" section, with the quadrupole positions remaining the same. The quadrupole strengths found are shown in Supplementary Table 2 in Supplementary Materials.
The on-axis dose distributions for the four different beam energies are shown in Fig. 4. The result is that higher energy beams require higher quadrupole strengths to focus them, and as such the on-axis dose distributions are more sharply focused for the higher energy VHEE than the lower energies. The entrance dose for the four beam energies are shown in Table 2, showing a reduction in dose by energy of 20% when increasing the beam energy from 100 to 250 MeV. This result also shows that the beam energy used changes the position of the maximum dose; this is due to the increasing penetration power and reduced scattering of VHEE with beam energy. The σ x and σ y of the dose distributions throughout the water phantom for the different beam energies are shown in Supplementary Fig. 1 in Supplementary Materials. This shows the differences in focusing in the x-and y-planes, as well as the reduced scattering with energy of the beam. The dose distributions in the transverse and longitudinal planes for the four different beam energies are shown in Fig. 5. In all cases, the dose distribution is different in each transverse plane as expected, with the strongest focusing achieved for 250 MeV. The cross sections of the dose distributions at the focal point are also different for each beam energy, as shown in Fig. 5. It is clear that the higher energy beams produced smaller beam σ s at the focal point; this effect can be attributed to the reduced penumbral spreading with energy for VHEE, as well as the higher focusing strengths used. The Gaussian σ x and σ y at the focal point for each of the four beam energies are shown in Table 2. This shows that by changing the beam energy from 100 to 250 MeV, the beam σ s at the focal point reduce by approximately 55%. Both the 200 MeV and 250 MeV VHEE have beam σ s approximately the same size or smaller as a typical 160 MeV proton Bragg peak in water 26 , with the 150 MeV beam σ x being only 3 mm larger. This implies that any of these beam energies could be clinically useful with regards to beam size at the target.
Scanning could be simplified if only one quadrupole was the variable, therefore the effect of changing the final quadrupole for each beam energy were investigated first analytically and then in terms of resultant dose distributions. The Twiss parameters 28 (see "Methods" section) at the start of the phantom were found for each of the beam energies, and used to predict where the focal point should be. The resultant focal positions in the xand y-planes were found for a range of different final magnet strengths for each of the beam energies. These were compared to the position of the maximum dose in the water phantom. The results are shown in Supplementary Table 3 in Supplementary Materials. It was found that the maximum of the dose distribution occurred at the same position in both the x-and y-planes, which was different to each of the focal positions predicted by Twiss  www.nature.com/scientificreports/ parameter analysis. To investigate this, TOPAS simulations were run for 250 MeV VHEE with the same magnet parameters but in a vacuum, and again in air, in order to compare the results for each, shown in Supplementary  Fig. 2 in Supplementary Materials. It was found that there was a difference of 1-2 cm in the position of the focused dose between the analytical predictions from the Twiss parameters from Elegant 24 and the full water phantom simulation. This is likely due to scattering of the beam in the water phantom producing a maximum in the dose at a depth less than that in a vacuum. A small difference in the results for the analytical prediction and the vacuum results is likely due to space charge effects or differences in the code of how the quadrupole magnets are simulated. The TOPAS simulations though, through many PDD simulations compared to experiments, have been show to be robust and accurate 29 .
Changing the strength of the final quadrupole has an effect on the width and position of the resultant on-axis dose distribution, as well as the shape of the dose distribution in the transverse planes. Figure 6 shows how the position of the maximum dose (i.e. the focal point) changes by changing the final quadrupole strength for each of the 4 beam energies. This could be used to cover a tumour with any length inside the patient, just by changing the final quadrupole strength. The relationship between the position of maximum dose (i.e. the focal point), ẑ (m) and final quadrupole strength, g 4 (T/m) was found to be approximately linear. Accounting for the beam energy, E (MeV) as well, this fit is extended by including quadratic terms, with the general formula where a 1 , a 2 , a 3 , a 4 and a 5 are constants obtained from a least square error fit. This readily enables us to predict intermediate energy behaviour. The fits for position of maximum dose as a function of final magnet strength and energy for the plots shown in Fig. 6 are given in Table 3. The R 2 value for the goodness of fit is 0.999. It is straightforward to convert ẑ to cm as all that is required is to multiply by 100. Changing the final quadrupole strength also affects the shape of the beam. Knowledge of how the beam shape changes with quadrupole strength could be used to shape the beam to the tumour at a particular location. The relationship between the beam σ x in m at the focal point as a function of the final quadrupole strength, g 4 (T/m) was found to be approximately quadratic, with the general formula where a 1 , a 2 and a 3 are constants obtained from a least square error fit, and with a similar relation for σ y . The fits for the beam σ x and σ y at the focal point as a function of the final quadrupole strength for the four beam energies are given in Table 4. For σ x and σ y in cm we simply multiply by 100. The fits produced for position of focal point as a function of final magnet strength, g 4 and σ x,y as a function of g 4 were used to produce a fit of σ x,y as a function of focal point position, ẑ . The results of these fits are shown in Fig. 6c, d. These results could be used to predict the beam size for a required depth of maximum dose. The half-width-half-maxima (HWHM)  www.nature.com/scientificreports/ of the on-axis dose distributions for each of the beam energies and magnet strengths are shown in Fig. 6. This shows that focusing is most efficient close to the surface, due to the higher focusing strengths used and reduced beam scattering. As a first step towards showing how focused VHEE could be used to cover a tumour region, a spread-out electron peak (SOEP) has been produced for 250 MeV VHEE, shown in Fig. 7 by changing the final quadrupole strength. Mathematica 30 was used to optimise for flatness over a region of 12-17 cm into the water phantom, i.e. corresponding to a target size of 5 cm. This was done by taking the on-axis dose for each of the different magnet strengths and multiplying by a variable weighting factor, then minimising the deviation from unity by summing the weighted doses over the region by changing each of the weights. This method mimics how multiple proton beams are combined using weighting factors to produce a spread-out Bragg peak (SOBP), with the idea that multiple SOEPs or SOBPs could be combined to treat a 3D tumour region. The resultant on-axis SOEP and SOBP covering the target region are shown in Fig. 7, with parameters shown in Table 5. The transverse dose distributions at the entrance, start of the target, end of the target and exit of the phantom are also shown in Fig. 7, showing the shape and intensity of the SOEP dose and SOBP dose through the water phantom. The colour of the dose distributions are all normalised to the maximum dose in the SOEP or SOBP. This shows the difference between the two types of beams; the proton SOBP has a higher entrance dose, but no exit dose and the dose is contained within approximately 1 cm transversely. The SOEP has a lower entrance dose, but the dose is spread out over a wider region at the entrance, and there is some exit dose. In both cases the target region is

Discussion
Focusing of deeply penetrating VHEE beams using realistic quadrupole magnets in Monte Carlo simulations has been shown here. The on-axis dose distribution and transverse dose distributions show a significant reduction in entrance dose and similar exit dose compared to unfocused VHEE, and the use of the quadrupole magnets means the position of the maximum dose can be chosen by simply changing the quadrupole strengths. The reduction in entrance dose is highest for 250 MeV VHEE focused asymmetrically, although the reduction in     www.nature.com/scientificreports/ entrance dose is only 20% from 100 to 250 MeV, compared to 47% for unfocused 250 MeV VHEE and 100 MeV focused VHEE. The quadrupole magnets shape both the on-axis dose distributions and the shape of the beam in both transverse planes. By careful selection of the magnet strengths, the shape of the beam at the focal point can be pre-selected with varying beam σ s in the x-plane and y-plane. The effect of changing the beam energy from 100 to 250 MeV results in the beam σ s at the focal point being reduced by approximately 55%, from σ x = 1.84 cm and σ y = 1.15 cm for 100 MeV to σ x = 0.84 cm and σ y = 0.53 cm for 250 MeV. For the 150 MeV, 200 MeV and 250 MeV VHEE, the beam σ s at the focal point are all approximately the same order of magnitude as the corresponding proton Bragg peak in air at the same depth (15 cm). This implies any of these beam energies could be appropriate for focused VHEE radiotherapy. Previous Monte Carlo studies in focused VHEE radiotherapy have mostly involved theoretical focusing by pre-selecting focusing strengths in the Monte Carlo simulation 31 with no attachment to realistic magnetic properties. This study goes a step further towards realistic focusing of VHEE radiotherapy by investigating the focusing strengths and number of realistic quadrupole magnets required to produce focusing in both transverse planes in a semi-analytical way, in contrast to previous studies in proton focusing (see for instance 21 )which take a simple FODO lattice using available magnets, rather than optimised magnet parameters. Two experimental studies 22,23 have been completed so far focusing VHEE beams at the CLEAR facility, with focal depths of approximately 5 cm into the water phantom and focusing in only one plane. This study shows for the first time that by changing the strength of only the final quadrupole magnet, the position of the maximum dose can be changed, enabling dose coverage of a whole region (> 20 cm) by changing the final quadrupole strength. This relationship is linear over the region of the water phantom for each of the four beam energies, so this could potentially be done in a continuous manner, which could be useful in a clinical setting when dose must be delivered in a short time frame. Predicting the position of the focal point from beam energy and final quadrupole strength is possible by extending this fit to include quadratic terms. This study shows for the first time that combining multiple focused VHEE beams produced simply by changing the final quadrupole strength allows for a spread-out electron peak (SOEP) to be produced over a target region on-axis, with a lower entrance dose than the corresponding on-axis spread-out Bragg peak (SOBP), a low exit dose, and with a similar flatness to the SOBP. This is a stepping stone towards 3D focused VHEE treatment planning, with the idea that multiple SOEPs could be combined to create a 3D target dose distribution.
The results in this study were produced by using four quadrupoles for asymmetrically focused beams and six quadrupoles for the symmetrically focused beams; this is not optimal for a clinical setting. In addition the quadrupole magnets are placed very close to the phantom and the beam size within the quadrupoles is not constrained.
Generally, clinically useful beams have been produced in Monte Carlo simulations in this study, and the magnet parameters required for an experiment to show realistic VHEE beam focusing have been found. Such an experiment would be possible at the CLEAR facility at CERN, with minimal changes required to account for the achievable beam parameters. This experiment would validate the transverse and longitudinal dose distributions of focused VHEE beams, the production of a spread-out electron peak by changing only the final quadrupole strength, as well as the potential for continuous scanning of the final quadrupole strength to cover a range of depths in a water phantom. Future research is required to translate this into a clinical setting, in particular how to reduce the size of the apparatus, and the possibility of combining the set-up with a dipole magnet steering system to produce a VHEE treatment planning system.

Conclusions
Monte Carlo simulations performed in this study show that it is possible to focus very high energy electron (VHEE) beams to produce high dose volumes inside a water phantom at depths required for deep-seated tumours, and that it is possible to change the shape of the dose distribution by changing only the quadrupole configuration. Symmetric and asymmetric focusing was investigated for 250 MeV VHEE, with lower entrance dose for the asymmetric case. The effect of beam energy on focusing was investigated for 100 MeV, 150 MeV, 200 MeV and 250 MeV asymmetrically focused VHEE, with the beam σ s all approximately 1 cm at the focal point for 150 MeV, 200 MeV and 250 MeV VHEE. Increasing the beam energy from 100 to 250 MeV decreased the entrance dose by a further 20%, from 45 to 25%. The on-axis dose distributions show reduced entrance doses and similar exit doses compared with unfocused VHEE beams. This study also shows that the position of the maximum dose can be changed by changing only the final quadrupole strength for each of the four beam energies, allowing for the possibility of producing spread-out electron peaks (SOEPs) similar to spread-out Bragg peaks (SOBPs) simply by altering the final magnet strength.The relationship between maximum dose position and final quadrupole strength is approximately linear, allowing for this to be performed potentially in a continuous manner during dose delivery. The position of the maximum dose can also be predicted by a straightforward two parameter quadratic fit involving beam energy and final quadrupole strength. Further research is required into how to translate focused VHEE into a clinical setting, including reducing the number of quadrupole magnets required and combining the magnets with a VHEE treatment planning system. Evaluation of beam focusing position using Twiss parameters. Particle motion in an accelerator is described by the Hill's equation 28 , where x is the particle's displacement from the central orbit, s, and x ′ = dx ds is the angle with respect to the central orbit, and K(s) is the focusing strength at a particular point along the central orbit (in practice K is piecewise linear). The general solutions to Hill's equation are where ǫ is the beam emittance, β(s) is the amplitude modulation due to changing focusing strength which combined with α and γ = 1+α 2 (s) β(s) make up the Twiss parameters. Also, ψ(s) is the phase advance. All particles in the beam travel along their own trajectories within an ellipse in phase space (x, x ′ ) which can be shown 28 to be of area

Methods
x ′ (s) = − α ǫ/β(s) cos ψ(s) − ǫ/β(s) sin ψ(s) www.nature.com/scientificreports/ Choosing the particle on the largest phase ellipse within a beam means that all the particles will stay within that ellipse, allowing for behaviour of the whole beam to be described by a single particle. This is shown in Fig. 9. The particle's transverse position and angle vary throughout the machine due to differences in focusing strength, and can be represented by a matrix. The modulations due to different types of magnets can be represented as a Transport Matrix, M, between two points, ( As in the literature 33 , the matrix element for a focusing quadrupole M QF of length L Q and focusing strength K is given by similarly where = √ |K|L Q . The matrix element for a defocusing quadrupole, M QD of length L Q and strength −K is similarly given by Note that the magnetic field gradient, g (T/m) of a quadrupole magnet can be found from the focusing strength, where E is the beam energy (MeV). Note here the speed of light has been approximated as 3 × 10 8 m/s and the beam travelling with velocity v is assumed to be sufficiently relativistic that β = v/c ≈ 1 . For a drift space, the matrix element, M d for a drift of length L (m) is given by The transfer matrix for the full four quadrupole setup can be found by multiplying the matrices, with a similar transfer matrix for the six quadrupole case.
The Twiss parameters α, β, γ (= (1 + α 2 )/β) characterise the beam electron trajectories. They are related to the magnet strengths and drift spaces the particles traverse. Our aim is to focus the beam at a particular depth, L, within the phantom. Equation (7) allows us to evaluate the beam position ( x, x ′ ) in phase space after progressing through all of the quadrupole magnets and a water phantom-all encompassed within the matrix M. We assume the water phantom is modelled by a drift space of length L p (no scattering included) and we can now propagate the Twiss parameters through all of the magnets using Eq. (12) and obtain their values at the phantom by (9) g ≈ EK/300   The positions of the focal points from the Twiss parameters for the 250 MeV beam were found using the Elegant code for the different focusing strengths shown in Supplementary Table 1 (asymmetric beam) in Supplementary Materials to produce the Twiss parameters at the entrance to the phantom, and then using the equation shown in the Twiss parameter analysis section, the corresponding focal lengths in the water phantom were found for each of the x-and y-planes. The focal lengths in the x-plane for each of the different final magnet strengths were then plotted with the positions of the maxima of the dose distributions from TOPAS 25 . These were found by running the simulation as described in the Monte Carlo simulation details section, but with the whole simulation in only a vacuum for the vacuum case, and only air for the air case, each with beam energy of 250 MeV and the focusing strengths shown in Supplementary Table 1 (asymmetric beam) in Supplementary Materials. 501 × 501 × 501 voxels were used here and the central voxel was taken to find the on-axis dose distribution to increase the precision of the maximum dose position.
SOEP. The spread-out electron peak corresponds to delivering a flat distribution of dose over the target region. We were able to achieve this by summing a series of n weighted electron dose profiles with various magnet strengths. Here we only needed to vary the final quadrupole magnet strength, g 4 . Over the target region, we minimised the following equation where D i (z) corresponds to the dose profile for a particular magnet strength, a i to a weighting factor (varied), and z min,max to the location of the minimum and maximum of the target region. We used the Mathematica code 30 to perform this minimisation. We characterised the flatness over the target region with where σ D is the standard deviation of the dose and D is the mean dose. For the 5 cm target, at a depth of 12-17 cm, we found no more than four focused beams were necessary to produce an on-axis field flatness of better than 99%. It is worth noting that a similar optimised SOEP can be achieved by varying the energy (rather than the magnet strength). Indeed, for protons we obtained a spread-out Bragg peak (SOBP) by varying the energy-but in this case we needed more than twice as many separate energies. The transverse dose at various depths through the SOEP and SOBP were produced by summing the weighted doses of each of the contributing beams and plotting the dose deposited in a 30 cm × 30 cm × 0.297 cm slice through the water phantom at depths from 0 to 30 cm in increments of 3 cm, as well as the dose at 17 cm (i.e. the edge of the target region).

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.