Numerical optimization of longitudinal collimator geometry for novel x-ray field

Objective. A novel x-ray field produced by an ultrathin conical target is described in the literature. However, the optimal design for an associated collimator remains ambiguous. Current optimization methods using Monte Carlo calculations restrict the efficiency and robustness of the design process. A more generic optimization method that reduces parameter constraints while minimizing computational load is necessary. A numerical method for optimizing the longitudinal collimator hole geometry for a cylindrically-symmetrical x-ray tube is demonstrated and compared to Monte Carlo calculations. Approach. The x-ray phase space was modelled as a four-dimensional histogram differential in photon initial position, final position, and photon energy. The collimator was modeled as a stack of thin washers with varying inner radii. Simulated annealing was employed to optimize this set of inner radii according to various objective functions calculated on the photon flux at a specified plane. Main results. The analytical transport model used for optimization was validated against Monte Carlo calculations using Geant4 via its wrapper, TOPAS. Optimized collimators and the resulting photon flux profiles are presented for three focal spot sizes and five positions of the source. Optimizations were performed with multiple objective functions based on various weightings of precision, intensity, and field flatness metrics. Finally, a select set of these optimized collimators, plus a parallel-hole collimator for comparison, were modeled in TOPAS. The evolution of the radiation field profiles are presented for various positions of the source for each collimator. Significance. This novel optimization strategy proved consistent and robust across the range of x-ray tube settings regardless of the optimization starting point. Common collimator geometries were re-derived using this algorithm while simultaneously optimizing geometry-specific parameters. The advantages of this strategy over iterative Monte Carlo-based techniques, including computational efficiency, radiation source-specificity, and solution flexibility, make it a desirable optimization method for complex irradiation geometries.


Introduction
As the field of radiotherapy explores new realms of treatment, including ultra-high dose rate therapy, microbeam therapy, spatially fractionated radiotherapy, and highly conformal stereotactic radiosurgery (Zhang and Mayr 2023), there is a greater demand for increasingly complex radiation sources and beam shaping devices in the clinical and preclinical spaces (Insley et al 2024).Traditional wisdom on x-ray source design has been upended by the advent of exciting technologies such as field emission carbon nanotube sources, x-ray tube miniaturization, complex target design (Insley et al 2024), x-ray focusing lenses (Bartkoski et al 2021), and beyond.With the implementation of novel sources into clinical practice, it has become clear that if the x-ray source deviates from conventional design principles, then the construction of other necessary accessory devices, including dosimetry tools, positioning apparatuses, and collimating structures, should be reassessed (Yin et al 2010).
Toward such an effort, a proof-of-concept analysis has been described in the literature for a novel conical target geometry optimized for high intensity and high directionality (Insley et al 2024).This nontraditional target features a 1.5 mm thick tungsten layer printed on a 1 mm thick conical diamond substrate with a 10°cone semi-angle and a 6 mm base radius.A conical fluence of 200 keV electrons converge and bombard the surface of the tungsten-coated diamond cone producing an x-ray field that is higher intensity and more forward-directed compared to other x-ray tubes of this size and energy.To modulate the x-ray field's intensity and focal spot size, the annular surface of carbon nanotubes used to generate the electron beam is segmented into concentric rings.Independent electrical switching of each ring adjusts the size of the electron bombardment area on the target (Insley et al 2024).Some of the main components of this tube are displayed in the upper half of figure 1 along with exemplary electron trajectories shown by semitransparent yellow arrows.
The high dose rates and directional coherence of the source make it optimal for low-penetration external beam therapeutic applications, such as preclinical small animal irradiation and intraoperative radiotherapy, where collimation of the x-ray field to precise field sizes is necessary for conformal dose delivery.However, the optimal collimator geometry that preserves both intensity and precision is made ambiguous by the complexities of the x-ray target geometry.The conical shape of the target and spread of the electron beam creates a highly distributed radiation source-essentially, this design approximates a line source while preferentially attenuating photons emitted at oblique angles along the tungsten-diamond structure.While pinhole collimators are regularly used in nuclear medicine to focus distributed sources to smaller regions (Bushberg et al 2012), the selfattenuation of the target and inherent heterogeneities in the underlying bremsstrahlung shape function (Omar et al 2020) differentiate this source's radiation field from that of a radionuclide point source distribution.
Attempts at optimizing collimator geometry found in the literature resort to iterative Monte Carlo fluence and dose calculations when simple analytical models of the radiation production and transport are not available (Lu et al 2012, Wang et al 2021, Clements et al 2022).The immense computational burden of Monte Carlo forces the designer to first assume a collimator shape (pinhole, parallel-hole, etc) and then optimize that shape's limited set of parameters (radius, thickness, etc) (Smith et al 2003, Lu et al 2012, Wang et al 2021, Clements et al 2022).While this type of simplification may be sufficient for simpler source geometries, it reduces the overall robustness of the design process for more complex targets.The aim of this work is to develop an algorithm that i and inner radius, r .
i The system is set with a constant distance between the phase space plane and detector plane, z 0 .
avoids the repeated use of Monte Carlo during optimization.Thus, fewer constraints can be placed on the collimator geometry, allowing for a more generic optimum to be found from a broader search of the parameter space.The following sections first derive a general analytical model for the photon flux resulting from a radiation source coupled to a cylindrically symmetrical collimator.We then describe a simulated annealing approach to optimize the collimator geometry.Finally, this approach is applied to the novel radiation source described previously (Insley et al 2024).

Theory
For this work, optimizations were focused on stereotactic irradiations akin to CyberKnife and GammaKnife treatments.These devices use circular collimators, multiple dose focal spots, and converging beam geometry to paint arbitrary dose distributions (Rad andOrton 2006, van de Water et al 2011).Adhering to those same principles, this collimator was restricted to cylindrical symmetry (i.e.circular cross sections), and the longitudinal shape of the collimator was optimized.The geometry of this formulation is shown in figure 1.This algorithm consists of two separate models: one describing the radiation emitted from the x-ray tube (the source model) and the other describing the attenuation through the collimator (the collimator model).Both models reference two planes: the phase space plane (the surface flush to the aperture of the x-ray tube) and the detector plane (some downstream parallel surface that may represent an isocentric plane in a treatment context).Coordinates in the phase space plane are described using ( ) ¢ ¢ x y , ,and coordinates in the detector plane are described using ( ) x y , .det det These planes are separated by a constant, arbitrary distance, z .0 The source model describes the radiation field using the following differential flux density equation: is the mean number of photons, N, originating with energy in the interval ( ) from a differential area segment in the phase space plane located at ( ) ¢ ¢ x y , with area ¢ ¢ x y d d and passing through a differential area segment in the detector plane located at ( ) x , 0 det with area dx dy .
det det y det does not appear explicitly in this formula because we have assumed cylindrical symmetry of the radiation field, thus allowing us to just calculate the photon flux along one arbitrary axis (e.g. the x det axis).With the implicit inclusion of y d , det the units of ( ) 1 2 The collimator model is an extension of an idea posed by Wang et al in which a preliminary optimization of several stacked conical collimator segments was performed (Wang et al 2021).To make the geometry even more generic and flexible, the collimator is split longitudinally into very thin slices with central holes, resembling a stack of washers.Each washer is assumed to extend well beyond the range of the radiation field and is thus characterized by two parameters: z, the position of the washer above the detector plane, and r, the inner radius.For a stack of washers that comprise the whole collimator, let z i and r i be the ith washer's position and inner radius, respectively.
The path between a photon's starting position in the phase space plane, ( ) ¢ ¢ x y z , , , 0 and its final position in the detector plane, ( ) x y , ,0 , det det is uniquely defined.The attenuation of that track through the ith washer is given by the Beer-Lambert exponential attenuation law: where N 0 is the number of photons leaving from the differential area segment at ( ) ¢ ¢ x y z , , , 0 N is the number of photons arriving at the differential area segment at ( ) x , 0, 0 , det m is the material-and energy-dependent linear attenuation coefficient, and t i is the photon's pathlength through the ith washer.This pathlength is formulated as the following piecewise function to model the material discontinuity between the central hole and the body of the washer: where z d is the thickness of the collimator slice. () q ¢ ¢ x y x , , det describes the angle of incidence of the photon at the plane of the i th washer, which causes a lengthening of the path through the washer.The cosine of this angle is given by ( ( )) ( ) x y x cos , , . 2 Finally, ¢ r i is the radius at which a given photon track crosses through the plane of the ith collimator slice, calculated using the photon's initial and final positions as In the calculation of t , i ¢ r i is compared to r i to see if the photon passes through the given washer's central hole or body.These two models can be combined for an arbitrary number of collimator slices to calculate the photon flux at the detector plane: x det is the flux of photons in units of inverse area that are passing through the detector plane as a function of distance from the central axis.The flux is assumed to be cylindrically symmetrical, so this profile is equivalent for any axis in the detector plane.If we modulate this equation to , cos det det then ( ) Y x det gives the photon fluence.Given the small field sizes and small phase space plane studied in this work, a small angle approximation is used (i.e. ( ) q » cos 1) such that the planar fluence and fluence are approximately equal.This allows for a significant reduction in computation time when this system is discretized later.This section will conclude with a formulation of relevant objective functions for the optimization of the collimator geometry.The collimator is optimized via the flux model, ( ) Y x , det which should be modulated towards high intensity and high precision.Collimators are typically designed for various field sizes, so precision should be defined within the context of a desired field size.With that in mind, the primary objective function used for optimization is given by: This objective function is composed of two main factors: ( ) Y 0 and ¥ .

F r F
col ( ) Y 0 is the flux at the center of the detector plane, which is the maximum of the flux profile.

F r F
col is defined using the function, det which represents the number of photons passing through an infinitesimal strip of thickness y d det out to a distance R along the x det axis.If r col measures the radius of the desired field size (e.g. 5 mm), then col is the fraction of photons within the desired field size.Accordingly, field size is defined at the detector plane.Maximizing the product of these two factors, ( ) gives a simultaneous optimum for both intensity and precision.
Alternatively, these factors could be of different relative orders, as in the following functions: G 2 and G 3 represent intensity-weighted and precision-weighted objectives, respectively.We also investigated a fourth objective function that prioritizes flatness within the desired field size: det col Together with the inverted precision fraction, this objective function should be minimized to create a flat field of width r .col Finally, these objective functions are combined into a single function to represent all three prioritiesintensity, precision, and flatness.This objective is given by: In the second term, the constant b is selected by the user to control the relative magnitude of the first and second terms, which in turn controls the prioritization between intensity/precision and flatness.For this study, b is held constant at 1.Each of the objective functions in equations (4)-( 8) are summarized in table 1 for reference.

Materials and methods
In this section, the practical methods for calculating the phase space model, ( ) , , det and performing the optimization will be explained.First, it is necessary to discretize each of our models, i.e.

N x y x E det det
The resolutions and ranges for each parameter are given in table 2. ) in the target and window.Additionally, a directional filter was used to perform Russian Roulette on particles not headed towards the phase space scoring surface, and particles with weight above 0.001 (i.e.'fat particles' that would disrupt the statistics of the final result (Shanmugasundaram and Chandrasekaran 2018)) were excluded.These physics settings were validated using published measurement data on the Sculptura electronic brachytherapy system (Badali et al 2019), which shares many similarities with the x-ray target described by Insley et al (2024), including its ultrathin tungsten-coated diamond target design and energy range.Using the same physics settings described above, simulations were constructed from the Sculptura specifications provided by Badali et al (2019).An x-ray spectrum (1 keV resolution) was calculated in Geant4 for the 100 kVp setting, which is the energy closest to that used by Insley et al (200 kVp).The first and second half-value layers determined in these simulations were within 3% of the corresponding values measured by Badali et al.Despite not having the full description of the Sculptura source including exact material specifications, this exercise provided evidence that the physics module, range cut, and variance reduction parameters were appropriate.
The phase space model, ( ) , , det was calculated on a surface flush to the aperture of the x-ray tube.To save CPU time and memory, standard deviation calculations were performed on a relatively small number of electron histories (n = 200 000).Error was then propagated from the uncertainty in each of the histogram's bins to the final flux calculation in equation (3): Here, equation (3) has been converted to its discrete form where k l m , , , and n are indices for the ¢ ¢ x x y , , , det and E bins, respectively.s Y;k is the standard deviation of Y in the k th x det bin.s , , , , , was approximated as the mean variance of all non-zero bins following this short simulation.The mean relative standard deviation from these 200 000 histories was calculated to be 91.18%.Assuming the uncertainty decreases proportional to the inverse square root of the number of histories, the mean relative standard deviation should be approximately 3.45% after 140 million histories.Propagating this mean standard deviation to the calculated flux using equation (9) without any collimation (i.e.= t 0 ) resulted in a maximum uncertainty of 0.2%.Performing these simulations on The University of Texas MD Anderson Cancer Center's Seadragon cluster, which features 80 CPUs per node at a minimum clock speed of 1.9 GHz, allowed for the completion of one billion histories in under 240 h.Thus, calculation of each phase space was performed at one billion initial electron histories for sufficient fidelity.
A few final notes on the generation of this phase space model are included here for completeness.The x-ray target is characterized in the literature for nine different focal spot size settings, ( ) ¢ ¢ f x y x E , , , det was calculated independently for each setting.In addition, ( ) ¢ ¢ f x y x E , , , det was calculated separately for source-to-detector distances of 70, 90, 120, 150, and 200 mm, as stated in table 2. Thus, the full data collection process resulted in 45 different phase space files, each occupying between 2.5 and 3.5 GB.Simulations were performed in parallel on separate nodes of our computing cluster, allowing for the completion of this stage within four weeks from the initial job submissions.
Following the calculation of each phase space, equation (3) was validated against Monte Carlo-derived photon flux calculations.The electron beam simulation was rerun with the same settings, but the phase space scoring surface was replaced by a detector plane 15 cm downstream of the x-ray target aperture.This detector plane scored photon surface current (number of particles crossing per area) in 0.1 mm radial segments.The results of these simulations were compared to the model calculations for three collimator configurations: a flat, 0.1 mm thick attenuator with no central hole, a 1 cm thick parallel-hole collimator with a 1 cm diameter hole, and a 1 cm thick parallel-hole collimator with a 0.5 cm diameter hole.All collimators in this study were assumed to be lead, though this can be easily changed via the attenuation parameter from equation (3).
With this validation method, we were able to assess the change in performance of our model as various memory-saving reductions were made to the phase space.The larger the histogram for ( ) , , det the more computationally intensive the optimization.Thus, we investigated cuts to the energy resolution and phase space plane range.Specifically, we tested re-binning the energy spectra from 1.0 to 10.0 keV resolution.To avoid completely losing attenuation fidelity, the 1.0 keV resolution files were used to calculate effective attenuation coefficients for each bin using the following equation: Here,  n is the new energy index once the phase space is re-binned.Attenuation coefficients for lead were acquired from the National Institute of Standards and Technology for each 1.0 keV energy bin (Hubbell and Seltzer 2022).
In addition, the range of ¢ x and ¢ y were reduced from [−8, 8] to [−6, 8] mm and [0, 6] mm, respectively.This range reduction removed less than 0.1% of the largest focal spot size's total phase space histogram.It is important to note that the ¢ x reduced range is asymmetric and different from the ¢ y reduction. () ¢ ¢ f x y x E , , , det describes photons that originate from anywhere in the phase space plane and pass through the +x det axis between 0 and 25 mm.More photons from the + ¢ x axis will contribute to the flux along the +x det axis than from the -¢ x axis.Thus, more of the phase space from the -¢ x axis can be omitted without significant consequence to the final flux calculation.Conversely, symmetry is still maintained along the ¢ y axis, so the data from the -¢ y axis can be lumped into the + ¢ y axis data without loss of accuracy.A phase space intensity map for the large focal spot size is included in figure 2 to display graphically where these reductions were made.This section will conclude with an overview of the optimization algorithm and relevant hyperparameters.Due to the piecewise formulation of the photon pathlength through a collimator slice in equation (2), discontinuities are introduced into the flux calculation of equation (3) and thus the objective functions of equations (4)-( 8).Hence, the objective functions are not differentiable on their entire domain and any optimization algorithms that rely on the gradient of the objective function cannot be used.Instead, simulated annealing was implemented according to the formulation described by Kirkpatrick et al (1983).As a brief summary, simulated annealing applied to the present optimization problem follows the following process: (1) Initialize the set of washer inner radii, r , i and calculate the objective function.
(2) Iterate through the complete set of washers in a random order.
(3) For each washer, adjust the inner radius by some small, constant value, r d .Whether the inner radius is increased or decreased (i.e. r d ) is determined by a 50/50 probability, unless a decrease would result in the inner radius being less than 0, in this case, increase the inner radius.T n is the simulated 'temperature' as a function of iteration, n.To be clear, n is iterated after each full cycle through all the washers.( ) ( ) a = T n T , n 0 where a and T 0 are hyperparameters that determine the rate at which the system cools.
(5) Repeat steps 2-4 with a different random washer order reach time.For some set number of initial iterations, known as the grace period, track the average DG for unfavorable transitions, D - G .Following this grace period, count the number of consecutive full iterations (i.e.full cycles through all the washers) where the objective function changes by less than bD - G , where b is some small positive number (i. i 1 4 col A single solution was selected from these 3 runs based on the solution with the most optimal objective function value.In a small fraction of runs, the early stages of optimization found a greater optimum than the final solution.These transient solutions can be lost due to the high probability of poor transitions early on.In such scenarios, the state of the optimization was saved and rerun later with a much lower temperature, allowing for the exploration of a different local optimum.
For comparison to the primary objective optimizations, the other objective functions, G , 2 G , 3 G , 4 and G combo were used for optimization of the medium focal spot size for = r 5 col mm only.Just one initial solution was used for these optimizations: = r r .
i col For G 2 and G 3 optimizations, only the order = n 2 was investigated.Finally, each of the optimizations mentioned was repeated for five source positions: = z 7, 0 9, 12, 15, and 20 cm.
In total, 285 optimizations were performed.Subsequently, it became necessary to determine how these collimators could be implemented practically.Ideally, any setting of the focal spot size, source position, and desired field size would be accompanied by a single, specific optimized collimator.If any of those settings changed during a treatment, then the collimator would automatically change its geometry accordingly.This implementation would have to be accomplished by a stack of high-resolution, variable iris collimators that can automatically adjust their inner radius, analogous to the collimator model used in the optimization process.Variable iris collimators have been described by Graves et al in the context of small animal irradiation (Graves et al 2007), but the level of precision needed for our purposes is still not available.
Instead, it is more likely that a few collimators should be selected and interchanged depending on the needs of the application.To investigate the set of optimized collimators further, several configurations were implemented into Monte Carlo simulations.This was done for only three of the nine focal spot sizes: the smallest, medium, and largest.For each of those three focal spot sizes, the optimized collimators from each source position (z 0 ) were modeled in Geant4.Then, photon surface current at the detector plane was calculated in the same manner as the validation simulations.A traditional, 10-tuple phase space was calculated first and used repeatedly in each simulation to save computation time.Only 10 8 electron histories were used for the phase space generation to keep the file size manageable.For each optimized collimator, the radiation source was moved to positions other than the position that the collimator was optimized for.Analyzing a collimator's response to different source positions allowed us to identify the properties of different collimator choices and select a diverse final set.A few of these final collimator choices as well as a parallel hole collimator for comparison will be presented in the results section for each focal spot size and field size.

Results
Figure 3 displays the results of the model validation simulations for the largest focal spot size.Each set of curves represents the photon flux at the detector plane, normalized to the maximum flux of each respective curve, as a function of position along the +x det axis.Each curve uses the same set of lines and markers: blue crosses for the Monte Carlo calculations, red-orange line for the model calculations (equation ( 3)), and yellow line for the model calculations after the phase space had been compressed according to the steps described previously.The three sets of curves correspond to the following setups: a single lead sheet with no central hole (Thin Attenuator), a 1.0 cm diameter parallel-hole collimator, 1.0 cm in thickness (1.0 cm Collimator), and a 0.5 cm diameter parallel-hole collimator, 1.0 cm in thickness (0.5 cm Collimator), as labeled.The same validation data was collected for the medium and small focal spot sizes to assess any change in model performance with x-ray tube setting.The mean percent errors for the large, medium, and small focal spot sizes across all 3 curves are 0.19%, 0.21%, and 0.20%, respectively.The maximum errors are 1.81%, 1.50%, and 7.08%, respectively.These maxima occurred around the smaller central radial bins for the 0.5 cm collimator where Monte Carlo statistical uncertainty was greatest.The small focal spot size's maximum error of 7% is significantly larger than the other focal spot sizes.However, its mean error is comparable to the larger focal spot sizes, and it was calculated that 99.33% of the small focal spot size's points fell within 2% of the Monte Carlo results.This indicates that this single 7% maximum error point is extremely uncommon.
Next, figures 4(a) and (b) give two examples of the G 1 objective function progress throughout optimization.The y-axis displays the value of G 1 in units of flux ( mm 2 ) as a function of iteration number on the x-axis.To be clear, iteration number describes the number of cycles through all the collimator slices.There are 100 collimator slices in the 1 cm thick stack, so each iteration involves 100 collimator adjustments.Figure 4(a) displays a typical objective function progression where the solution at termination was the most optimal solution found.Figure 4(b) gives an example of a rare optimization run in which the most optimal solution was found very early but was lost due to the high temperature.In these runs, the optimization state was saved and rerun with a much lower temperature to investigate a potentially deeper local optimum.
Figures 4(c)-(e) display the progression of each factor of the combination optimization to show the relative magnitudes of each component.While the intensity factor, ( ) Y 0 , is the largest in magnitude, this optimization saw relatively little change in its value compared to precision, i which ranges from 60 to 70 mm.The x-axis gives the radius measured from the x-ray tube's longitudinal axis, r .
i The bottom row contains the resultant photon flux at the detector plane with normalized flux units on the y-axis and distance along the x det axis on the x-axis.The y-axis is normalized by the maximum flux of all the source positions for a given focal spot size.Figure 6 presents the same data for = r 0.5 col mm. Figure 7 shows the optimizations on G , 2 G , 3 G , 4 and G combo (precision-weighted, intensity-weighted, flatness-weighted, and combination, respectively) for the medium focal spot size and = r 5 col mm.Note that the flux units are no longer normalized to the maximum of all the source positions, but rather to the maximum that was found in the corresponding G 1 optimization.This normalization aids in comparing the various objective functions.
Finally, figures 8 and 9 display the Monte Carlo calculations for two selected optimized collimators (one chosen for 'Intensity' and the other chosen for 'Precision') and a parallel-hole collimator (chosen for 'Flatness') for = r 5 col mm and = r 0.5 col mm, respectively.Information on the selected collimators and the parallel-hole collimators for each focal spot size and field size is provided in table 4. For each flux profile shown, the y-axis is normalized to the maximum from each focal spot size and field size (i.e.all the profiles along one row were normalized by the same value).The five curves in each graph correspond to the five source positions: 7, 9, 12, 15, and 20 cm above the detector plane.

Discussion
Validation of equation ( 3) is an important first step in ensuring the robustness of this optimization.The data in figure 3 and the error statistics presented for the small, medium, and large focal spot sizes demonstrate, both Table 4. Monte Carlo-simulated collimator features.This table summarizes which collimators were used in the Monte Carlo-calculated flux profiles of figures 8 and 9 for each focal spot size and field size.In the columns for Intensity and Precision, the listed values are the source positions from which the selected collimators were optimized, i.e., 90 in the first row indicates that the Intensity collimator for the smallest focal spot size and = r 5.0 col mm was optimized when the source was 90 mm from the detector plane.The last column displays the radii of the parallel-hole collimators used for each focal spot size and each field size.graphically and quantitatively, the high level of accuracy our model achieves across three disparate collimation geometries.The Thin Attenuator simulation, where the radiation field is transported through a single 0.1 mm thick lead plate with no opening, isolates the exponential attenuation model from the geometric considerations.Equation (3) uses the total linear attenuation coefficient, thereby assuming narrow beam attenuation.At 200 keV in lead, the photoelectric cross section is ∼6 times greater than that of scattering (Compton and Rayleigh combined) (Attix 1986).This reasoning is complimented by the results presented in figure 3, in which the model calculations adhere tightly to the Monte Carlo results.The 1.0 cm Collimator and 0.5 cm Collimator curves  demonstrate the performance of the model's geometric mechanisms.Both the full phase space and reduced phase space calculations provide results indistinguishable from the Monte Carlo calculations, indicating that equation (3) performs well overall and the reductions made to the phase space were appropriate.
In addition to the accuracy of the radiation transport model, this optimization relies on the robustness of simulated annealing.Figure 4 provides evidence that the optimizer is working as designed.In both figures 4(a) and (b), the objective function oscillates rapidly towards the beginning of the optimization and then steadily climbs towards a stable plateau later on.It is clear from this figure that optima are consistently achieved by the optimizer.Beyond the data presented in these graphs, the purpose of running each optimization three times with three different initial solutions was to test the dependence of the optimizer on starting point.In all of the optimizations, at least two of the three runs outputted the same final solution, providing evidence that the annealing hyperparameters were chosen appropriately, if not conservatively.Runs that did output alternative solutions demonstrate that the objective function can have multiple disjointed local minima, indicating that the optimization is nonconvex.Moreover, direct downhill optimization algorithms will be dependent on the initial solution and thus less robust.Despite the consistency of the optimizer across varying initial solutions, this method of selecting the best output from a set of optimizations is not exhaustive or systematic.When combined with a slow annealing schedule, this initialization strategy appeared to avoid shallower local minima.However, future work will investigate more intelligent means of selecting the optimization starting point.These alternative strategies may include (1) using a previous optimization's output as an initial solution, or (2) developing a convex, differentiable surrogate function that can be optimized by direct methods such as gradient descent, and then using that output as the initial solution for simulated annealing.In addition, alternative optimization algorithms, especially those that are population-based such as genetic algorithms, particle swarm optimization, and ant colony optimization, may offer broader searches of the parameter space (Wu et al 2019).While these population methods typically suffer from worse run times than their single-solution counterparts  (Wu et al 2019), future work will investigate if they offer any advantages in optimization robustness.Currently, the size of the phase space files and the computational cost of calculating the objective function make the singlesolution simulated annealing method a reasonable starting point.
Moving on to the optimization solutions, figures 5 and 6 display the G 1 -optimized collimators for 5 and 0.5 mm radius fields for the smallest, medium, and largest focal spot sizes.The vast majority of these optimizations resulted in conical (pinhole) collimators and some parallel hole collimators, which makes sense given the distributed nature of the radiation source.For = r 5 col mm, the smallest and medium focal spot sizes produced a series of collimators that steadily open to larger radii as the source was moved further from the detector plane.From the perspective that the source size becomes less significant and the radiation field more parallel as the source moves further away, it follows that the collimator opening should approach the desired field size as z 0 increases.As the focal spot size increases and the desired field size decreases, occlusion of the source makes collimation more difficult.Hence, for the largest focal spot size, the optimizer is less consistent in balancing trade-offs between intensity and precision.Thus, in figure 6, the profiles with the greatest flux are not when = z 70 0 mm.Instead, the optimizer had to employ more severe collimation when the source was closer to minimize dose leakage, then, as the source moved further away, the collimators were able to open up and increase the amount of the source that was in view.
The investigation of alternative objective functions in figure 7 elucidates some additional expected conclusions.High intensity and low leakage of flux outside the desired field size is supported by conical collimators.This finding is demonstrated by the optimizations on objective functions G , 1 G , 2 and G .
3 The size and focal length of the collimator allows different weightings towards either intensity or precision, though most of the time G 2 and G 3 produced equivalent solutions to G .
1 In general, these objective functions cause a widening of the penumbra.In traditional 3D conformal radiotherapy, the severe gradients shown in figures 5-7 would be unacceptable for treatment planning and dosimetry.However, the advent of inverse planning and the desire for highly conformal, heterogeneous dose distributions in modern treatment modalities such as intensitymodulated radiotherapy and stereotactic radiosurgery has made non-uniform fields advantageous, especially those of flattening filter free, TomoTherapy, and CyberKnife treatments (Xiao et al 2015).Alternatively, figures 7-9 show that parallel-hole collimators can be used to maximize uniformity within the desired field size.In the G 4 optimizations, flux intensity was preserved while flux leakage worsened.There is some strange behavior in the collimator shape of figure 7 for = z 70 0 mm.This G 4 -optimized collimator takes the form of a parallel-hole collimator in its lower half but has significant jaggedness in the upper half.This may be an overcollimation artifact given the thickness of the lead collimator for a 200 kV source.Despite the lack of intuition in the upper half of these collimators, the resulting flux profile appears reasonable.
The combination optimization shown in figure 7 provides a method of combining the three field shaping priorities-intensity, precision, and flatness.This objective function produced collimators similar to the G 1 optimization with slightly greater flatness within r .col Since the prioritization between intensity/precision and flatness is controlled by the user-defined regularization parameter, b, in equation (8), the user can tune the optimization to get different feasible solutions between G 1 and G .
4 The multi-objective optimization demonstrated here provides a means for developing any number of various collimator goals and tuning their relative importance to the final design.For this particular x-ray field, the optimized collimators were similar across various formulations of the objective function.
With the understanding that conical collimators can be used to prioritize high intensity or low flux leakage while parallel-hole collimators can produce flatter fields, three choices of collimators were selected for Monte Carlo simulation: one that prioritized intensity, one that prioritized precision, and a parallel-hole collimator for flatness, as described in table 4.These results are shown in figures 8 and 9 for = r 5 col mm and = r 0.5 col mm, respectively.For the smallest focal spot size, the choice of collimator was significantly less consequential.The difference in flux profile between intensity and precision profiles is much smaller than what's observed for the other focal spot sizes, owing to the insignificant radiation source size compared to the collimation radius.The profiles for the large and medium focal spot size, however, show relatively disparate flux profiles for each collimator selection, providing a wide range of trade-offs between intensity, precision, and flatness.
One limitation of this work was noted by Wang et al (2021) in the development of their optimization algorithm.Optimizations in this study were performed on photon flux profiles, which ignore the broadening of dose distributions caused by scatter and a long secondary electron range.Essentially, a 1 cm field measured by photon flux may result in a dose profile wider than 1 cm.While this critique is more pertinent to megavoltage radiotherapy where scattered photons and secondary electrons are more penetrating, there should still be concern over the translation from flux to dose.The present algorithm maintains an efficiency and robustness advantage over iterative Monte Carlo-based solutions.Therefore, it may be prudent to use this method for a preliminary optimization where the overall geometry of the collimator is determined (pinhole, parallel hole, etc), and then Monte Carlo dose calculations can be used to perform a local optimization of that geometry's parameters.This would combine the efficiency and robustness of the present algorithm with the clinically relevant dose calculations of Monte Carlo-based optimizations.
As a final note to this discussion, optimization was performed for 10 and 1 mm field sizes for exemplary purposes.Larger and intermediate field sizes can be readily optimized with the same phase space files and optimization hyperparameters.However, = r 0.5 mm col approaches the lower limit of what the phase space's resolution can handle.For smaller field sizes, e.g.0.1 mm, ¢ ¢ x y d , d ,and x d det in table 2 and r d in table 3 would need to be decreased by at least one order of magnitude.The spatial ranges of the phase space could be decreased accordingly to reduce the size of the phase space file.However, the Monte Carlo calculation of ( ) ¢ ¢ f x y x E , , , det would suffer from worse statistical noise and thus need to be run for many more histories.Future work could involve increasing the efficiency of this optimization algorithm to handle larger phase space files or reduce computation time, thus allowing us to extend the range of field sizes for which using this algorithm is a reasonable approach or even relax some of the stated symmetry requirements.

Conclusion
An optimization algorithm was developed for designing cylindrically symmetrical collimators for the novel x-ray tube described by Insley et al (2024).Advantages of this algorithm include its analytical radiation transport model, which precludes the need for iterative Monte Carlo approaches.Computational efficiency is greatly increased, resulting in the ability to reduce constraints on the collimator shape, increase the number of free parameters, and increase the range of exploration in the collimator parameter space.This study was able to demonstrate the accuracy of our algorithm's radiation transport model compared to that of Monte Carlo, run optimizations for a wide range of focal spot sizes and source positions, investigate multiple formulations of the objective function, and analyze the performance of the optimized collimators with Monte Carlo simulations.As radiation sources increase in complexity and specificity, this algorithm provides a computationally efficient and robust method of designing custom collimating apparatuses without making unnecessary assumptions about the radiation source or the optimal collimating shape.This work provides a clear starting point for future efforts that seek to increase the source, collimator, objective, or optimization complexity.

Figure 1 .
Figure 1.Radiation transport model diagram.The blue cone at the top of the figure represents the novel x-ray tube, courtesy of Insley et al (2024).Electrons are emitted from the upstream CNT cathode and accelerated by a 2-stage voltage configuration.The green cone displayed at the bottom of the tube is the conical x-ray target receiving the tube's 200 keV electron beam.The parameters outlined in the diagram below the target are those described in equations (1) and (2).Photon tracks are defined by starting position, ( ) ¢ ¢ x y , , ending position, ( ) x , 0 , det and energy, E. Collimator slices are defined by their position above the detector plane, z ,i and inner radius, r .iThe system is set with a constant distance between the phase space plane and detector plane, z 0 .
det was then calculated as a four-dimensional histogram using the Monte Carlo code Geant4 via its wrapper, TOPAS (v.3.8.1)(Perl et al 2012, Faddegon et al 2020).Photons are generated from direct simulation of the x-ray tube electron beam (please see Insley et al 2024) for a full description of the simulation geometry).The default physics list, including Geant4's G4EMStandardPhysics_option4 module for electromagnetic interactions, was used with a range cut of 0.1 μm for all particles.The variance reduction parameters established by Ali et al for orthovoltage x-ray tube simulations were employed (Ali and Rogers 2007).This included bremsstrahlung particle splitting (

( 4 )
Recalculate the objective function for each change in washer radius.If the objective function value has improved, accept the change.Otherwise, accept the change according to the probability ( ) is the change in objective function values before and after the radius change, ( ) When this count has reached some maximum value, D , max terminate the optimization.The hyperparameters outlined above were determined based on the methodology outlined byPark and Kim (Park and Kim 1998) and summarized in table3.With the optimization process fully established, this section will conclude with a summary of all the x-ray tube and field size settings plugged into the algorithm.The x-ray tube was characterized in the literature for nine focal spot sizes(Insley et al 2024).Optimization was run using objective function G 1 for each focal spot size for two different field sizes: time, optimization was run thrice for three different initial solutions, all parallel hole collimators:

Figure 2 .
Figure 2. Phase space intensity map of the large focal spot size.The color of each pixel represents the phase space model, ( ) ¢ ¢ f x y x E , , , , det integrated over energy and x .detThe red bounding box shows the section of the phase space that is kept after the spatial reduction.

Figure 3 .
Figure 3. Model validation.The performance of equation (3) is demonstrated in three separate scenarios by the three sets of curves displayed above: a single lead sheet with no central hole ('Thin Attenuator'), a 1.0 cm diameter, 1.0 cm thick parallel-hole lead collimator ('1.0 cm Collimator'), and a 0.5 cm diameter, 1.0 cm thick parallel-hole lead collimator ('0.5 cm Collimator').For each curve, the Monte Carlo-calculated flux is shown in blue crosses, the model calculations are shown by the red-orange line, and the model calculations using the reduced phase space are shown by the yellow line.The yellow line is plotted on top, so places where the red-orange line is not visible is because of very strong agreement between 'Model' and 'Model (Reduced)'.

Figure 5
Figure 5 displays the optimized collimators and resultant photon flux at the detector plane for the smallest, medium, and largest focal spot sizes, respectively, for = r 5 col mm.Each figure contains three of the five source positions described previously: = z 70, 120, 0 and 200 mm.The top row of charts in each figure contains crosssections of the optimized collimators, colored according to source position, in which the y-axis describes each washer's position above the detector plane, z ,i which ranges from 60 to 70 mm.The x-axis gives the radius measured from the x-ray tube's longitudinal axis, r .iThe bottom row contains the resultant photon flux at the detector plane with normalized flux units on the y-axis and distance along the x det axis on the x-axis.The y-axis is normalized by the maximum flux of all the source positions for a given focal spot size.Figure6presents the same data for = r 0.5

Figure 4 .
Figure 4. Optimization progress examples.Each graph displays the progression of the objective value or a single factor in the objective function as a function of iteration.In (a) and (b), the G 1 objective function is maximized during optimization, thus, both progressions trend upward.The run depicted in (a) was for the smallest focal spot size, = r 5 col (flatness) throughout the optimization, respectively.

Figure 5 .
Figure 5. G 1 -optimized collimators and photon flux for = r 5 col mm.Top row shows longitudinal cross-sections of the optimized collimators.Bottom row shows the resultant photon flux at the detector plane.Source positions of each optimization are denoted by color: blue for 70 mm, red for 120 mm, and green for 200 mm.In the medium focal spot size collimator plot, the 70 mm collimator is only slightly narrower than the 200 mm collimator, and is therefore difficult to see in the plot.

Figure 7 .
Figure 7. Objective function comparison.This figure is organized like figures 5 and 6 with the top row showing collimator cross sections, the bottom row showing flux profiles, and color indicating the source position.However, columns now organize the different objective functions: primary (G 1 ), Intensity-Weighted (G 2 ), Precision-Weighted (G 3 ), Flatness-Weighted (G 4 ), and Combination (G combo ).These optimizations were performed for the medium focal spot size with = r 5 col

Figure 8 .
Figure 8. Monte Carlo profile calculations for = r 5 col mm.The profiles for the 'Flatness' column are created by the parallel hole collimator.Each row of graphs are normalized by the maximum from that row.

Figure 9 .
Figure 9. Monte Carlo profile calculations for = r 0.5 col

Table 2 .
Flux model parameters.Listed are the various parameters' resolutions and ranges for the discretization of the photon flux model in equation (3).The z i used in the discrete version of equation (2) is defined as the bottom of each washer.z d was chosen to be thin enough so that this choice of z i definition is insignificant to the calculation of equation (3).The square brackets on the ranges indicate inclusivity, while rounded brackets indicate exclusivity.