Calculation of absorbed dose in radiotherapy by solution of the linear Boltzmann transport equations

Over the last decade, dose calculations which solve the linear Boltzmann transport equations have been introduced into clinical practice and are now in widespread use. However, knowledge in the radiotherapy community concerning the details of their function is limited. This review gives a general description of the linear Boltzmann transport equations as applied to calculation of absorbed dose in clinical radiotherapy. The aim is to elucidate the principles of the method, rather than to describe a particular implementation. The literature on the performance of typical algorithms is then reviewed, in many cases with reference to Monte Carlo simulations. The review is completed with an overview of the emerging applications in the important area of MR-guided radiotherapy.

conditions of energy conservation also apply. The change in fluence in the direction of interest can therefore be expressed in terms of the fluence in other directions. This model is then considered for all points of interest in the patient volume, a number of directions of interest and a number of ranges of particle energy. The resulting set of simultaneous equations is solved to yield the fluence distribution in the patient, and hence the absorbed dose.
For a derivation of the LBTE from first principles, see Boman (2007). Jörres (2015) also carries out a theoretical treatment of the problem for the more mathematically inclined reader, and the work of Vassiliev (2017) is an invaluable resource for those wishing to study the subject in detail. The paper by Hensel et al (2006) contains the most intuitive and transparent treatment of the subject in terms of physics and this review initially follows their approach. We make the following definitions: Ω γ a unit normal in the direction of interest, subscripted according to radiation type. r the position of interest. E γ the photon energy of interest. E e the electron energy of interest. ρ c (r) the density of atomic cores at position r. ρ e (r) the electron density at position r.
The photon fluence at position r, with direction Ω γ and energy E γ . Φ e (r, Ω e , E e ) The electron fluence at position r, with direction Ω e and energy E e . σ C,γ E γ , E γ , Ω γ · Ω γ The differential Compton scattering cross-section of a photon travelling initially with energy E γ in direction Ω γ and finally with energy E γ and direction Ω γ . σ C,e E γ , E e , Ω γ · Ω e The differential Compton scattering cross-section of a photon travelling initially with energy E γ in direction Ω γ and giving rise to an electron travelling with energy E e and direction Ω e . σ M E e , E e , Ω e · Ω e The differential Møller scattering cross-section of an electron travelling initially with energy E e in direction Ω e and finally with energy E e and direction Ω e . σ Mott E e , Ω e · Ω e The differential Mott scattering cross-section of an electron travelling with energy E e , initially in direction Ω e and finally in direction Ω e . σ tot M (E e ) The total Møller scattering cross section for an electron travelling initially with energy E e . σ tot Mott (E e ) The total Mott scattering cross section for an electron travelling initially with energy E e . Note that the prime is used to denote initial energy and direction for consistency with other works referred to in this review, e.g. Lewis and Miller (1984), Hensel et al (2006) and Vassiliev et al (2010).

Transport equation
For maximum accuracy, all types of photon interaction should be modelled. However, this increases the complexity of the problem, so it is customary to omit or simplify some of the interactions. For example, the model of Vassiliev et al (2010) includes pair production, but assumes that both of the charged particles are electrons. In that instance, partial coupling of photons and electrons is also considered, whereby photons produce electrons but not vice versa. In order to illustrate the concepts, we confine ourselves to Compton interactions, which for megavoltage beams of low energy are predominant. In Compton interactions, the incoming photon scatters from an atomic electron, resulting in a change of energy and direction (figure 2). The electron is ejected at an angle of up to 90° to the original direction of photon transport. The transport equation for photons is then (Hensel et al 2006): where the integration is over initial photon energy, E γ , and initial photon direction, Ω γ . This is saying physically that the increase in particles travelling in a particular direction with a particular energy (i.e. the component of ∇Φ γ (r, Ω γ , E γ ) in direction Ω γ ) depends on the increase due to photons travelling in a different direction but which Compton scatter into the direction of interest, and the loss due to photons which scatter out of the direction of interest, also due to Compton interactions. The increase in photons is described by the integrals over initial energy and direction and the reduction in photons is given by the total scattering cross section. Total scattering cross section is used because it is simply the loss of photons that we are interested in. The direction that the photons acquire after scattering is not important for the particular point of interest that is being considered: it is their loss from the direction and energy of interest that matters. Note that the left-hand side refers to the gradient of the fluence, while the right-hand side refers to an integral of the fluence, so this is an integro-differential equation, which is not trivial to solve, hence the time taken for the concept to come into routine practice in radiotherapy.

Multiple scattered photons
Solution of equation (1) can be simplified by explicitly considering the multiple scattering events that contribute to the integral term. Hensel et al (2006) describe the division of photon fluence into 0-times, 1-times, 2-times, …, N-times scattered photon fluences, with the scattering of each providing the source for the next: where M denotes the number of scattering events that have been undergone. The variables r, Ω γ and E γ have been suppressed for clarity. Photons have a large mean free path in tissue, so there are only several scattering events in the patient, and consequently N can be chosen as only 2 or 3 (e.g. M = 0…2). In general, an even simpler approach is taken, wherein N is taken as 1, so that the photons are divided up simply into unscattered and scattered fluences. The unscattered photon fluence can then be calculated analytically using ray tracing (Siddon 1985). The collided photon fluence, Φ (coll) γ , is then calculated using equation (1), but with the uncollided photon fluence, Φ (unc) γ as a fixed source (Vassiliev et al 2010): Figure 2. Compton interaction. The photon interacts with a bound electron, which, provided the energy of the photon is high with respect to the binding energy of the electron, can be considered to be free. There are specific kinematics which describe the directions of interaction.
Both the uncollided and collided fluences scatter to provide electron sources in the electron transport equation.

Transport equation
As with the photon equation, the gradient of electron fluence is the difference between increase in electron fluence due to electrons scattered into the direction of interest and loss of electron fluence due to electrons scattered out of the direction of interest. The sources of electrons scattered into the direction of interest are described by four terms which are very similar to one another and are based on the physical scattering processes. Each of these terms is now considered separately.

Compton scattering events
This is of relevance for all photons where the electron is ejected with the energy of interest and in the direction of interest. All initial energies, E γ , and directions, Ω γ are considered, each giving rise to an electron with final energy, E e , and direction, Ω e : The integration is shown as being over 4π steradians of solid angle, but note that electrons are only ever ejected in the forwards direction, so that the integration can in practice be carried out over a reduced range of angles. See figure 2.

Møller scattering events
These are Coulomb interactions between a free electron and the bound electrons in the medium. The bound electrons are considered to be free, so long as the energy imparted to the bound electron during the inelastic scattering event is much higher than the binding energy of the bound electron (figure 3). Then: Ω e · ∇Φ e (r, Ω e , E e ) = ρ e (r)ˆ∞ 0ˆ4πσ M E e , E e , Ω e · Ω e Φ e r, Ω e , E e dΩ e dE e .
The differential scattering cross section in this case is for Møller scattering. Some of the Møller scattering events impart sufficient energy to the secondary electron that these then cause appreciable ionisation of their own. However, the direction of the ejected electron is still in the direction of interest, and the subsequent changes in direction of the Delta ray are separate interactions, so the delta rays can be included in the Møller scattering term of the LBTE.

Mott scattering events
Electrons also scatter elastically, without transferring energy, in the presence of the heavy, positively-charged atomic nuclei (figure 4). This type of scattering accounts for almost all of the changes in direction of the travelling electrons: Note that in this case, there is no integral over energy as the interaction is elastic, so that the final energy of the scattered electron is always the same as the initial energy.

Transport equation in full
The complete transport equation for electrons is (Hensel et al 2006): Physically, this states that the gradient of electron fluence is the difference between increase in electron fluence due to electrons scattered into the direction of interest and loss of electron fluence due to electrons scattered out of the direction of interest. The first three terms have been described in the previous sections, and the last two terms describe the loss of electron fluence due to Møller and Mott scattering. Bremsstrahlung is a relatively insignificant effect, so is here neglected, but if considered, would also be included with the last two terms, and would also form a photon source.

Calculation of absorbed dose
The principle of dose calculation is to solve the above equations for photon and electron transport and then to use the calculated fluence to determine the absorbed dose. The dose can be calculated from the stopping power (energy loss per unit distance travelled): by integrating over all fluence directions (Larsen et al 1997): where ρ (r) is the physical density of the medium. However, the only process to actually deposit dose in the medium at the point of interest is Møller scattering, with or without Delta ray production. All of the other processes transfer energy from the primary particle to a secondary particle, which then continues to transport through the medium. (These then in turn give rise to Møller scattering at some other point.) The absorbed dose is therefore given by Hensel et al (2006): , Ω e · Ω e Φ e r, Ω e , E e · (E e − E e ) dE e dE e dΩ e dΩ e . (10) Figure 3. Møller (inelastic) scattering. An incoming electron interacts with a bound electron by Coulomb scattering, the distance of closest approach governing the scattering angles. If the energy of the incoming electron is high with respect to that of the binding energy of the bound electron, the bound electron can be assumed to be free. The released electron may have sufficient energy to cause appreciable ionisation, in which case it becomes a so-called Delta ray. In other words, the dose deposited is calculated by integrating all Møller scattering events, at all initial energies, all final energies, all initial directions and all final directions. Each of these events deposits energy E e − E e .

Spherical harmonics
As equation (1) is difficult to solve for practical geometries, it is common to employ approximations. These are now described.
In order to conveniently describe the macroscopic differential scattering cross sections, the cross sections are usually expanded in Legendre polynomials, P l (µ 0 ), where µ 0 = Ω · Ω , the cosine of the scattering angle. Hence, describing a general scattering source as: the scattering cross section becomes (Gifford et al 2006, Vassiliev et al 2010: where the 'scattering moments' σ l , are given by: The angular fluence of equation (11) is also expanded in spherical harmonics: where Y l,m Ω are spherical harmonic functions, and φ l,m (r, E ) are the spherical harmonic moments of the angular fluence. The equations given above are exact if the integrals and sums are enumerated fully, but in the interests of computation time and storage, it is normal to truncate the sum such that 0 l L. Then by using the Legendre addition theorem, the complete expression for the scattering source is: The value of this formula in Legendre polynomials and spherical harmonics is that it replaces the integral over orientation in equation (11) with a discrete summation, which is much easier to compute. It also provides a very good analytical approximation to the exact formulae when relatively few terms of the expansions are retained. Retention of only some of the terms allows for a short computation time.

Fokker-Planck approximation
Scattering of electrons is forward-peaked, i.e. the majority of interactions involve only a small change in direction of the incident electron, and a correspondingly small loss in energy. Consequently, it is possible to carry out Taylor expansions of the scattering integrals in the LBTE equations and then retain only the terms relating to small directional changes, with only a small loss in accuracy. This is known as the Fokker-Planck approximation, in which electrons are considered to continuously change direction and energy as they move through the medium (Larsen et al 1997). The derivation of this result (Chandrasekhar 1943, Pomraning 1992) is beyond the scope of this review, but the results are stated here for completeness (Larsen et al 1997, Hensel et al 2006. Equation (7) becomes: In other words, the Compton source remains, but the integrals in equation (7) representing electron scattering, together with the total scattering cross section, are replaced with the Fokker-Planck approximation. In equation (16), the coefficient T M (r, E) is defined as: with a similar expression for Mott scattering, except that there is no energy transfer during the elastic Mott scattering process, so the energy integral disappears (Hensel et al 2006):

Continuous slowing down approximation
The Fokker-Planck approximation neglects large-angle scattering events, which limits its accuracy. Therefore, in the continuous slowing down approximation (CSDA), energy losses are constrained to be small, but changes in direction are allowed to be larger. Consequently, in this model, the electrons are considered to scatter at discrete positions, with a well-defined change of direction at collisions, while the energy is lost continuously in between collisions. This allows larger angle scattering events to be considered, with a correspondingly small loss in accuracy. The continuous slowing down approximation simplifies equation (7) to Larsen et al (1997): where the expression, Σ E e , Ω · Ω , in the second term is defined as: and there is no corresponding expression for Mott scattering because Mott scattering does not involve energy transfer. The differential scattering cross section in the second term is now the only part to be integrated over energy and this is separated from the angular integral. S (r, E e ) is defined in equation (8) and this term provides the description of energy transfer. In other words, the second and third terms of equation (19) describe the angular behaviour of the scattering events and the fourth term describes the energy dependence. Thus, the value of the approximation is in its separating and removing integrals, which makes the problem much easier to solve computationally.
In both of the above types of approximation, there are also versions which involve a combination of true Boltzmann integrals as in equation (7), combined with these approximations. This allows large angle scattering to be carried out more thoroughly, while creating a tractable mathematical problem.

Small-angle scattering in proton and ion therapy
The LBTE can also be solved for the case of proton and ion transport with application to proton and ion therapy treatment planning. In proton transport, the elastic Coulomb interactions with the atomic nuclei and the electrons in the media mean that there are a large number of small-angle scattering events. Uilkema (2012) produced a simple LBTE model using the Fokker-Planck approximation to model these small-angle scattering events.

Phase space and discretisation
The position of a point of interest in the patient is given by three coordinates and the space of all such points is a 3D space. However, at each point in space, there are photon and electron fluences at a range of energies and in a range of directions, and each of these quantities forms an additional variable. The set of values of these variables can be thought of as a point in a space with more dimensions than just the three needed to express the spatial position, and this multi-dimensional space is commonly referred to as the phase space.
The equations given above are in their nature continuous, so that they apply for all directions and the integrals are also over all energies and angles. The method of solving the above equations for fluence and dose in practice involves discretising the equations. Specifically, this means satisfying the LBTE for certain discrete directions and energies and then discretising the integrals. The discrete directions and energies are referred to as the quadrature. This concept will be familiar to the reader from Simpson's rule for approximating a continuous integral: the integral is approximated by a series of function evaluations, weighted by appropriate weighting factors.
In the context of LBTE, however, it is reserved for the expression of fluence, as distinct from the discretisation of integrals described in equation (11).
The simplest method of discretisation is referred to as the discrete ordinates method (or the S N method for historical reasons). The scalar fluence is here approximated as a weighted sum of fluence in the discrete directions: where Ω γn and w n are the directions and weights of the quadrature, respectively, and Φ γn are the corresponding fluences in these directions. We use the photon fluence here, but the concept is also applicable to the electron equation. The value of N is referred to as the order of the quadrature. The accuracy of this method depends upon the intelligent choice of the quadrature directions and weights to suit the particle transport being considered. Alternatively, the angular fluence can be expressed in spherical harmonics, additional to the expression of scattering integrals as discrete sums as in equations (11) to (14). For simple geometries, the discrete ordinates and spherical harmonics methods are identical, but for more complex geometries, they give different results.
The distribution of fluence over space is handled by a voxelated model or a triangulated mesh, incorporating a model of fluence within each voxel. For example, the simplest model is to assume that fluence is equal across a voxel, with discontinuities at the voxel boundaries. A more sophisticated method is to assume that the fluence varies linearly across a voxel or tetrahedral finite element, again with discontinuities in both fluence and its gradient at the boundaries. This latter approach can then be solved using the linear discontinuous Galerkin method Hill 1973, Lewis andMiller 1984). The concept of this method is that the overall solution is constructed from a series of weighted approximate functions, each relating to a particular finite element. The main advantage of this method is that the solution to a particular element does not depend on the neighbouring elements. Consequently, it is compact in its implementation, and it can be applied near to boundaries without any special treatment.
If the physical reality is to be represented as accurately as possible, a large number of directions should be considered at each location in space. However, in practice, a limited number of directions are considered in order to keep computation times to a minimum. The set of directions chosen can be tailored to the energy and type of radiation being modelled. For example, if large scattering angles are expected, giving rise to radiation travelling in all directions, the directions can be chosen to be evenly distributed over the 4π space. If, on the other hand, protons or ions are being considered (see section 5.4), scattering angles are small, so the particles travel mostly forwards. In this case, the directions chosen for solution can be focussed around the direction of particle travel, with fewer directions chosen in the remaining space Adams 2002, Sanchez andMcCormick 2004). The different solid angles represented by the different directions are accounted for by weighting factors, one for each direction, with the choice of directions being referred to as the quadrature. Improved methods of discretization of space and direction have been investigated by Kópházi and Lathouwers (2015) (figure 5).
One of the disadvantages of discretised directions, is that the solutions tend to favour the basis directions, so that ray effects occur. Several works have investigated the mitigation of ray effects (Gifford et al 2006).
In terms of energy, the multigroup theory is commonly used (Lewis and Miller 1984). In this approximation, the energies are grouped into discrete groups, so that the integrals become sums. This allows various performance enhancements to be made. The energy range of the particles under consideration is divided into G energy ranges, indexed by g (1 < g < G) such that energy decreases as g increases. The boundaries of the energy ranges are denoted by E g and the upper energy bound of E 1 , generally denoted E 0 , is chosen to be sufficiently high that all particles are considered, and E G is close to zero. The angular fluence in each energy range can then be described as: If the energy distribution within each energy range, or spectral weighting function, is denoted by where f is normalised appropriately so that equation (22) holds. Using this concept, the integrals over energy in sections 2-5 become a discrete summation combined with an integral over each energy range. For example, the stopping power from equation (8) becomes: The cross sections for each energy range can be combined in the same manner as fluence in equations (22) and (23): Thus, the method provides a powerful means of discretising the energy component of the LBTE, and is consequently used universally.
In practice, the discretisation of space, orientation and energy are not entirely independent. For example, Vassiliev et al (2010) use an angular quadrature order which is adaptive according to the energy group, with a higher order for higher energy groups. The rationale for this is that particles of higher energy have longer path lengths, so that a greater range of scattering angles is produced.

Numerical framework
The practical solution of the LBTE is described in detail by Lewis and Miller (1984) and is summarised here. The transport equations have to be solved over a spatial grid, which is here taken to be Cartesian for simplicity, although as described above, a triangulated mesh is commonly more efficient. Particles are taken to be travelling according to discrete angular ordinates n, given by direction cosines µ n , η n and ξ n . As with the example above, we take the photon equation (1) as an illustration. Rewriting equation (1) in Cartesian coordinates and omitting the energy dependence, the problem to be solved is given by: where Q(r) is the generalised scattering source as defined in equation (11), i.e. the contributions from other discrete ordinates to the one being considered. The volume is divided into voxels indexed by i, j, k, bounded by x 1/2 , x 3/2 , …x I+1/2 , y 1/2 , y 3/2 , … y J+1/2 , z 1/2 , z 3/2 , …z K+1/2 , where it is assumed that the fluence is flowing from the direction of low index to the direction of high index. The cross sections are constant over each voxel. Integrating equation (26) has the effect of turning the derivatives into integrals of differences between fluence values: +ξ n´i dx´j dy Φ n x, y, z k+1/2 − Φ n x, y, z k−1/2 +ρ e (x, y, z) σ tot C,γ´i dx´j dy´k dzΦ n (x, y, z) =´i dx´j dy´k dzQ (x, y, z) where ´i dx is used to denote the integral between x i−1/2 and x i+1/2 , etc. We then define the voxel dimensions: Figure 5. Illustration of the discretisation of directions. This particular scheme is described as hierarchical sectioning of the angular sphere into patches. The surface of each octant of the sphere forms a spherical triangle, the midpoints of whose sides can be joined to create smaller triangles. A hierarchy of successively smaller triangles results from recursive application of this scheme. Reproduced from Kópházi and Lathouwers (2015), with permission.
integrals over the faces of the voxel: and integrals over the volume of the voxel: so as to simplify the following equations. If we substitute these expressions into equation (27) and divide by ∆x i ∆y j ∆z k , we have: So far, we have made no approximation, but in order to relate the fluence at the surface of the voxel to that at the centre, it is necessary to apply some sort of relationship. This is often the diamond difference relationship: and: depending upon which direction is being considered. By substituting these last three equations into equation (34), the fluences at three of the surfaces of the voxel can be eliminated, giving: The value of this formula is that it provides a relatively simple expression for the directional fluence at the centre of each voxel in terms of the inflowing fluence from the adjacent cells.

Matrix solutions
Equation (38) provides a relationship between the fluence at the centre of each voxel and the adjacent surfaces, which in turn are governed by the fluence at the centre of the adjacent voxels. The essence of the problem is therefore that the fluence at one voxel depends on the fluence at other voxels. Expressing this in a general form leads to a set of simultaneous equations: which can be solved by matrix methods. However, the linear system becomes rapidly larger with the number of variables to be solved. For example, if there N × N × N voxels in three dimensions, with both photons and electrons to be considered, and D directions and E energies to be considered, then to express the matrix becomes 2DEN 3 × 2DEN 3 , which is extremely large (Boman et al 2005). However, it can be seen from equation (38) that the radiation flowing from one voxel only ever flows into an adjacent voxel, usually a downstream voxel, i.e. with forward scatter and a loss of energy. Consequently, the matrix is sparse, and sparse matrix solution methods can therefore be used to make the solution tractable. For example, Hensel et al (2006) generate a tridiagonal system and solve it with the Thomas method. As computer power increases, this approach may eventually become relevant, but at present it is not feasible for solution of practical problems involving complex scatter conditions and tissue heterogeneities.
There are many iterative matrix methods in the literature for the solution of these linear systems, depending upon the exact form of the equations leading up to them. Most of the methods use some form of the Gauss-Seidel iteration scheme (see Press et al (1989)), in which matrix A of equation (39) is decomposed into the sum of a lower triangular matrix, a diagonal matrix and an upper triangular matrix: It follows that: and the Gauss-Seidel iteration scheme operates by successively computing: where i is the number of the iteration. Various acceleration schemes are available to improve the convergence.
In practice, many of the implementations of LBTE use a triangulated mesh for computation. For example, Gifford et al (2006) describe the Attila implementation to be based on tetrahedral elements, and the solution to equation (38) therefore involves finding 4DEN unknowns, where N is the number of such elements. In this case, the lower triangular matrix of equation (42) is composed of 4 × 4 blocks, each relating to the four unknown fluences of a tetrahedral element. Jia et al (2012) also suggested that the LBTE could be solved efficiently using graphics processing unit (GPU) technology for maximum computational performance.

Source iteration
In practice, the set of equations (38) is solved iteratively by a process known as source iteration. This involves a series of 'sweeps', in which each voxel of the grid is solved individually, going in the direction of the radiation transport, and from regions of high particle energy to regions of low particle energy. The solution for one voxel is used to provide a starting point for the solution of the next voxel. For example, if the incoming fluence at the surface of the calculation volume is known and some starting distribution of fluence is estimated, equation (38) can be used to calculate the dose at the centre of the most superficial voxel. Equations (35)-(37) can then be used to find the dose at the deepest end of that voxel. Equation (38) can then be applied to find the fluence at the centre of the next most superficial voxel, and so on.
A simple physical interpretation of this method is that the first complete sweep through the phase space corresponds to determining the contributions of unscattered particles, the second sweep corresponds to determining the contributions of first-scattered particles, the third sweep the contributions of second-scattered particles, and so on. Mathematically, the source iteration method is equivalent to expressing the entire set of equations as a matrix and then solving this large matrix by Gauss-Seidel iteration. The exact order of the sweeps is important for the convergence of the method, and it is important to sweep in the direction of the particle transport. Consequently, it is necessary to carry out a sweep for each of the quadrants in which the directional ordinates lie. Various types of acceleration are available for improvement of convergence (Lewis andMiller 1984, Press et al 1989).

A 1D example
Let us take the simplest possible case of a broad volume, whose upper edge is irradiated uniformly with a unit fluence of monoenergetic photons, and let us neglect electron transport. Also suppose that the photons only travel in the downwards direction, so that there is only one dimension of photon transport, i.e. the direction cosine µ n in section 7.1 is unity and the other direction cosines η n and ξ n are zero. A line of unit-sized voxels down the centre of the volume allows for calculation of fluence and dose (figure 6). These assumptions mean that the phase space of particles is unidimensional, which makes its solution easier to calculate and visualise. We would like to solve the LBTE for photon fluence.
The fluence, Φ 1/2 , across the upper edge of the volume is unity and the relationship between fluence entering a voxel and fluence exiting the same voxel is given by equation (34) with η n and ξ n set to zero: where we have dropped the y-and z-dimensions, and also dropped the subscript n as the order of the quadrature is 1. All the quantities we are now dealing with are total or integral fluences passing from one voxel to the next. We define: which is the integral of scattered photons joining the positive x-direction, and: (Note that the contribution a forms part of b since b describes the photons which are lost from the beam, but some of those scatter in the forwards direction to create a.) Substituting equations (44) and (45) into (43) gives: which gives us the linear system of the form:  Note the near-diagonal form of the matrix. This can be inverted by a standard matrix inversion: The result of this simple exercise is that fluence falls off as: in which Φ 1/2 is unity. The fluence at the centre of the voxels can then be calculated from equation (35): Using the source iteration approach, we begin with an initial estimate of the fluence distribution, which can be taken to be zero. In a real implementation, an initial ray tracing would be carried out to give a realistic starting estimate of photon fluence. The fluence at the centre of the first voxel, Φ 1 is calculated from equation (38). The 'difficulty' is that we do not know the value of Q i (x) because it depends on the value of Φ (x), for which we are trying to solve. However, we have an initial estimate of zero, so the value of Q i (x) is taken to be zero at this stage. We then use equation (35) to calculate Φ 3/2 . Alternating between equations (35) and (38) enables the entire fluence distribution to be determined. As this is in the absence of the scattering term, Q i (x), the solution can be thought of as the unscattered radiation, and this forms the basis for the estimate to be used for the next iteration. Repeating the process several times, using the value of Q i (x) determined from the previous iteration, yields the solution. For values of a = 0.1 and b = 0.5, the problem above takes five iterations to give a result the same as equations (49) and (50) to within four significant figures.

A 2D example
The above example is now extended to two dimensions, four directions and two energies. Thus, at each pixel, photons are considered to be travelling either up, down, left or right, which simplifies the form of equation (34). The photons have an energy of one unit or two. For simplicity, the photons are assumed to scatter either directly forwards, 90° left or 90° right, and are assumed to lose energy in the process, from two units to one. Equation (46) in this case has the form: where the superscripts 1+, 1−, 2+ and 2− are used to denote fluences with one or two units of energy, directed in the positive or negative directions. This is described in figure 7. Equation (51) is the equation describing fluence with one unit of energy, directed in the positive x-direction, but there are corresponding equations describing fluence with two units of energy, and for the other three directions of transport. The equations for two units of energy do not have any of the three a/2 terms since interactions only ever cause a loss in energy: The resulting linear system from these equations is as follows: (53) where a/2 is denoted simply by a, and b/2 by b to keep the form compact. The first eight rows of the matrix are boundary conditions, notably that the Φ 2+ i=1/2 = 1, and all of the other inflows are zero. Although equation (53) actually only solves the LBTE for one voxel, it gives some insight into the form of the matrix, with couplings between the energies and directions.
Using source iteration, the method of solving this problem is to begin with the fluence of highest energy, and sweep from the boundary, where the value is known, in the direction of transport. Once the distribution of Figure 7. Particle interactions assumed for the 2D example. Three high-energy photon streams can contribute low-energy photons to the direction of particle transport directed vertically downwards. The same scenario applies to the other directions of transport. high-energy fluence is known, the scattering sources for the low-energy fluence are defined. The sweep of lowenergy fluence can then take place, once for each of the four directions. As this problem only involves photon fluence, further iterations are not needed. Figure 8 gives an example of the results for a 10 × 10 grid of pixels for a = 0.1 and b = 0.5. The high-energy fluence component falls off with increasing pixel number, while the lowenergy component in the same direction increases in value then decreases due to the high-energy source. The low-energy component in the lateral direction increases in value across the grid.

Results of research implementations
There have been several practical implementations of the LBTE in the radiation therapy context and accuracy studies have been conducted for all of these. One of the earliest implementations of LBTE solution was the Attila solver (Wareing et al 2001), which was subsequently applied to external beam radiotherapy (Gifford et al 2006, Vassiliev et al 2008. Gifford et al (2006) evaluated this solver for simple heterogeneous media by comparing with Monte Carlo simulations and found an agreement of better than 2%. Realistic clinical cases, including prostate and head-and-neck were then considered by Vassiliev et al (2008) and the solver was found to compare with the EGSnrc Monte Carlo code typically to within 3% and 3 mm.
The Attila solver was then rewritten, giving roughly an order of magnitude increase in speed, and the resulting algorithm was described by Vassiliev et al (2010). In comparisons with the EGSnrc Monte Carlo code, the solver was found to agree to within 2% of local dose or 1 mm for a heterogeneous slab geometry. For a breast case, the agreement was 2% of local dose or 2 mm (Vassiliev et al 2010).
Other studies using different algorithms found similar promising results. For example, Boman et al (2005) and Boman (2007) investigated the performance of their algorithm for both the forward and inverse problems in very simple demonstrative geometries and found reasonable agreement with Monte Carlo simulations. Similarly, Hensel et al (2006) compared their calculation with tabulated Cobalt-60 depth dose curves and found agreement to within several percent. Simple heterogeneous geometries were also investigated. Akpochafor et al (2014) used an implementation of LBTE to model photon beam data, with application to quality assurance during treatment planning system commissioning. Several other works focussed on improving the accuracy of LBTE solvers (Olbrant et al 2009, Jörres 2015, and work is ongoing in this area, also outside of the main radiation oncology field.

Clinical implementations
The Acuros algorithm of Vassiliev et al (2010) was adopted and further developed by Varian for the Eclipse treatment planning system (Varian Medical Systems, Palo Alto, CA). Consequently, the most widespread implementation of dose calculation using the LBTE was that of the Acuros XB algorithm in Eclipse, and nearly all of the most recent publications related to this implementation (Ojala 2014, Park 2014, Rana 2014. Many of the studies compared Acuros XB with the earlier Analytical Anisotropic Algorithm (AAA) (Varian Medical Systems) (Van Esch et al 2006). Although the AAA algorithm used a 2D convolution in the style of a pencilbeam algorithm, the dose deposition kernels were adjusted for each location in the patient in both the depth and lateral directions according to the local tissue density. The resulting inhomogeneity correction allowed for photon attenuation along the incident direction and subsequently in the direction perpendicular to it, in the same manner as a convolution-superposition algorithm. Van Esch et al (2006) noted that the most significant approximation was the splitting up of the dose calculation into depth-dependent and lateral components, which was not expected to model oblique scattering events accurately. The use of a discrete number of angular sectors was also a limitation, expected to affect the accuracy near heterogeneous interfaces.

Absorbed dose to medium and absorbed dose to water
The implementation of Acuros XB allowed the choice of calculating dose to water or dose to medium, the latter being normal for Monte Carlo simulation methods. Referring to equation (9), fluence was calculated using the characteristics of the medium, but then the values of S (r, E e ) and ρ (r) could be chosen to be that for water or medium. The application of the latter quantities was a post-processing step, so did not require recalculation of fluence. As absorbed dose to water and absorbed dose to medium differed only by the values of S (r, E e ) and ρ (r), absorbed dose to water in medium could be defined without reference to a particular cavity. However, in practice, determination of absorbed dose to water in medium would require the definition of some sort of small water-filled cavity, and the concepts of chamber perturbation factors, as used in dosimetry protocols, would be applicable (Bouchard and Seuntjens 2013).
This effect was examined by , who found that the difference between reporting absorbed dose to water or absorbed dose to medium made a difference of several percent in clinical cases for prostate, lung and breast. The effect was patient-specific, even within treatments of the same anatomical area.
Other studies investigating the impact of absorbed dose to water versus absorbed dose to medium were those of Park et al (2016) and Mampuya et al (2016). The latter focussed on stereotactic body radiation therapy (SBRT) for lung, one of the most challenging techniques for accurate calculation of absorbed dose due to the small fields and low density of surrounding lung tissue. They found that values of dose to 2%, 50%, 95% and 98% of the internal target volume and planning target volume were up to 1% higher when calculating absorbed dose to medium, compared to calculating absorbed dose to water. Fogliata et al (2011bFogliata et al ( , 2011c evaluated the results of the Acuros XB algorithm against measurements, and also against the AAA in water for several energies, including flattening-filter-free (FFF) beams and found results within 1% of measurements for open fields. The same authors also investigated the accuracy of the algorithm in heterogeneous media (Fogliata et al 2011a) and found results mostly within 3% and 3 mm of Monte Carlo simulations (figure 9). They also found the calculated results to be within several percent of measurements for small stereotactic arc fields of side 1-3 cm (Fogliata et al 2012a). Similar studies were conducted by other groups (Hoffmann et al 2012, Bush et al 2011, Han et al 2011,Bueno et al 2017, Zavan et al 2018.

Comparisons of simple fields
The Acuros XB algorithm was compared against Monte Carlo simulation and a comprehensive set of measurements for 6 MV and 18 MV beams by Alhakeem et al (2015). The measurements were taken using radiochromic film and metal oxide silicon field effect transistors (MOSFETs), allowing precise measurements in the interfaces of various materials. In the region of electronic equilibrium, all methods agreed within 3%, with the exception of AAA. In the interface regions, the film and MOSFET results agreed more closely with each other than with Monte Carlo simulation, due to the finite size of the dose calculation grid with the latter. However, the measurements in interface regions were generally within around 3% of the Monte Carlo simulation results, with rather larger differences of up to 13% being observed between Acuros XB and Monte Carlo simulation. Hirata et al (2015) also compared Acuros XB, AAA and pencil beam calculations against measurements for 4 MV beams in heterogeneous phantoms. Other studies of single fields in heterogeneous media are reviewed by Kan et al (2013c).

Verification of complete treatment plans
As well as evaluating individual fields, Fogliata et al (2012b) proceeded to show the impact of Acuros XB for lung cancer treatments using conformal beams, intensity-modulated radiation therapy (IMRT) beams and volumetric modulated arc therapy (VMAT) beams created using RapidArc, by comparing with AAA. The handling of heterogeneities with Acuros XB appeared to be more accurate than AAA, and calculation time was slower for conformal and intensity-modulated beams, but a factor of 4 faster for RapidArc beams (Fogliata et al 2012b). This latter effect was due to the need in Acuros XB to only solve the directional part of the problem once, this then being used for all control points.
Complete plans were evaluated by Han et al (2012) for head and neck. They used the Radiological Physics Center's head and neck audit phantom to compare Acuros XB against measurements and found satisfactory results. Again, a comparison with AAA was also simultaneously made. The same group also used a thorax audit phantom for a comparison of Acuros XB against measurements in lung, as described in figure 10 (Han et al 2013). Measurements for head and neck treatment plans using RapidArc were also made by Kan et al (2013b), who found that Acuros XB was slightly more accurate than AAA in the regions adjacent to inhomogeneities in the nasopharynx region. These authors also found that the mean dose to the planning target volume (PTV) was lower with Acuros XB than with AAA by about 1% (Kan et al 2013a). Other authors studied pelvic (Lloyd and Ansbacher 2013, Ojala et al 2014, Koo et al 2015, breast (Guebert et al 2018) and thoracic (Kroon et al 2013, Soh et al 2018, Padmanaban et al 2014 applications. As with Monte Carlo simulation, the need for correct assignment of tissue types was recently observed (Fogliata et al 2018).
The good agreement of Acuros XB with Monte Carlo simulation was noted by the recent study of Hoffmann et al (2018), in which a variety of complex treatment plans were recalculated using both Acuros XB and a Monte Carlo code SciMoCa (Scientific RT GmbH, Munich, Germany). The agreement was almost entirely within 2% and 2 mm throughout the patients. This served to illustrate clearly that the deterministic and stochastic methods were converging very closely in their accuracy and performance.

Stereotactic body radiotherapy
SBRT has been gaining in importance in recent years and use of an accurate dose calculation is important due to the small PTV with steep dose gradients. For the case of lung SBRT, the PTV margin is frequently composed of low-density lung material, so that differences in dose occur between the PTV and the gross tumour volume (GTV) or internal target volume (ITV) (Lacornerie et al 2014, Lebredonchel et al 2017, so that alternative methods of prescription have been sought, avoiding the PTV (Bedford et al 2018). The accurate calculation of dose in this situation is therefore particularly important. Accordingly, Liu et al (2014) recalculated a large number of lung SBRT plans with Acuros XB and found differences in dose distributions compared to AAA, resulting in differences in conformity and heterogeneity. Khan et al (2014) found that for the same prescribed dose, the number of monitor units required was 2% higher when using Acuros XB compared to AAA. This effect was also observed by Liang et al (2016). The reduction in PTV dose with Acuros XB compared to AAA was found by Huang et al (2015) for flattening filter-free beams to be mainly related to the low-density part of the PTV overlapping with lung. Other studies found better agreement with phantom measurements when using Acuros XB than when using AAA (Kan et al 2012, Stathakis et al 2012, Tsuruta et al 2014, Muralidhar et al 2015. Differences between AAA, a convolution-superposition algorithm (Pinnacle 3 , Philips Radiation Oncology Systems, Madison, WI) and Acuros XB were also noted by Zhen et al (2015) for SBRT treatment of thoracic spine.

Emerging conclusions
The previous sections show that there is a large collection of evidence supporting the accuracy of the LBTE solvers. The general agreement is that this method is slightly more accurate than the AAA algorithm. A number of practical factors affect the performance of a dose calculation method, such as determination of tissue density and choice of stopping powers. The method used for handling these factors can easily make a difference of 1%-2% in the results of a dose calculation, regardless of the inherent accuracy of the algorithm used, and this should be borne in mind when interpreting these results. Nevertheless, the outcome is positive, and those studies which have compared LBTE solvers with measurements have also reported good results. Comparison with Monte Carlo simulation also shows that current LBTE solvers are comparable in accuracy to the stochastic approach.

Treatment in a magnetic field
With rapidly increasing use of MR-guided radiotherapy, either using a Cobalt-60 source (Mutic and Dempsey 2014) (2017)), it is important to validate a modern dose calculation algorithm for prediction of doses in a magnetic field. The presence of an external force means in principle that charged particles are acted upon by the Lorentz force, so that they leave one particular direction and gradually adopt another. This affects the coupling of the different directions in the LBTE, resulting in an additional term in the left-hand side of the continuous slowing down approximation (St. Aubin et al 2015): (54) where ∇ v is the gradient of the velocity, v. The derivation of the additional term is not straightforward, but can be achieved by using concepts of statistical mechanics, in which a system of particles with a probability distribution of position and momentum is acted upon by the magnetic field, so that the distribution changes (Vassiliev 2017). The additional term has a similar form to the first term, the streaming operator, and can be considered as a streaming operator in momentum space. The acceleration, a, is dependent upon the Lorentz force: Where q is the electronic charge, γ is the relativistic factor and m 0 is the electron rest mass. The magnetic field strength is given by B. It follows that the force term in (54) is given by: Figure 10. Dose-volume histograms calculated using AAA, Acuros XB with dose to water (AXB_Dw,m), and Acuros XB with dose to medium (XB_Dm,m). Reproduced from Han et al (2013), with permission.
The gradient on the right-hand side of this equation is with respect to velocity: where x, ŷ and and ẑ are unit vectors in the orthogonal Cartesian directions. This can be converted into a gradient with respect to position by using the chain rule:  where v, µ and ϕ are the spherical coordinates of the velocity. This allows the expression of the magnetic field term as: where the relativistic particle momentum is given by:  (54) by using the multigroup method for energy discretisation and discrete ordinates for angular discretisation. The angular fluence in the scattering source Q (r, E, Ω) was expanded using spherical harmonics and the macroscopic scattering cross section expanded in Legendre polynomials as described in section 5.1. The continuous slowing down term of equation (54) was handled by using a diamond difference approximation (see section 7.1). Furthermore, by expanding the magnetic field term in spherical harmonics, the explicit dependence on the angular derivatives in µ and ϕ could be calculated analytically, leading to a soluble system of equations. A discontinuous finite element method was used for spatial discretisation and the whole LBTE in the presence of magnetic fields was solved by source iteration. The results of this study showed that the LBTE could be solved in the presence of a magnetic field to give dose distributions which were within 2% and 2 mm of doses calculated by EGSnrc for a variety of field orientations with respect to the beam direction (figure 11).
Although this study was very successful, it also highlighted that if the magnetic field strength was high in comparison to the total scattering cross section, the source iteration could diverge. This was because the magn- etic field term was treated as a scattering source, and the total of the effective scattering sources became larger than the total scattering cross section. This was likely to occur for high-strength magnetic fields in lung tissue. Consequently, a further work treated the magnetic field term as a streaming term rather than as a scattering term, i.e. as the first term in equation (54) rather than the third (St. Aubin et al 2016). This new formalism was solved by using a discontinuous finite element method for discretisation of the space-angle domain. The results agreed within 2% and 2 mm of doses calculated with GEANT4 Monte Carlo simulations, even for strong magnetic fields and low densities. Further work by Zelyak et al (2018aZelyak et al ( , 2018b confirmed this behaviour and further showed that parallelism in the handling of the angular component of fluence was possible to improve computation speed. These studies also established that the correct order of sweeping in source iteration was necessary. For stability of solution, it was necessary that the sweeps should proceed in the direction of particle transport, but when the magnetic field was added, this became more complex. To overcome this, Yang et al (2018) used a finite element method for angular discretisation, with the finite elements respecting the octants of the angular space. Combined with Cartesian voxels for spatial discretisation, this yielded stable solutions which agreed with GEANT4 simulations to within 2% and 2 mm for a 1.5 T magnetic field.

Conclusions
Application of LBTE solvers to radiation therapy is a much younger field than that of Monte Carlo simulation, which has been the subject of research for many years. However, progress so far with LBTE solvers has been very promising, with one of the major treatment planning vendors offering a highly successful LBTE algorithm for photon calculation. Recent developments in the use of LBTE for calculation of absorbed dose in a magnetic field shows that the area of research is continuing to develop rapidly, and it is expected that LBTE will continue to form an important and significant role in future treatment planning. The iterative nature of the solvers suggests that they may find application in adaptive radiation therapy, where a previous solution may be iteratively refined to suit a new patient geometry without the requirement to fully recalculate dose.