text:
# Algorithm and performance of a clinical IMRT beam-angle optimization system

## 1 Introduction

The beam-angle optimization (BAO) in intensity-modulated radiation therapy (IMRT) aims to find, for each patient-specific case, a set of beam angles that can produce a three-dimensional (3D) dose distribution that conforms most closely to the prescribed one. In general, selecting the best set of beam angles that can achieve this goal is a highly non-intuitive process. This is already true in 3D conformal radiation therapy (CRT), where each beam has a preset profile inside its beam’s-eye-view aperture and only the beam weights are optimized, but especially more so in IMRT. The optimal distribution of the total dose in IMRT is the result of a complex interplay among a large number of beamlet intensities from multiple beam angles and therefore there is no a priori way of predicting the dependence of the objective function (OF) on the beam angles. The optimizer has to balance the requirement of creating a prescribed dose distribution at the target with that of giving as little dose as possible to the surrounding normal tissue and organs at risk (OARs). Stein et al(1997) provided a well-known counter-intuitive example of the necessary compromise that often has to be struck to achieve this balance between target coverage and OAR sparing: sometimes the beams have to go directly through an OAR to achieve an optimal distribution.

The complex dependence of the dose distribution on the beam angles and their beamlet intensities gives rise to a multiple local minima problem. The fact that the OF has a highly non-convex dependence on its variables has so far precluded any deterministic approach to its minimization and necessitated a brute-force assail to this problem, namely a global search method. This approach requires one to perform an IMRT optimization on each one of a large number of possible combinations of beam angles. Two main hurdles have to be overcome, however, in making this highly-iterative approach practical and can be used routinely in a clinical setting. The first hurdle is the substantial amount of time that is required to perform an IMRT beam-profile optimization for a single set of beam angles. Depending on the complexity of the case and the algorithm used in the optimization, the required time already ranges from several minutes to a few hours on a typical workstation. The second hurdle is the huge size of search space involved in the beam-angle search. To rigorously search for, e.g., the best 5 gantry angles out of 180 possible directions would require one to search over $`C_{180}^5(10^9)`$ possible combinations. A BAO system that is based on an exhaustive search would therefore, most likely, be too slow for utilization in a clinical setting.

The recent past has seen a limited but steady flow of publications on this topic. A few papers have provided theoretical analysis of the BAO problem (Bortfeld and Schlegel 1993, Jackson et al1994, Braunstein and Levine 2000). Bortfeld and Schlegel (1993) considered the BAO for the 3D CRT, as opposed to IMRT, and concluded that, for dose-based criteria, the optimum beam configuration for multiple-beam irradiation tends to be an even distribution over the $`2\pi `$ angular range for the gantry.

Rowbottom et al(1998, 1999a, 1999b) have proposed several approaches for BAO as applied to the 3D CRT. Dose-based OF was used in the optimization of the beam weights and simulated annealing algorithm was used to select the beam orientations in a reduced search space of coplanar (Rowbottom et al1998) and non-coplanar (Rowbottom et al1999a) beams. Some a priori criteria were also used in the beam-orientation selection: the beam orientations were constrained to be separated by at least $`35^{}`$ and parallel-opposed beams were disallowed. In addition, a single-beam cost function (Oldham 1998), based largely on the geometry of the target, was used to rank the beam orientations in the search space according to their likelihood for being a member of the optimal set and only the top 10% of the set were searched in the BAO. Rowbottom et al(1999b) also attempted a more intelligent search using an artificial neural network technique. The results of this search, however, seem to be inferior to those that were obtained using their earlier simulated annealing method.

Pugachev and Xing recently proposed several approaches to BAO for the IMRT (Pugachev et al2000, 2001, Pugachev and Xing 2001, 2002). In their earliest proposal, a non-iterative filtered back projection (FBP) method, using simple exponential attenuation for parallel X-ray beams in their dose calculation, was used to solve the inverse problem in IMRT while the beam orientations were optimized using a simulated annealing method. This fast FBP approach, however, seems to preclude the use of OFs other than a dose-based OF. In Pugachev et al(2001), an iterative method (Xing and Chen 1996) was used to optimize the beam profiles and non-coplanar beams were used in the search space. A typical prostate case in this study required about 100 hours to be optimized by searching over 5000 combinations of beam angles. The use of the beam’s-eye-view dosimetrics (BEVD) (Pugachev and Xing 2001, 2002) allowed them a more efficient sampling of the search space: The sets of beams that have low BEVD score, i.e. those that are not likely to produce an optimal dose distribution, are immediately rejected without going through the beam-profile optimization. It was claimed that the BEVD-guided sampling can improve the speed by a factor of 10 while maintaining the same capacity to locate the global minimum as that of an unguided sampling.

It is interesting to note that most of previous approaches to BAO used the dose-based OF (Bortfeld and Schlegel 1993, Rowbottom et al1998, 1999a, 1999b, Pugachev 2000, 2001, Pugachev and Xing 2002). Das et al(2003) proposed an IMRT beam orientation selection method that is based on target equivalent uniform dose (EUD) maximization. The objective is to maximize the EUD on the target while simultaneously imposing DV constraints on the OARs. These two objectives were combined into a single OF by a simple scale factor and it would be interesting to compare the performance of this combined OF with that of a strictly DV-based or EUD-based OF.

The approach that we take in this paper is to use a fast dose calculation engine that significantly reduces the time required for the IMRT optimization of each given set of beam angles. This algorithm relies on the use of a compact, albeit approximate, dose kernel (Wu et al2003). The main purpose of this paper is to show that this fast algorithm is sufficiently accurate to be used as a dose engine in a BAO. Although simulated annealing has been used in several previous works on BAO, in this paper we show that the combination of simulated annealing, in particular the FSA algorithm, and the fast dose engine that we proposed recently can provide a practical, i.e. fast and useful, BAO system. In addition, by focusing the improvement on the speed of the dose engine in our system, we are also able to retain all other aspects of normal, i.e. clinical, IMRT optimization in our BAO. Thus our BAO system allows a DV-based OF, or indeed any other types of OF, to be used for the IMRT optimization for each given set of beam angles generated during the BAO. The promise of a BAO system is to provide better alternatives to plans produced by a planner, which in clinical settings are commonly generated by a DV-based optimization. It is therefore important to use a DV-based, instead of a dose-based, optimization in a clinical IMRT BAO system so that a direct comparison can be effected. If the BAO plan is deemed by the planner to be significantly better, it can be used for treatment; otherwise, the plan produced manually by the planner is retained. Thus the BAO can be used to automatically check if the best plan of the planner can still be improved by beam-angle adjustment. Starting with an already acceptable plan, the BAO is executed with the same DV constraints and parameters and its plan will be used if it is significantly better and can be produced within an allowable time.

It should be noted that no a priori assumption is made regarding the selection of beam orientations in this work. Although incorporation of some criteria, such as constraints on the angular separation between adjacent beams, is likely to improve the performance of the system, we feel that such topics require separate careful studies. These additional beam-selection criteria will be incorporated in further studies.

## 2 Methods and materials

In this work, we propose an efficient BAO method that combines a fast dose calculation engine, which we call the table lookup (TLP) method (Wu et al2003), with an effective global search method based on the fast simulated annealing (FSA) algorithm. We have recently implemented this method in a program that interfaces with our in-house clinical IMRT system and studied its performance and potential.

The following subsections will detail the components of our BAO implementation. Section 2.1 discusses the flow diagram of our BAO system and the algorithm that we use to select successive sets of beams from the search pool (the FSA algorithm). The dose engine that we use as a workhorse in the BAO (the TLP method) is presented in section 2.2 while section 2.3 describes the DV OF that we use to guide the IMRT optimization. Details on how the TLP method is used to accelerate our IMRT optimization and discussion on its accuracy are provided in section 2.4. The cases that were used to test our BAO system in this work are outlined in section 2.5.

### 2.1 Beam angle selection with fast simulated annealing

Figure 1 displays the flow diagram of our BAO system which is interfaced with an in-house IMRT system. The IMRT system, in which the planner sets general optimization parameters including the number and energy of beams whose directions are to be optimized, is in turn interfaced with a commercial treatment planning system (TPS) Pinnacle3 (Philips Oncology System, CA, USA). This in-house IMRT BAO program, which is external to the TPS, communicates with the TPS through an application programming interface, provided by the TPS vendor, called PinnComm. This allows the BAO and the IMRT programs to, e.g., request the TPS to perform a dose calculation.

The BAO is performed using an FSA algorithm. Prior to entering the SA loop, the planner defines a search pool of beams. We use 180 equispaced beams for a coplanar search pool and 216 beams (72 equispaced beams at couch angle of 0, 45, and 315) for a non-coplanar search in this work. For each beam in the search pool, the TPS is called by the BAO program to calculate an accurate dose for a uniform intensity profile. The jaws are automatically set for each beam to sufficiently cover the target with some margins. The kernel is then extracted from the TPS using the TLP method (section 2.2). Note that we do not optimize the number of beams in this study.

#### 2.1.1 Fast simulated annealing: temperature lowering schedule.

The BAO selects the optimal set of beams from the search pool. Successive sets of beam angles are generated during the BAO by using the FSA algorithm (Szu and Hartley 1987). The efficacy of this global search method has been proven in many fields, including IMRT and beam-angle optimization (Webb 1989, 1991, Mageras and Mohan 1993, Xiang and Gong 2000, Pugachev 2000). The FSA is known to be much more efficient than the classical SA (CSA) algorithm (Kirkpatrick et al1983). To converge to a global minimum with probability one, the CSA requires the annealing temperature to be lowered with an inverse-logarithmic schedule: $`T=T_0/\mathrm{ln}(1+i)`$, where $`i`$ is the iteration number (Geman and Geman 1984). On the other hand, the FSA allows the temperature to be lowered with a much faster schedule (Szu and Hartley 1987):

$$T=\frac{T_0}{1+i}.$$
(1)
The SA temperature $`T`$ is used in determining whether a new set of beam angles is to be accepted or rejected. Note that temperature in this context should be understood only as an optimization parameter and it does not have any relation, except in analogy, to the temperature of any physical system. A move generated by the FSA is accepted with a Metropolis probability (Metropolis et al1953):

$$p=\mathrm{min}[1,e^{\mathrm{\Delta }E/T}],$$
(2)
where $`\mathrm{\Delta }E`$ is the change of the OF (score). Thus a downhill move is always accepted while an uphill move is accepted with a probability $`e^{\mathrm{\Delta }E/T}`$. The Metropolis probability allows the optimization to simulate the behavior of a physical system in a thermal equilibrium while the SA temperature schedule controls the lowering of temperature such that the system is effectively at a thermal equilibrium at any time. As in a physical system, this recipe ensures that local minima are properly sampled and the global minimum is reached as the temperature is lowered to zero.

It should be noted that $`\mathrm{\Delta }E`$ in equation 2 is the difference between the optimal score (final value of objective function) of the current trial set of beams and that of the last accepted set of beams. Each of these scores is the output of a beam-profile optimization performed on the corresponding set of beams. Details of the beam-profile optimization is not needed by the FSA algorithm whose task is to generate subsequent sets of beams and accept or reject each set based on the value of its optimal score.

#### 2.1.2 Fast simulated annealing: sampling distribution.

The relative efficiency of the FSA over the CSA is related to the sampling distribution that it uses. At each BAO iteration, FSA samples the changes to the current beam angles $`𝐱=(\theta _1,\theta _2,\mathrm{},\theta _n)`$ from a long-ranged Cauchy-Lorentz distribution:

$$p(\mathrm{\Delta }𝐱)[(\mathrm{\Delta }𝐱)^2+W(T)^2]^{(n+1)/2},$$
(3)
where $`n`$ is the number of beam angles (the dimension of the vector x). In contrast to FSA, the CSA uses a short-ranged Gaussian sampling distribution: $`p(\mathrm{\Delta }𝐱)\mathrm{exp}[(\mathrm{\Delta }𝐱)^2/2W(T)^2]`$, where $`W(T)`$ is the temperature-dependent width of the sampling distribution. Search-space regions farther away than $`W(T)`$ from the current position are practically not sampled by the CSA distribution and therefore CSA requires a larger number of iterations than FSA to thoroughly sample the entire search space.

Several different choices for the definition of $`W(T)`$ are available in the literature (Mageras and Mohan 1993) and, in general, the “best” choice for an application can only be found by trial and error. For the FSA work reported in this paper, we have used a simple form $`W(T)^2=\alpha T^2`$. The value of the range parameter $`\alpha `$ was chosen such that even at the lowest temperature ($`T=1`$ in this paper) each beam direction still has 1% probability to jump to another direction “farthest away” from it (e.g., from 0 to 180 gantry angle). Thus at any temperature, the current trial set of beams always has a finite probability to jump to any other point in the search space. The ability to sample regions beyond the local neighborhood is imperative in the simulated annealing method to avoid being trapped in a local minimum.

In generating the random move for the set of beam angles $`𝐱=(\theta _1,\theta _2,\mathrm{},\theta _n)`$ we use the approach where the orientations of all beams in the current set are simultaneously updated. This has been shown (Mageras and Mohan 1993) to be more effective than the alternative approach, in which one beam is picked randomly from the current set and its orientation is moved by a random distance which is picked from the sampling distribution (Das et al2003).

### 2.2 Fast dose-calculation method

The total dose at a given voxel, $`D_i`$, is related linearly to the beamlet weights $`W_j`$ by the dose kernel $`K_{ij}`$:

$$D_i=\underset{j}{}K_{ij}W_j.$$
(4)
Note that in this equation the beamlet index $`j`$ comprises also the beam index; the index $`j`$ runs over all beamlets in all beams. The dose kernel $`K_{ij}`$ is a sparse matrix. At each dose voxel, significant contributions to its dose come only from a limited number of beamlets in each beam, the largest contribution being from the beamlet that directly traverses the dose voxel. The sparse nature of the kernel matrix can be shown numerically by, e.g., plotting the magnitude of the elements of the square matrix $`K^tK`$ (Crooks and Xing 2001). This also gives a measure of the overlap between beamlets.

The fact that the kernel matrix is sparse has been used to devise a fast, albeit approximate, table-lookup (TLP) dose calculation method, that is useful for IMRT optimization (Wu et al2003). This efficient method was proposed to improve the performance of applications that involve a great number of dose calculations. In the present work we have used this method as the fast dose calculation engine for our IMRT BAO. In the TLP method, each dose voxel is assigned one, and only one, beamlet from each beam, namely the beamlet whose raytrace intersects the voxel. The dose to the voxel is initially calculated using an accurate method, e.g., the superposition/convolution (SC) method, for a uniform intensity profile for all beams. The TLP kernel is then calculated using the following formula:

$$K_{ij}=\{\begin{array}{cc}D_i/W_j,\hfill & \mathrm{if}\mathrm{voxel}i\mathrm{is}\mathrm{on}\mathrm{the}\mathrm{path}\mathrm{of}\mathrm{beamlet}j,\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}$$
(5)
Once the kernel is calculated, dose calculation in subsequent iterations is done by looking up the kernel. Note that in addition to improving the speed of the calculation, the TLP method also greatly reduces the number of kernel elements that need to be saved and read, and therefore the time needed to access them. This number, for each beam, is equal to the number of dose voxels because of the one-to-one dose-beamlet correspondence (Wu et al2003). The use of this fast dose-calculation method in our IMRT optimization will be further described in section 2.4.

### 2.3 Objective function

For each trial set of beam angles generated during the BAO, an IMRT optimization is performed to optimize the beam profile. This IMRT optimization is implemented by using an objective function, or score, that provides a measure of “distance” between the current dose distribution and the desired distribution. Starting from an initial distribution, the objective function is iteratively minimized by the optimizer.

Several types of objective function (OF) exist in the literature and implemented in clinical IMRT optimization systems. Among these are the dose-based, DV-based, and the biology-based functions. The dose-based OF is the most straightforward conceptually and the easiest to implement (Webb 1989): it uses a penalty which is a quadratic function of the difference between the actual dose and the desired dose, which is usually the prescription dose for the tumor and zero for the organs at risk (OARs).

The DV-based objective function is among the ones that are most useful in IMRT optimization (Bortfeld et al1997, Spirou and Chui 1998, Wu 2000). It has the advantage of being more intuitive than the biology-based OFs and more flexible than, since it is the superset of, the dose-based OF. The DV-based optimization typically requires more iterations to achieve convergence than the dose-based optimization. Furthermore, it requires the dose-volume histogram (DVH) to be computed after each dose calculation which adds to the total optimization time. In general, however, DV objective requires less adjustment of the optimization parameters to achieve a suitable plan (Wu 2000). Although multiple local minima exist in a DV-based IMRT optimization (they are absent in dose-based optimization), it has been shown (Wu and Mohan 2002, Llacer et al2003) that they are not expected to affect the outcome of optimization using gradient techniques in any clinically significant way.

For the purpose of establishing our notation for later reference in the rest of this paper, here we present a brief description of our implementation of the DV OF. There are two types of constraint in a DV-based optimization, each is specified by a set of dose-volume values $`(D_a,V_b)`$. The first type of constraint (upper constraint) requires the volume receiving dose greater than $`D_a`$, which is denoted by $`V(D_a)=V_a`$, to be less than $`V_b`$. This constraint is implemented by searching for the dose value $`D_b`$ such that $`V(D_b)=V_b`$. A penalty is charged if $`V_a>V_b`$, which is equivalent to $`D_a<D_b`$ since the DVH is a non-increasing function of dose. The upper constraint for a DV-based objective function therefore has the following form:

$$f_1=p_1\underset{i}{}H(D_bD_i)H(D_iD_a)(D_iD_a)^2,\mathrm{upper}\mathrm{constraint}.$$
(6)
Here $`H(x)`$ is a Heaviside step function defined by

$$H(x)=\{\begin{array}{cc}1,\hfill & x>0,\hfill \\ 0,\hfill & x0.\hfill \end{array}$$
(7)
An upper constraint is used to limit the dose to an OAR or the hot-spot regions in a target. The index $`i`$ in the summation runs over all voxels in the OAR, while $`p_1`$ is the penalty prefactor. The constraint’s objective is to reduce the volume of the OAR that receives dose greater than $`D_a`$; it therefore penalizes voxels that have dose $`D_i>D_a`$. Note, however, that voxels that have dose greater than $`D_b`$ are not penalized since the total volume of such voxels is already lower than $`V_b`$.

A lower constraint is used to limit the cold-spot regions in a target. This constraint requires $`V_a`$, the volume receiving dose greater than $`D_a`$, to be greater than $`V_b`$. The penalty is imposed if $`V_a<V_b`$, or equivalently if $`D_a>D_b`$:

$$f_2=p_2\underset{i}{}H(D_aD_i)H(D_iD_b)(D_iD_a)^2,\mathrm{lower}\mathrm{constraint}.$$
(8)
In an IMRT optimization, DV constraints are assigned by the planner to the volumes of interest (VOIs). Each VOI can have more than one constraint; these constraints control the shape of the DVH curve for the corresponding VOI. Since there is no specific desired dose for each VOI in a DV-based optimization, in contrast to dose-based optimization, the quality of the DVH of one or more VOIs can be sacrificed slightly by the optimizer to improve those of some other VOIs. The importance of each DV constraint, relative to the other existing constraints, is set by the relative value of its penalty prefactor.

### 2.4 Intensity optimization for each set of beam angles

For each trial set of beam angles in the BAO, an iterative IMRT optimization using the Newton’s method is performed to optimize the beam profiles by minimizing a DV-based objective function. At each iteration during the IMRT optimization each beamlet weight $`W_j`$ is updated to $`W_j+\delta W_j`$. The size of this update $`\delta W_j`$ depends on the current values of dose and the kernel: $`\delta W_j=F[\{D_i\},\{K_{ij}\}].`$ The procedure for computing $`\delta W_j`$ can be found in Wu and Mohan (2000) and will not be repeated here. Our IMRT optimization program is interfaced with a commercial 3D TPS. For a given set of beams, we start from a uniform intensity profile for all beams and compute an accurate dose $`\{D_i\}`$ using the dose calculation engine (which is based on the SC method) of the TPS. Kernel $`\{K_{ij}\}`$ is then calculated for the iterative optimization using the TLP method (equation 5).

#### 2.4.1 Intensity Optimization with Superposition Convolution (IO-SC) Method.

Note that dose and kernel enter independently in the computation of the beamlet weight update $`\delta W_j`$. This allows different levels of approximation to be used for dose and the kernel. Ideally both dose and kernel are computed accurately using, e.g., an SC method. However, for an IMRT system that is interfaced with a commercial TPS, the kernel is usually not directly available from the TPS, although the dose is. In the IO-SC method (Wu and Mohan 2000), $`\delta W_j`$ is obtained using an accurate dose and an approximate kernel at each iteration. The accurate dose, that incorporates scatter effects, is obtained from the SC dose engine of the TPS. The approximate kernel, on the other hand, is obtained from the accurate SC dose using equation 5. Although scatter effects are neglected in this approximate kernel, they are partially accounted for in the calculation of $`\delta W_j`$ since a full-scatter dose is used. The TLP method of constructing an approximate kernel has been found to be useful in practice (Wu et al2003) and this IO-SC method has been used for IMRT routinely in our clinic. The main limitation of the IO-SC method is the fact that the dose has to be computed at each iteration using the SC dose engine from the TPS. The total optimization time is then strongly dominated by the lengthy time spent on repeated dose calculations.

#### 2.4.2 Intensity Optimization with Table Lookup (IO-TLP) Method.

For the purposes of BAO, we introduce a further approximation to the computation of $`\delta W_j`$ which provides a significant speed up in the calculation time. Two different methods of dose calculation are used in the IO-TLP method. An initial dose is calculated accurately using the SC method of the TPS. A kernel is then generated from this accurate dose using equation 5. In subsequent iterations, instead of computing a fully accurate dose using the SC method, dose is calculated by looking up the kernel, i.e., using equation 4. Thus both dose and kernel are approximate in the IO-TLP method as opposed to only the kernel being approximate in the IO-SC method. Note that scatter effects are not totally absent in the IO-TLP method because a full-scatter dose is always used as the initial dose for the beam-profile optimization.

The TLP method provides a substantial speed advantage over the SC method. Since it is based mostly on memory lookup operation, each TLP dose calculation, which is the dose engine in the IO-TLP method, is typically at least two orders of magnitude faster than the corresponding SC calculation, which is used in the IO-SC method. In comparing the speed of these two optimization methods, we need to consider two separate situations: (1) a single IMRT optimization for a given set of beams; and (2) the BAO in which a large number of sets of beams (1000 in this work) are optimized consecutively. In the single IMRT case, the IO-TLP method uses the TLP dose engine to calculate the total dose at each iteration except the first one where SC is used to compute the kernel. The IO-SC method, on the other hand, uses the SC dose engine at each iteration. The speed ratio between the IO-TLP and IO-SC methods in this case is given roughly by the total number of iterations needed to converge the optimization, which is usually about 10–20 in our IMRT system. In the case of BAO, the kernel for each beam in the search pool is computed, using the SC method, prior to any IMRT optimization. During the beam-profile optimization of each given set of beams, no SC dose calculation is needed because the kernels have all been computed and can be quickly read from memory by the TLP method. Except for the time overhead spent in the initial calculation of the kernels, which is small relative to the total BAO time (table 2) the speed ratio in this case is given directly by the speed ratio between the TLP and the SC methods. Thus our BAO system is at least two orders of magnitude faster than it would be if we used the SC method as the dose engine for beam-profile optimization.

The feasibility of using this approximation for a BAO is demonstrated in Fig. 2 where we compare the optimal scores obtained using the IO-SC method with the corresponding scores obtained using the IO-TLP method. This plot has been generated using 100 randomly selected sets of 3 beams that were used to plan a typical prostate case. Each score shown in Fig. 2 is the square root of the optimal value (according to the respective optimization method) of the objective function which is the sum of contributions from all DV constraints applied to the VOIs in this case (Wu and Mohan 2000). It can be seen that the IO-TLP scores are well correlated with the IO-SC scores: “good” or “bad” sets of beam angles according to the IO-SC method are predicted similarly by the IO-TLP method. In general the IO-TLP scores are lower than the IO-SC scores. This is due to the fact that less scatter components are included in the IO-TLP method than in the IO-SC method. Dose “leakage” outside the target, due to scatter, is slightly underestimated by the IO-TLP method and this gives rise to the undercounting of the score and the finite variance of the correlation. It should be noted that the systematic score undercounting in the IO-TLP method is practically irrelevant to a BAO since only the relative ranking of different sets of beams is needed. The presence of a finite variance in the correlation indeed introduces some limitation to the predicting power of the IO-TLP method. It is the purpose of this paper, however, to show that even with this limitation the IO-TLP method can still be used as an effective, and efficient, IMRT workhorse in a BAO system.

It can be seen in figure 1 that our BAO algorithm simply serves as an additional “outer loop” to a normal IMRT optimization (Wu and Mohan 2000). The beam-profile optimization component of the BAO, however, uses a different dose-calculation engine than the one that is used in normal IMRT optimization. The IO-TLP method is used as the beam-profile optimizer during a BAO, because of its speed, while a normal IMRT optimization (for a given set of fixed beam angles) should be performed using the IO-SC method, because of its accuracy.

It should be noted that although the IO-TLP method is used during the BAO, a final beam-profile optimization using the IO-SC method is performed on the optimal set of beam angles output by the BAO. This is a part of the finalization step shown in the last box in figure 1. This final optimization of the BAO plan facilitates its comparison with a human-prepared plan, which is also optimized using the IO-SC method. In the following section, all of the displayed results have been produced using a IO-SC-based optimization.

### 2.5 Materials

Two cases have been used in this work to test the practicality of our BAO system: a prostate case and a head-and-neck case. Some of the system parameters of these two cases are displayed in table 1. The prostate case represents the cases with relatively small size and simple arrangement of the target and the OARs. The head-and-neck (HN) case represents the other end of the spectrum namely those with relatively large size and with multiple targets that are located close to sensitive OARs. For the prostate case we compare the BAO results for 3, 5, and 7 coplanar beams with that of the corresponding number of equispaced coplanar beams. We also examine the dependence of the improvements that can be obtained from BAO on the criteria used to guide the optimization. This is done by studying the results of two different sets of criteria: a “soft” and a “hard” set. For the HN case we compare the BAO results for 7 and 9 non-coplanar beams with that of 9 equispaced coplanar beams.

The computational time required to perform a BAO on a given case generally scales with the system size. This is determined principally by the product of the total number of beamlet intensities that need to be optimized and the total number of voxels for which dose needs to be computed during the optimization. Note that table 1 lists the typical size of the intensity matrix (which is the size of the smallest rectangle that encloses the target in the BEV of the beam) since the specific size generally varies with the orientation of the beam. Furthermore, the number of beamlet intensities that are optimized is generally less than the size of the intensity matrix because the beamlets that do not pass through any voxel in the target can be completely turned off. Although increasing the number of optimization criteria is expected to increase the total optimization time, the increase is minimal because once the dose is calculated the value of the OFs can be computed very quickly. The size of the bitmap vector for the HN case in table 1 is about 4 times that of the prostate case and the size of its intensity matrix is about $`4\times 3`$ times larger. Thus, for the same number of beams, the total IMRT optimization time for the HN case is expected to be about 48 times that of the prostate case.

## 3 Results

### 3.1 Prostate case

Figures 4 and 5 display the improvements that can be obtained from BAO for a prostate case. Each of the results shown in these figures has been obtained from a 1000-iteration FSA run and was completed in about 1–2 hours on a 400 MHz workstation (table 2). The initial SA temperature $`T_0`$ was set to an arbitrary large number ($`T_0=1000`$) which allowed practically all random moves generated by the SA sampling distribution to be accepted. This is an important practical requirement of the SA algorithm: The initial temperature has to be high enough to allow random sampling of the search space during the initial iterations. Subsequently, the temperature was lowered using a $`T_0/(1+i)`$ schedule. Figure 3 displays a typical progression of the score with increasing iteration number. The optimization is stopped after 1000 iterations (final temperature is 1) because practically the temperature is already too low to allow a further significant change in the score to occur. Note that, based on our experience as shown in figure 3, the optimization can generally be terminated after a few hundred iterations without incurring a significant cost in the final score.

The optimal beams for this prostate case have been chosen from 180 coplanar 18-MV photon beams (at zero couch angle) which are distributed uniformly over the $`360^{}`$ range for the gantry angle with a $`2^{}`$ interval. To facilitate the comparison, each plan shown has been normalized to 72 Gy at 90% of the target volume. The DV-criteria that were used for the IMRT optimization are listed in table 3.

Figure 4a compares the DVHs of a set of 3 equispaced beams and those of a set of non-equispaced beams selected by the BAO. It can be seen that the BAO beams provide significant improvements over the equispaced beams. The target receives a more uniform coverage while both of the OARs, the bladder and rectum, receive significantly less dose. Since in a 3-beam equispaced set each beam is separated from its neighbors by a relatively large $`120^{}`$ interval, it is important to examine if better results can be obtained from equispaced beams by uniformly rotating the complete set. This is done in figure 4b where we compare 4 different sets of equispaced beams, each equispaced configuration is rotated by $`30^{}`$ from the previous one. While their target coverages are similar, they spare the OARs rather differently. It is interesting to observe that the set that provides better protection for the bladder does so at the expense of the sparing of the rectum, or vice versa. In addition, all of these plans are still inferior to that of the BAO beams shown in figure 4a. These comparisons illustrate the importance of non-equispaced beam-angle selection in IMRT.

DVH comparisons between equispaced and BAO plans for the same prostate case using 5 and 7 beams are presented in figures 5a and 5b respectively. These plans have been optimized using the same set of DV criteria as that used for the 3-beam plans in figure 4. Although it is expected that the importance of BAO will diminish as more beams are used (Stein et al1997), these cases show that BAO can still provide rather significant improvements over the equispaced beams. We have also performed similar comparisons as in figure 4b between different sets of 5 and 7 equispaced beams. As is expected, the difference among different 5- and 7-beam sets is smaller than that for the 3-beam set and it decreases as we increase the number of beams from 5 to 7.

### 3.2 Head-and-neck case

Due to its anatomical position, head-and-neck tumors are good candidates for the application of non-coplanar IMRT. Limited experience with non-coplanar HN plans and the fact that selection of non-coplanar beam directions for IMRT is a highly non-intuitive process also make global search methods like the FSA, which is implemented in our BAO system, particularly attractive.

Figures 6 and 8 compare the DVHs for a HN case of 9 and 7 non-coplanar beams, respectively, with those of a benchmark 9 equispaced coplanar beams. The DV criteria used to optimize these HN plans are listed in table 5 while the gantry and couch angles selected by our BAO algorithm for these sets are tabulated in table 4. Each of the BAO results for this HN case has been obtained from a 1000-iteration FSA runs. The optimal beams have been chosen from 216 possible directions: the search pool for this case consists of 3 possible couch angles ($`0^{}`$, $`45^{}`$, and $`315^{}`$) with 72 possible gantry angles distributed uniformly with $`5^{}`$ interval for each couch angle. In principle, finer resolution for the couch angle may be used. However, this will also increase the memory requirement since our BAO system, as it is currently implemented, requires the kernel for each beam to be available for quick access in the system memory. For this HN case, the size of the kernel is about 3 MB for each beam and therefore it requires about 650 MB of memory to store all the kernels. We have used a 750-MHz Sun workstation to perform the BAO for this HN case. The 1000-iteration BAO run for the 9 non-coplanar beams required 17 hours to complete (this includes about 2 hours of computing time required for kernel extraction) while the 1000-iteration run for the 7 non-coplanar beams required 13 hours (kernels were read from disk for this case). The relatively long optimization time for this HN case, compared to the corresponding time for the prostate case, is due to the substantially larger size of the system. A side-by-side comparison of these two systems is presented in table 1. Note that, as pointed out previously in the discussion regarding figure 3, the number of iterations used can be substantially reduced without significant reduction in the optimal score. This has the potential of shortening the total optimization time by a factor of two. Further reduction in time can be obtained from using a more selective search, e.g. not allowing the beams to be too close to each other, inclusion of prior knowledge, and utilization of parallel processors. These possible avenues for improvements are currently being investigated.

Figure 6 compares the DVHs for a set of 9 non-coplanar BAO beams with those of a 9-beam equispaced coplanar set. The BAO set provides both better dose uniformity at the GTV and CTV and better protection of the parotid glands. This better sparing of the parotid glands can also be seen in the isodose contours shown in figure 7. Figure 8 shows that a plan of similar quality can also be obtained using 7 non-coplanar beams, although in this case the quality of target coverage is lower than that of the 9 non-coplanar beams. The dose prescriptions for these plans are shown in table 6. The plans have been normalized to best satisfy all of these prescriptions. These have been achieved by normalizing the DVHs to 60 Gy at 95% of the CTV volume.

## 4 Discussion

### 4.1 Dependence on criteria

A noteworthy feature of a DV-based optimization is that further improvement will no longer be attempted by the optimizer once the criteria are satisfied. In other words, once the contribution of a DV criterion to the OF reaches zero, it becomes irrelevant to the optimizer. Consider, e.g., a single criterion which is assigned to the rectum: $`<6000`$ cGy at 5% of its volume. Once this criterion is satisfied, e.g., when $`D_5`$(rectum) has been lowered to 5900 cGy, the optimizer will no longer attempts to improve the dose to the rectum, even though it is still possible to reduce its dose further.

We have found that the improvement on plan quality that is output by our BAO system generally depends on the DV criteria that we assign to the VOIs. Figure 9 compares two of the top 5 plans generated by the BAO (for the prostate case presented in section 3.1) using several “soft” criteria listed in table 7, with the corresponding plan for a set of 3 equispaced beams. These plans have been generated for the same prostate case as that used in figures 4 and 5. The BAO sets are both predicted, by the IO-TLP method which is used during the BAO, to have lower score compared to the equispaced beams (55 and 75 for plans shown in figure 9a and 9b, respectively, compared to 186 for the equispaced beams). However, upon reoptimization using the IO-SC method, the BAO sets produce higher score than the equispaced beams (186 and 212 for figure 9a and 9b, respectively, compared to 162 for the equispaced beams).

The large score increase encountered in going from the IO-TLP to the IO-SC method for the plan shown in figure 9a is largely due to the additional score accumulated by the DVH of the rectum. In this case, the IO-TLP method minimized the score by creating a steep gradient of the DVH curve just below the 6000 cGy dose. Note that even though below 6000 cGy the BAO DVH for the rectum is higher, and therefore is worse, than the corresponding DVH of the equispaced beams, the contribution of this DVH to the OF can be smaller if its DVH at 6000 cGy is lower than the corresponding value for the equispaced beams. This is specifically what happens in this case due to the steep gradient of the BAO DVH curve. In going from the IO-TLP to the IO-SC method, however, this steep gradient cannot be sustained because of the inclusion of more scatter components in the IO-SC method. This “dose leakage” substantially increases the partial volume of the rectum that receives dose greater than 6000 cGy and, consequently, produces a large increase in score.

The BAO plan in figure 9b satisfies the criteria listed in table 7 as well as, or indeed better than, that of the equispaced beams. In this case, the higher score of the BAO plan is due to the contribution from the “normal tissue” (NT) region. The NT region encompasses generic normal tissue surrounding the target that does not belong to any particular OAR (Wu and Mohan 2000). For the prostate case that we use in this work, the NT is a 1 cm extension of the target that is not part of any OAR and it is assigned an upper constraint of 6000 cGy at 5% of its volume. The “dose leakage” from the target to its surrounding NT, for the case shown in figure 9b is substantially larger than the allowance set by the assigned DV criterion and it therefore introduces a significant contribution to the OF on going from the IO-TLP to the IO-SC method. This is due to the small angular separation of two of the 3 beams in the BAO set ($`238^{}`$ and $`254^{}`$) that provides an elongated dose distribution along the direction of these beams.

In contrast to this relatively large uncertainty in the prediction of the IO-TLP method for cases with “soft” DV criteria, in our experience we found that its prediction is much more reliable for cases with more demanding criteria, such as listed in table 3 and displayed in figures 4 and 5. Although we only present this dependence on criteria using a prostate-case example in this paper, similar behavior has also been observed for the HN case. One possible reason that may be given for this observation is that the degeneracy of the solution, i.e. the multiminima of the OF, is larger for “soft” criteria (Wu and Mohan 2002). A careful study, however, is necessary in order to provide a definitive answer to this question. From the point of view of clinical applications, the fact that BAO is more useful for “hard” cases is not expected to eliminate its practical relevance. It is well known that “soft” cases can be planned satisfactorily with equispaced beams or, indeed, class solutions. On the other hand, for “hard” cases, e.g. those of re-irradiation or of complicated geometry, class solutions are non-existent (almost by definition) and equispaced beams may not provide the optimal solution. For these cases the fast BAO described in this paper provides an attractive alternative to a manual trial-and-error approach.

### 4.2 Further speed up?

Several schemes have been proposed recently to rank beam directions according to their individual potential for creating a desirable IMRT distribution such as the single-beam cost function (Oldham et al1998), beam’s-eye-view dosimetrics (BEVD) (Pugachev and Xing 2001), and the figure-of-merit parameter (Chuang et al2003). These schemes have also been studied as a possible way to intelligently accelerate a BAO. If a particular beam direction can be reliably predicted to always give rise to high values of objective function, independent of the directions of the other beams in the set it is a member of, it can be quickly removed from the search space. It is perhaps important to point out here, in connection with the OF of the optimization, that these schemes are expected to have less significance in a DV-based optimization than in a dose-based optimization for which they have been proposed and studied. Individual ranking of beam directions is useful in a dose-based optimization where each VOI is assigned a specific prescribed dose. In contrast, in a DV-based optimization each VOI is effectively assigned a specific range of acceptable doses. In this case, several different distributions of the total dose may have equally good value of the DV OF (in other words, the OF possesses multiple local minima and degenerate solutions). It is therefore less meaningful to attach a specific geometric “shape”, and therefore potentially “good” or “bad” beam directions for creating this shape, to the optimal distribution.

No a priori knowledge has been used to guide the selection of beam angles in this work. The use of reliable prior knowledge is expected to further improve the speed and potential of our BAO system. However, the complexity of the BAO problem requires judicious selection, and therefore careful systematic study, of what “knowledge” is best to be used for each given case. As has been discussed in this paper, this systematic study is especially more pertinent to the DV-based IMRT optimization that we use in our BAO system, because of the presence of the multiple local minima of the OF, than to the dose-based optimization that is commonly used in other BAO systems. In our BAO system, the beam angles are selected from the preset values in the search space. Some prior knowledge can be incorporated in the initial selection of beams that are to be included in the search space, e.g., more and less beams can be positioned in favorable and unfavorable directions, respectively. Alternatively, the beam angles selected during the optimization can be restricted to have a minimum separation from each other. Prior knowledge can also be used to reduce the number of beams in the search pool. Although large search pools are used in this study, 180 and 216 beams for prostate and HN case respectively, in a clinical setting one may not need to use such large numbers. Using a reasonable spacing of 10 between adjacent beams in a coplanar pool, e.g., allows us to reduce the size of the search pool to 36. This will significantly decrease the number of possible combinations of beam angles to be sampled during the BAO. In the current version of our BAO system, each beam in the trial set is given a uniform initial profile prior to the beam-profile optimization. Instead of this uniform intensity distribution, the initial beam profile can be set to a non-uniform one, chosen based on some prior knowledge. This has the potential of reducing the number of iterations needed for the beam-profile optimization and, consequently, the total BAO time.

## 5 Conclusion

In this paper we have demonstrated that a combination of a fast table-lookup dose calculation method and the FSA algorithm can be used in a practical BAO system. This method has the advantage of being able to produce clinically significant improvement within a reasonable time. We showed, e.g., that the beam angles of a 7-beam prostate plan can be optimized within about 2 hours on a moderate 400 MHz workstation. A DV-based OF is used to guide the IMRT optimization in our BAO system. This provides a direct link to our clinical IMRT system and facilitates the comparison between the BAO plans and those produced by human planners. Since the beam-angle optimizer in our system simply wraps around the IMRT beam-profile optimizer, other OFs that are available to our IMRT optimizer, such as EUD, can also be used for BAO in our system. Our BAO system can also search for the optimal set over non-coplanar beams. This feature may be particularly helpful to planners because of the insufficient experience with non-coplanar beams in treatment planning.

Acknowledgements

This work is supported in part by grant CA74043 from the National Cancer Institute and a research grant from Philips Radiation Oncology Systems.

## References


format: text

text:
# Large Scale Structure in the ELAIS S1 Survey

## 1 Introduction

Theories of structure formation were strongly constrained by the statistical measurements of clustering in some of the early galaxy redshift surveys. Surveys of infrared galaxies, in particular, were able to rule out the then standard Cold Dark Matter Model (Efstathiou et al., 1990; Saunders et al., 1991). Present day redshift surveys such as the 2dFGRS (Colless et al., 2001), SDSS (York et al., 2000) and, in the far-infrared, the Point Source Catalog Redshift survey, PSC-z (Saunders et al., 2000) are now able to provide definitive measurements of the galaxy clustering in the local Universe.

Despite this success, we have always known that galaxies are biased tracers of the matter distribution and yet we have a poor observational or theoretical understanding of this bias, although it is assumed to be related to the process of galaxy formation and evolution. To understand bias and, by inference, galaxy formation, we need to better understand the clustering of different galaxy types and the evolution of this clustering with redshift.

In this paper we attempt to provide an estimate of the clustering of infrared galaxies a factor of ten deeper (in redshift) than those seen in the IRAS surveys. To do this we provide the first estimate of clustering from any of the extragalactic ISO surveys. This is thus the first estimate of clustering from galaxies selected at 15$`\mu `$m. We have used part of the ELAIS survey (Oliver et al. 2000) as this probes the largest volume of any of the ISO surveys. We measure the projected clustering by calculating the angular correlation function, we then discuss the constraints this places on the three dimensional clustering using Limber’s equation.

## 2 The ELAIS S1 Survey

The European Large-Area ISO survey (ELAIS, Oliver et al. 2000) was the largest Open Time programme on ISO. This project surveyed 12 square degrees, larger than all but the serendipitous surveys, making it ideal for clustering studies. The main survey bands were at 6.7, 15, 90 and 170 $`\mu `$m. Of these bands the 15$`\mu `$m catalogues contain the greatest density of galaxies (see e.g. Gruppioni et al., 2002; Serjeant et al., 2000), and provide the best statistics for clustering. The final analysis of the 15$`\mu `$m data using the Lari method for one of the ELAIS fields (S1) has recently been completed (Lari et al., 2001) and this is the sample that we use in this analysis.

The S1 field is located at $`\alpha (2000)=`$ $`00^h34^m44.4^s`$, $`\delta (2000)=`$ $`43^{}28^{}12^{\prime \prime }`$, covering an area of $`2\times 2`$ square degrees. The 15µm survey is made from 9 raster observations, each one of $`40^{}\times 40^{}`$. The central raster S1\_5 has been observed three times. Using the Lari method we have obtained a sample of 462 sources to $`5\sigma `$ in the flux range 0.45-150 mJy (Gruppioni et al., 2002).

## 3 Selection function

Besides the galaxy catalogue itself, the selection function is the most important ingredient in the calculation of clustering statistics. Errors in the selection function will invalidate the answer, whereas errors in the weighting scheme will usually make the answer more noisy.

A selection function is required for each source list that is being investigated. The selection function, $`\varphi `$, is defined as the expected number density of sources as a function of $`𝐫`$ (which might be two or three dimensional), in absence of clustering; i.e., the expected number of galaxies $`d𝒩`$ in a volume $`dV`$ is $`d𝒩=\varphi (𝐫)dV`$. With this definition, $`\varphi (𝐫)𝑑V=𝒩`$. The selection function is used to simulate a catalogue with no clustering.

To be selected from the ELAIS S1 catalogue sources had to exceed a signal-to-noise threshold, $`\sigma _{min}`$. The signal to noise of a detected source $`i`$, is $`\sigma _i(𝐫_𝐢)=S_i/N(𝐫_𝐢)`$ where $`S_i`$ is the signal of the source and $`N(𝐫_𝐢)`$ is the noise at the position of the source. Had this source been in a different part of the survey, $`𝐫`$ it would have had a different signal-to-noise. We can define then a mask, $`M_i(𝐫)`$, which represents the detectability of each object as a function of position as follows

$$M_i(𝐫)=\{\begin{array}{cc}0\hfill & \mathrm{if}\sigma _i(𝐫)<\sigma _{min}\hfill \\ 1\hfill & \mathrm{if}\sigma _i(𝐫)\sigma _{min}\hfill \end{array}$$
(1)
where $`\sigma _{min}=5`$. The un-normalised selection function can then be written as

$$\varphi ^{}=\underset{i}{}M_i(𝐫)$$
(2)
which can be normalised

$$\varphi =\varphi ^{}\frac{𝒩}{\varphi ^{}𝑑V}$$
(3)
### 3.1 Building the masks

In the full ELAIS S1 region there are 9 independent noise maps $`N(𝐫)`$, corresponding to 9 independent sub-catalogues. Note that the central noise map is less noisy and the corresponding sub-catalogue deeper, because the ISO data were already combined (Gruppioni et al., 2002). We constructed a selection function as follows: for each source in the sub-catalogue we calculate the hypothetical signal-to-noise ratio (defined as the peak flux over the rms value) at each point in the raster. Where these exceed the extraction signal-to-noise threshold $`\sigma _{min}`$ (equation 1), the value of the selection function at that position is incremented (equation 2).

The 9 individual selection functions are then combined into a single one. Figure 1 shows the final image. In the overlap region only one selection function was used and the final catalogue excludes sources in that region that arose from the other sub-catalogues. Sources with stellar counterparts have also been removed (see Gruppioni et al., 2002) from the catalogue and excluded from the calculation of the selection function. We end up with a catalogue of 329 sources.

The selection function so obtained is then used to generate the random catalogues with no clustering, essential to properly calculate the two point correlation function.

## 4 The Angular Correlation Function

Correlation functions are widely used to study the distribution of sources in surveys and to derive large scale properties of galaxies. The two-point spatial correlation function is defined so that

$$dP=n^2[1+\xi (r)]dV_1dV_2$$

is the joint probability of finding a source in a volume element $`dV_1`$ and another source in a volume element $`dV_2`$. The function $`\xi (r)`$ is the excess probability of finding an object compared to a random distribution of objects.

Similarly, one can define the two point angular correlation function so that

$$dP=n^2[1+w(\theta )]d\mathrm{\Omega }_1d\mathrm{\Omega }_2$$

is the joint probability of finding a source in a solid angle element $`d\mathrm{\Omega }_1`$ and another source in a solid angle element $`d\mathrm{\Omega }_2`$. These two statistics are related by Limber’s equation (Peebles 1980).

In order to calculate the angular correlation function of mid-IR sources we use the Landy & Szalay estimator (Landy & Szalay, 1993)

$$w(\theta )=\frac{[\mathrm{DD}]2[\mathrm{DR}]+[\mathrm{RR}]}{[\mathrm{RR}]}$$
(4)
where \[DD\], \[RR\] and \[DR\] represents the normalized number of galaxy-galaxy, random-random and galaxy-random pairs with angular separation in $`(\theta ,\theta +d\theta )`$.

Errors in the calculation of the angular correlation function are dominated by Poisson noise. The error in each bin can be estimated using the following expresion (Baugh et al., 1996b):

$$\delta w(\theta )=2\sqrt{\frac{1+w(\theta )}{DD}}$$
(5)
where, in this case, $`DD`$ is the total number of galaxy-galaxy pairs (not normalized). Errors calculated using this equation are comparable to the errors obtained from a bootstrap resampling technique (Ling, Barrow, & Frenk, 1986).

A second source of errors comes from the to finite size of the sample. In order to correct for this effect we use the random sample to calculate the integral constraint as (e.g. Infante et al. 1994)

$$\delta =\frac{N_{rr}(\theta )}{N_{rr}(\theta )(1+w(\theta ))}$$
(6)
and divide the calculated correlation function by this factor, $`\delta =0.945`$.

Figure 2 (top) shows the obtained angular correlation function, calculated using 200 realizations with 2000 random sources each, in intervals of $`\mathrm{log}\theta =0.08`$ degrees. Random catalogues have been built using the selection function calculated in previous section. Data points have been fitted by an exponential law of the form $`w(\theta )=A\theta ^{1\gamma }`$ resulting in $`A=0.014\pm 0.005`$ and $`\gamma =2.04\pm 0.18`$ (where $`\theta `$ is measured in degrees).

The mean number of objects in each field is $`n=35.25`$ (excluding the central raster S1\_5, with 71 sources), with a standard deviation of 6 objects. Since S1\_5 field reaches deeper flux limits, it could in principle be subject to clustering variations that would affect the whole clustering estimation. We repeat the above calculation removing those sources with fluxes fainter than the flux limit excluding S1\_5: a total of 27 sources with fluxes lower than 0.7 mJy are removed. By calculating again the selection function and the angular correlation function we obtain then an estimate of the clustering from sources detectable over most of the S1 field. In this case the integral constraint is slightly larger than the previously obtained, $`\delta =0.970`$. The angular correlation function is shown in figure 2 (bottom) and is fitted by the same exponential law previously calculated. Although the correlation function is now slightly larger at scales $`0.1`$ deg, the errors are also larger and the overall correlation function is well fitted by the previous function.

## 5 Spatial correlation function

In case of small angles, both $`w`$ and $`\xi `$ can be approximated by power law shapes, and the spatial correlation function can be written as (e.g., Phillipps et al., 1978)

$$\xi (r,z)=\left(\frac{r}{r_0}\right)^\gamma (1+z)^{(3+ϵ)}$$
(7)
where $`r_0`$ is the comoving correlation length at $`z=0`$ and $`r`$ the comoving distance. The parameter $`ϵ`$ is the clustering evolution index and is interpreted as follows. A value $`ϵ=0`$ corresponds to stable clustering in physical coordinates, i.e., galaxy clusters remain unchanged and clustering changes due to the expansion of the Universe, while $`ϵ=3\gamma `$ corresponds to clustering fixed in comoving coordinates, i.e., clustering does not change with time and the galaxy clusters expand with the Universe.

If $`w(\theta )`$ is parametrized as $`w(\theta )=A_w\theta ^{1\gamma }`$, then Limber’s equation becomes (e.g., Phillipps et al., 1978)

$$A_w=Cr_0^\gamma \frac{D_\theta ^{1\gamma }g^1(z)(1+z)^{(3+ϵ)}(dN/dz)^2𝑑z}{((dN/dz)𝑑z)^2}$$
(8)
where $`D_\theta `$ is the angular diameter distance, $`g(z)`$ is the scale factor multiplied by the element of comoving distance

$$g(z)=\frac{c}{H_0}[(1+z)^2(1+\mathrm{\Omega }_0z)^{1/2}]^1$$
(9)
and

$$C=\pi ^{1/2}\frac{\mathrm{\Gamma }[(\gamma 1)/2]}{\mathrm{\Gamma }(\gamma /2)}$$
(10)
The only unknown quantity in equation 8 is the redshift distribution of the sources. A distribution given by

$$\frac{dN}{dz}=\frac{3𝒩\mathrm{\Omega }_s}{2z_c^3}z^2\mathrm{exp}\left[\left(\frac{z}{z_c}\right)^{3/2}\right]$$

where $`z_m=\sqrt{2}z_c`$ is the median redshift of the survey, provides good fits to the distribution of optical galaxies (Baugh, 1996a, and references therein) and as been widely used to invert Limber’s equation. However, the redshift distribution of mid-IR sources is more extended and the equation above is no longer valid. Figure 3 (top) shows the redshift distribution of the ELAIS sources, obtained from optical spectroscopy (La Franca et al, 2004, in prep; Perez-Fournon et al, 2004, in prep) and photometric redshifts (Rowan-Robinson et al., 2003). Instead of the optical $`dN/dz`$ we use

$$\frac{dN}{dz}z\mathrm{exp}\left[\left(\frac{z}{z_c}\right)^{3/4}\right]$$
(11)
which fits reasonably well the distribution of ELAIS sources.

We can now integrate equation 8 in order to obtain $`r_0`$. Assuming $`\mathrm{\Omega }_0=1.0`$ and $`ϵ=0`$ and using a median redshift of $`z_m=0.1`$ in equation 11 , we obtain $`r_0=4.3h^1`$ Mpc with a 95% confidence level in the range $`3.8h^1`$ Mpc to $`4.7h^1`$ Mpc.

An alternative redshift distribution of the ELAIS sources has been presented by Pozzi et al. (2004), computed from the luminosity function fit of galaxies on ELAIS S1 and S2 (see figure 3 – bottom). When using this redshift distribution as input to equation 8 we obtain similar results of $`r_0=4.2h^1`$ Mpc with a 95% confidence level in the range $`3.6h^1`$ Mpc to $`4.8h^1`$ Mpc.

## 6 Discussion

We use a sample of 329 15µm galaxies detected in the ELAIS S1 survey covering a region of 4 square degrees of sky to determine the angular correlation function of the galaxies. We measure $`w(\theta )`$ up to scales of 1 degree, corresponding to $`7`$ Mpc at the mean redshift of $`z=0.14`$. The resulting correlation function is well fitted by a power law $`w(\theta )=A\theta ^{1\gamma }`$ with $`A=0.014\pm 0.005`$ and $`\gamma =2.04\pm 0.18`$. Assuming a redshift distribution we invert the angular correlation using Limber’s equation in order to determine the correlation length $`r_0`$. We find a value of $`r_0=4.3h^1`$ Mpc with a 95% confidence level in the range $`3.8h^1`$ Mpc to $`4.7h^1`$ Mpc.

Table 1 lists the correlation lengths obtained from several optical and infrared surveys. While optical surveys show a significant clustering with $`r_05h^1`$ Mpc, IRAS survey shows a lower value. The data described in this paper is consistent with this trend: that mid infrared selected sources show a smaller correlation length. This is expected because optical surveys favour elliptical galaxies which are more strongly clustered than spiral galaxies, while infrared surveys preferentially select spirals and star-forming galaxies. We note that the ELAIS selected galaxies appear to have a marginally steeper two point correlation function than the optical and the IRAS surveys. By fixing $`\gamma `$ to the lower value allowed by the fit $`\gamma =1.86`$ we obtain $`r_0=4.1_{0.5}^{+0.2}`$, bringing the clustering amplitude closer to that seen by IRAS.

It is interesting to note that the optical correlation function at $`z0.2`$ (Zehavi et al., 2002; Norberg et al., 2001) is consistent with that at $`z0`$ (Maddox et al., 1990), i.e. there is no apparent evolution in the correlation function over this albeit small redshift range. Likewise the infrared galaxy correlation function estimate from this work at $`z0.14`$ is consistent with the infrared galaxy correlation function estimated from IRAS galaxies at $`z=0`$ (Jing et al., 2002). This apparent lack of evolution may be because evolutionary effects are small or that evolution in mass clustering is compensated by evolution in galaxy bias. It will be interesting to see if this apparent non-evolution in clustering of different galaxy population mixes is still found in deeper surveys (e.g. SWIRE, Lonsdale et al. 2003) as this might imply a striking conspiracy in the evolution of the bias in different galaxy types.

By performing a systematic analysis of all density peaks in the redshift distribution of field galaxies, Elbaz & Moy (2004) recently found an excess of ISO selected galaxies over the whole range of density peaks. This suggests than infrared galaxies are more strongly clustered than optical galaxies, in apparent contradiction to our results. Since the ISO-CAM surveys are deeper than ELAIS, it is plausible that there is an evolution of the clustering over the redshift interval spanned by these two ISO surveys. Infrared galaxies would then become more clustered towards higher redshift, while the clustering of optical galaxies changes little. This would agree qualitatively with hierarchical pictures of structure formation. Such models predict that star formation and galaxy formation is driven by merger rate that would be a function of environment. Star formation would thus have occurred first and vigorously in the denser (more clustered) regions of the Universe, taking more time to initiate in lower density (less clustered) regions (see also Elbaz & Cesarsky, 2003). The strong evolution in luminosity function of infrared galaxies (e.g. Pozzi et al. 2004) might then be coupled with an evolution in their clustering. Optically selected galaxies sample regions where past star formation activity was high as well as those where current activity is high and are thus less sensitive to these effects. Finally this apparent contradiction may simple be that the deep ISO-CAM surveys do not sample a sufficiently large volume to be representative of the rest of the Universe.

## 7 Acknowledgements

Eduardo Gonzalez-Solares was supported by PPARC grant PPA/G/S/2000/00508 and EC Marie Curie Fellowship MCFI-2001-01809. This paper is based on observations with ISO, an ESA project, with instruments funded by ESA Member States (especially the PI countries: France, Germany, the Netherlands and the United Kingdom) and with participation of ISAS and NASA. Some of this work was supported by EEC Training Mobility Research Network “Probing the Origin of the Extra-galactic background light” (POE) HPRN-CT-2000-00138. The author also wants to thank D. Elbaz for useful comments on the discussion of the results presented.


text:
# Detection of Sagittarius A∗ at 330 MHz with the Very Large Array

## 1. Introduction

The central radio-bright region of our Galaxy, known as the Sagittarius Complex, is composed of three major components: The supernova remnant (SNR) Sagittarius A East (Sgr A East (catalog )), the Sagittarius A West (Sgr A West (catalog )) H II region, and Sagittarius A (Sgr A (catalog )), recently established as our Galaxy’s central massive black hole (e.g., Ghez et al. 2000; Eckart et al. 2002).

Models attempting to explain the emission from Sgr A (catalog ) fall into three broad classes. Emission is modeled as arising from thermal sources, such as a low-temperature accretion disk, from non-thermal sources such as a jet (e.g., Melia & Falcke, 2001 and references therein), or from a mixture of the two. Such models are constrained primarily by the observed spectrum of Sgr A (catalog ), but large gaps in frequency coverage exist. For this reason, filling in such gaps, as with the recent near-IR detections (Ghez et al. 2003, submitted; Genzel et al. 2003), and extending the range of frequencies over which the source is observed is important in order to place additional observational constraints on these models.

In 1976, Davies et al. attempted to observe Sgr A (catalog ) at 410 MHz, but reported no detection above a level of 50 mJy. The first Very Large Array (VLA) images of the Galactic Center (GC) at 330 MHz (Pedlar et al. 1989, Anantharamaiah et al. 1991) represented a major improvement in sensitivity and resolution over previous meter wavelength images. Pedlar et al. did not detect Sgr A (catalog ) down to a 5$`\sigma `$ level of 100 mJy. Though Anantharamaiah et al. had higher sensitivity, (rms $`3`$ mJy beam-1) they also reported a non-detection. The absence of any detections below 1000 MHz was interpreted as thermal absorption obscuring Sgr A (catalog ), indicating the low frequency properties of Sgr A (catalog ) are dominated by extrinsic absorption, not intrinsic emission.

New low frequency detections of Sgr A (catalog ), including the first detection of the source below 1 GHz (610 MHz; Roy et al. 2003), and this 330 MHz detection suggest that the low frequency properties of Sgr A (catalog ) and the implied line of sight optical depths need to be re-examined.

## 2. Observations and Data Reduction

The Galactic Center was observed at 330 MHz using the Very Large Array in its A configuration; Table 1 summarizes the observations. The visibility amplitudes were calibrated with observations of 3C 286 (catalog 3C), while the visibility phases were calibrated with observations of the VLA calibrators J1751$``$253 (catalog VLA) and 1714$``$252 (catalog VLA).

In this letter we will summarize only salient details of the image processing. Full details of the processing including further details on issues such as radio frequency interference excision, astrometry, bandwidth smearing, and extensive details pertaining to calibration and ionospheric compensation will be provided in Nord et al. (2003).

Observations were taken in spectral-line mode, both to reduce the impact of bandwidth smearing and to enable excision of radio frequency interference. Excision was done manually by searching for visibilities with anomalously high amplitudes, relative to surrounding visibilities in the $`u`$-$`v`$ plane as well as by identifying the Fourier components that gave rise to systematic “ripples” in the image. In general, low-frequency observations with the VLA must take the non-coplanarity of the array into account. Here this was not done because the field of view required to image Sgr A East and West is sufficiently limited ($`3\mathrm{}`$) that the array can be assumed to be nearly coplanar. Several iterations of phase-only self-calibration followed by imaging were performed in order to improve the dynamic range of the images and compensate for ionospheric-induced phase distortions. Because data are being combined from different epochs, each epoch was imaged separately to produce an image with an acceptable dynamic range before being combined. The resulting self-calibrated $`u`$-$`v`$ data from the separate epochs were combined and phase and amplitude calibrated against the 1998 image in order to account for any offset in flux density scales and coordinates grids.

In order to obtain maximum resolution, the super-uniform weighting scheme (Briggs, Schwab, & Sramek, 1999) was used. Super-uniform weighting111UVBOX$`=3`$ and ROBUST$`=1`$ within the AIPS task IMAGR. reduces the weight, both globally and locally, given to the central, densely sampled region of the $`u`$-$`v`$ plane, maximizing the resolution at the expense of increased noise and sidelobe level. In this case the increase in noise and sidelobe level is an acceptable trade off, as the analysis is limited more by the $`u`$-$`v`$ coverage and resulting image fidelity than by the thermal noise level. Figure 1 shows the super-uniform image of the Sgr A (catalog ) region.

## 3. Analysis

### 3.1. Predicted Characteristics of Sgr A (catalog ) at 330 MHz

The expected size of Sgr A (catalog ) was calculated from the measurements from 5 to 45 GHz by Bower et al. (2003, in preparation). They find the scattering disk of Sgr A (catalog ) to be $`(1.39\pm 0.02\times 0.65_{0.14}^{+0.11})`$ milliarcseconds at $`\lambda =1`$ cm. Scaling these dimensions by $`\lambda ^2`$ as expected from interstellar scattering gives an expected 330 MHz diameter of $`(11.5\pm 0.2\times 5.4_{1.1}^{+0.9})`$ arcseconds. The position angle of the scattering disk of Sgr A (catalog ) is constant at $`80\mathrm{°}`$ over the range of 1.4 to 46 GHz. We have assumed that its orientation remains constant below $`1.4`$ GHz. Sgr A (catalog ) has a time averaged flux density at $`1.4`$ GHz of $`0.5`$ Jy and a time averaged spectral index over the centimeter radio regime of $`\alpha 0.3`$ ($`S_\nu \nu ^\alpha `$) (Zhao, Bower, & Goss, 2001). Assuming no change in the spectral index below $`1.4`$ GHz, the expected 330 MHz flux density is roughly 0.4 Jy.

### 3.2. The Detection of Sgr A (catalog )

Figure 2 shows slices through the assumed major and minor axes of Sgr A (catalog ). Gaussian fits with baseline subtraction to these slices results in a measured source size of $`(17\stackrel{}{\mathrm{.}}6\pm 4)\times (15\stackrel{}{\mathrm{.}}9\pm 5)`$, a peak intensity of 80 $`\pm `$ 15 mJy beam-1, and a total flux density of 330 $`\pm `$ 120 mJy. The source deconvolved with the $`6\stackrel{}{\mathrm{.}}8\times 10\stackrel{}{\mathrm{.}}9`$ beam gives a source size of $`(16\mathrm{}\pm 4)\times (11\mathrm{}\pm 5)`$. If the source is fit in the image instead of fitting slices,, the position angle is then a free parameter instead of assumed to be 80°. A 2 dimensional Gaussian fit of the region yields a position angle of $`35\mathrm{°}\pm 35`$. However, this fit is contaminated by flux density to the south of the source as shown in the slice on the right in Figure 2.

The location and flux density of the source agree with the extrapolations of §3.1, and the size of the source agrees to within $`1.2\sigma `$ along the assumed minor axis and $`1.1\sigma `$ along the assumed major axis. However, an alternate interpretation of this emission is that it originates from the Sgr A East SNR.

The absorption in Figure 1 is caused by the Sgr A West H II region, first detected in absorption by Pedlar et al. (1989). The Sgr A West H II region is composed of a dense region known as the “mini-spiral” (Lo & Claussen 1983, Ekers et al. 1983), known to be orbiting Sgr A (catalog ) (Roberts and Goss 1993), and a more diffuse ionized region $`80\mathrm{}\times 60\mathrm{}`$ in size, known as Sgr A West (extended) (Mezger & Wink 1986; Pedlar et al. 1989). Both components of Sgr A West are detected in thermal (free-free) absorption against the Sgr A East SNR with 330 MHz optical depths ($`\tau _{330}`$) $`>1`$ across a large portion of the H II region (Pedlar et al. 1989). The position of Sgr A (catalog ) is offset from the mini-spiral (Zhao & Goss 1998), therefore in order for Sgr A East to be absorbed at the position of Sgr A (catalog ) most of the absorption must arise from Sgr A West (extended). Estimates of the optical depth along lines of sight near Sgr A (catalog ) based on higher frequency measurements ($`\nu =5`$ GHz, Yusef-Zadeh & Morris 1987) suggest that 330 MHz optical depths near Sgr A (catalog ) are $`\tau _{330}10`$. Moreover, any gap through the Sgr A West H II region would have to occur with the size, shape, and location necessary to allow Sgr A East to be detected with the expected size, shape, and flux density of Sgr A (catalog ). Given the agreement in observed and predicted source properties, the likelihood that Sgr A East is detected can be discounted. We conclude that we have detected Sgr A (catalog ) at 330 MHz.

## 4. Implications

### 4.1. Previous Non-detections

Previous non-detections at similar frequencies can be attributed to the limited sensitivities of the observations. Davies et al. (1976) observed with the Mk I-Defford interferometer. This two-element interferometer has a north-south resolution of $`\theta 2\stackrel{}{\mathrm{.}}4`$ at 408 MHz. The correlated flux density that would have been measured by this interferometer is $`S_{\mathrm{corr}}=S_0e^{(\theta _s/\theta )^2}`$, where $`\theta _s`$ is the scattering diameter and $`S_0`$ is the intrinsic flux density (Roy et al. 2003). The expected scattering diameter of Sgr A (catalog ) is 4$`\stackrel{}{\mathrm{.}}`$1. Their upper limit, $`S_{\mathrm{corr}}<50`$ mJy, implies that they would not have detected any source weaker than 0.9 Jy, a value that is well above the flux density that we have observed at 330 MHz.

Pedlar et al. (1989) observed the Galactic center with the VLA. At the time of their observations, only a few ($`10`$) of the 27 VLA antennas were outfitted for 330 MHz observations. In addition to a higher thermal noise level (rms sensitivity $`16\text{mJy beam}\text{-1}`$), the number of baselines was substantially lower, resulting in less complete $`u`$-$`v`$ coverage. Sgr A (catalog ) would have been only a 3$`\sigma `$ source in a crowded and confusing field.

In 1999, Anantharamaiah, Pedlar, & Goss reanalyzed data from Anantharamaiah et al. (1991) in which the GC was observed at 330 MHz with all antennas of the VLA. In their Figure 2, there is a hint of a detection of Sgr A (catalog ) at a level of $`40`$ mJy beam-1, which is consistent with this detection. The authors do not comment on this possible detection.

### 4.2. The Spectrum of Sgr A (catalog )

Figure 3 shows the observed spectrum for Sgr A (catalog ), including our flux density measurement. This spectrum is not constructed from simultaneous measurements and as Sgr A (catalog ) is variable, simultaneous measurements may show differences.

There are at least two models of the low frequency spectrum of Sgr A (catalog ) that are suggested by our detection. Either the $`S_\nu \nu ^{0.3}`$ spectrum observed across the centimeter radio band continues toward low frequencies and the source is unobscured by intervening thermal gas, or the spectrum of Sgr A (catalog ) rises below $``$ 1 GHz and is obscured by significant optical depth.

Assuming a $`S_\nu \nu ^{0.3}`$ spectrum, a 330 MHz optical depth of $`\tau _{330}<0.1`$ is required to explain the measured Sgr A (catalog ) flux density. Even a flat spectrum below 1.4 GHz permits an optical depth of only $`\tau _{330}0.4`$. Therefore, the free-free optical depth along this line of sight is much less than the $`\tau _{330}`$ 3–5 previously assumed to explain low frequency non-detections (Pedlar et al. 1989).

A conclusion of §3.2 is that lines of sight to Sgr A East that pass near Sgr A (catalog ) have large 330 MHz optical depths. If the line of sight to Sgr A (catalog ) were similarly obscured, its intrinsic flux density would have to be increasing at low frequencies in order to agree with the observed flux density. The optical depth has a frequency dependence of $`\nu ^{2.1}`$. In order to agree with the observed flux density, the intrinsic flux density would have to rise exponentially ($`S_\nu e^{\nu ^{2.1}}`$) below $``$ 1 GHz.

Low 330 MHz optical depths near Sgr A (catalog ) are also indicated by the 22 GHz, $`0.1\mathrm{}`$ resolution image of the Sgr A (catalog ) region by Zhao and Goss (1998), which shows emission measures due to ionized gas within $`0.2\mathrm{}`$ of Sgr A (catalog ) are nearly two orders of magnitude lower than in the nearby mini-spiral. Furthermore, an exponential rise in intrinsic flux density to almost exactly cancel the large (factor of $`10`$ between 1000 and 330 MHz) increase in free-free optical depth at 330 MHz would be quite improbable .

We conclude that it is likely that there is little or no free-free absorption along the line of sight to Sgr A (catalog ), and that the observed flux density is intrinsic to the source. Our results imply a large optical depth towards the Sgr A East SNR, but a small optical depth towards Sgr A (catalog ). This could be explained by a localized clearing of the ambient gas accomplished either through the direct influence of the black hole, or through the stellar winds of the central cluster. It is also possible that the ionized medium in the GC region is clumped, with the line of sight to Sgr A (catalog ) passing through a gap. Our data cannot differentiate between these possibilities.

### 4.3. Emission Mechanisms

Previous modeling of the spectrum of Sgr A (catalog ) did not have to confront its flux density at low frequencies, as Sgr A (catalog ) was thought to be unobservable below 1 GHz (e.g. Mahadevan 1998). Furthermore, the 410 MHz non-detection of Davies et al. (1976) is often incorrectly treated as an upper limit. Our 330 MHz result and the 610 MHz detection by Roy et al. (2003) require these models to be reevaluated. Because we cannot exclude the possibility of free-free optical depth toward Sgr A\*, our reported 330 MHz flux density should be considered a *lower* limit.

Advection dominated accretion flow (ADAF) models (e.g. Yuan, Markoff & Falcke 2002) predict flux densities nearly two orders of magnitude lower than the 330 MHz measurement. The Relativistic Jet plus ADAF model of Yuan et al. (2002) does not explicitly make predictions below $`1.4`$ GHz, but would appear to closely match the measured 330 MHz flux density. It should be noted that in the jet models, low frequency emission arises from farther out along the jet, probing regions farther from the event horizon.

## 5. Conclusions

We have presented a new high resolution, high sensitivity image of the Sgr A (catalog ) region at 330 MHz, constructed from observations at multiple epochs with the Very Large Array. We report a detection of Sgr A (catalog ), with a flux density of $`330\pm 120`$ mJy. Previous non-detections of Sgr A (catalog ) at comparable frequencies are attributable to the poor sensitivity and/or $`u`$-$`v`$ coverage of the interferometers used.

We argue that only small amounts of free-free absorption ($`\tau _{330}0.1`$) exist along the line of sight to Sgr A (catalog ), and therefore that Sgr A (catalog ) is observed through a region of low density in the Sgr A West (catalog ) H II region. It is possible that the intrinsic spectrum of Sgr A (catalog ) rises below $`1`$ GHz. However, the intrinsic rise in flux density would have to be balanced by free-free absorption to result in the observed value. Therefore, this flux density measurement is a *lower* limit to the intrinsic flux density of this source as the optical depth along the line of sight remains unknown (e.g., Anantharamaiah et al. 1999).

This detection at 330 MHz, combined with the recent detection at 610 MHz (Roy et al. 2003), expands the frequency range over which Sgr A (catalog ) has been detected and allows for low frequency observational constraints on Sgr A (catalog ) emission mechanisms. Future observations at and below 240 MHz with the Giant Metrewave Radio Telescope (GMRT) and the Low Frequency Array (LOFAR; Kassim et al. 2001) may expand the frequency range at which Sgr A (catalog ) is detectable to even lower frequencies.

We thank G. Bower for the use of unpublished results pertaining to the scattering size of Sgr A (catalog ), C. Brogan for help with data reduction and imaging, and A. Cohen and L.J. Rickard for helpful discussions. Basic research in radio astronomy at the NRL is supported by the Office of Naval Research. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.