Continuous energy adjoint transport for photons in PHITS

Adjoint Monte Carlo can be an efficient algorithm for solving photon transport problems where the size of the tally is relatively small compared to the source. Such problems are typical in environmental radioactivity calculations, where natural or fallout radionuclides spread over a large area contribute to the air dose rate at a particular location. Moreover photon transport with continuous energy representation is vital for accurately calculating radiation protection quantities. Here we describe the incorporation of an adjoint Monte Carlo capability for continuous energy photon transport into the Particle and Heavy Ion Transport code System (PHITS). An adjoint cross section library for photon interactions was developed based on the JENDL4.0 library, by adding cross sections for adjoint incoherent scattering and pair production. PHITS reads in the library and implements the adjoint transport algorithm by Hoogenboom. Adjoint pseudo-photons are spawned within the forward tally volume and transported through space. Currently pseudo-photons can undergo coherent and incoherent scattering within the PHITS adjoint function. Photoelectric absorption is treated implicitly. The calculation result is recovered from the pseudo-photon flux calculated over the true source volume. A new adjoint tally function facilitates this conversion. This paper gives an overview of the new function and discusses potential future developments.


Introduction
The Monte Carlo method can be inefficient for solving radiation transport problems with a small tally relative to the source.The inefficiency is that many particle histories during the course of a simulation fail to cross the tally region, and hence make no contribution to the calculation result.Large source, small tally geometries can require long computing times to converge the calculation result to a desired level of precision.
Adjoint Monte Carlo is one strategy for improving the efficiency of large source, small tally transport problems.In adjoint Monte Carlo, radiations are generated within the tally region and transported backwards through the geometry, undergoing any physical interactions in reverse.They contribute to the calculation result if they cross the source region.The efficiency benefit is conferred by solving a problem where the 'adjoint source region' is relatively much larger than the 'adjoint tally region', and thus a larger fraction of particle histories score [1].
The adjoint method is a promising strategy for environmental radioactivity transport problems.In these problems it is often necessary to calculate gamma ray dose rates at precise locations, where radioactivity dispersed widely within the environment contributes significantly to the result.For example, the source may be natural potassium, uranium and thorium series radioisotopes within soil, or e-mail: malins.alex@jaea.go.jp radioactive fallout from a nuclear accident or test.In such cases decay photons emitted hundreds of metres from the detector are important [2].
Hoogenboom previously outlined a method of adjoint Monte Carlo for photons with continuous energy representation [3,4].One benefit of continuous energy representation over a multigroup scheme is that energy spectra can be evaluated with fine resolution, allowing the simulation of gamma spectrometer responses.Continuous energy representation is also important for accurately calculating radiation protection quantities, such as ambient dose equivalent (H * (10)).However, there have only been a few implementations of adjoint Monte Carlo for photons with continuous energy representation in popular transport codes [5,6].This paper outlines the addition of a continuous energy photon adjoint function to the Particle and Heavy Ion Transport code System (PHITS) [7].There are three main components to the code developments.First, an adjoint cross section library was developed adding adjoint incoherent scattering and pair production cross sections.Second, a new adjoint function was added to PHITS, allowing the transport and collision of pseudo-photons by coherent and incoherent scattering, and implicit treatment for photoelectric absorption.Finally, an adjoint tally function was developed to recover the forward calculation result from the pseudo-photon flux tallied over the real source volume.

Overview of forward Monte Carlo photon transport
Monte Carlo radiation transport codes sample the Boltzmann transport equation The subject of the equation is the emission density of photons with energy E [eV], in the direction Ω, from point r in space, i.e. χ(r, E, Ω) [cm −3 eV −1 sr −1 ].There are two contributions to χ(r, E, Ω): photons emitted by radiation sources, S (r, E, Ω) [cm −3 eV −1 sr −1 ], and photons emerging from collisions, given by Photons emerging from collisions at r first have to be transported to the collision site.This process is described by the transport operator, T , which operates on the emission density χ(r , E , Ω ).The result of the operation is the collision density ψ(r, E , Ω ) [cm −3 eV −1 sr −1 ]: The integral here is over all path lengths L = |r − r | [cm] by which a photon can travel from r to r.The factor Σ t (r, E ) [cm −1 ] is the total macroscopic cross section for a photon with energy E to undergo an interaction at r. Assuming there is no directional dependence of the cross sections within the transport media, and that polarization effects can be discounted, the interaction cross sections are independent of a photon's direction of travel.The exponential term in Eq. ( 2) describes the attenuation faced by photons traversing towards the collision site, where The outcome of collisions is described by the collision operator The collision operator has intentionally been expended here into various factors, corresponding with steps or decisions taken in Monte Carlo sampling of the Boltzmann equation.The integral is over all possible incident photon energies and incident directions of travel to the collision site (note dΩ is a solid angle).The first factor is the probability that a scattering rather than an absorption interaction occurs: Here Σ s (r, E ) [cm −1 ] is the total macroscopic cross section for all scattering interactions.Monte Carlo codes with implicit treatment of the photoelectric effect, i.e. photoelectric interactions are modelled as pure absorption processes into E = 0 eV, adjust the photon weight by P s (r, E ) upon each collision.Note PHITS offers more advanced treatments of the photoelectric effect in forward mode, including simulation of luminescence photons and ejected electrons [8], but the implicit treatment is outlined here as that is what is used in the adjoint mode.
The second factor in Eq. ( 4) is the probability that the scattering interaction occurs with an atom of type A amongst the various nuclides in the material: The sum of numerator over all atom types in the material is the denominator, i.e.Σ s (r, E ) The third factor is the probability for a type j scattering interaction to occur upon collision with the type A atom: Again the sum of the various possible numerators in this equation is the denominator, Σ s,A (r, E ) The macroscopic interaction cross sections are related to microscopic interaction cross sections, σ j,A (E ) [b], by Here n A (r) [cm −3 ] is the number density of atoms of type A in the material at r.The factor of 10 −24 ensures Σ j,A (r, E ) has units of cm −1 .The final factor of the collision operator (Eq.( 4)) is the normalized probability density function (pdf) for the possible outcomes of the scattering interaction: The first factor on the right hand side of this equation is the differential microscopic cross section.Monte Carlo codes sample Eq. ( 1) by simulating photon histories.A history is spawned by sampling the source term S .The photon is transported through the geometry to a collision site by sampling T .The weight of the photon history is multiplied by the factor P s prior to sampling the collision itself.The collision nuclide and the type of interaction are selected according to P A and P j,A , respectively.The outcome of the interaction is determined by sampling f j,A .The process of transport followed by collision repeats until either the photon leaves the simulation geometry, or its energy or weight fall outside the bounds of interest, upon which the history is terminated.
The result of a Monte Carlo simulation is typically calculated using a track length tally for the photon flux combined with a chosen response function.This corresponds to solving the integral: Here η φ (r, E , Ω ) is the response function of choice, and φ(r, E , Ω ) [cm −2 eV −1 sr −1 ] is the photon flux, related to the collision density by

Adjoint Monte Carlo formalism
The PHITS adjoint function follows the formalism outlined by Hoogenboom [3,4] and Gabler et al. [5].The adjoint Boltzmann transport equation is: where superscript pluses and diacritic tildes denote adjoint quantities and operators, respectively.The first quantity, ζ + (r, E, Ω), is the adjoint analogue of the forward emission density, i.e. the density of pseudo-photons emitted at r.The operators T and C are analogues of the forward transport and collision operators, respectively.Eq. ( 12) is sampled with the Monte Carlo method by simulating pseudo-photon histories.A pseudo-photon history is spawned by sampling an initial position, energy and direction from the 'adjoint source term', η φ .This is in fact the tally response function for the corresponding forward transport problem.For a typical response function, such as a flux calculation, a pseudo-photon is spawned at a random location within the forward tally volume.The initial energy is chosen randomly from within the energy range of interest at the forward tally, and the initial direction of travel is also random.
The adjoint transport operator is basically the same as the forward transport operator: The result of the transport operation is the pseudo-photon collision density, ξ + (r, E , Ω ).A pseudo-photon is transported to a collision site by randomly selecting a path length through the transport media from the exponential distribution in Eq. ( 13).
Although similar in form, the adjoint collision operator is different from the forward collision operator in a number of respects: Each of the factors in this integral are described in turn.
The first factor, P + s (r, E ), is the adjoint analogue of the forward scattering probability (P s ): The denominator here is the same as the forward case, i.e. the total macroscopic cross section for both absorption and scattering (Eq.( 5)).However the numerator, , is a new quantity.It is the total adjoint macroscopic cross section for the scattering interactions.Note P + s is not a probability per se, as it can take values greater than 1.The weight of a pseudo-photon is multiplied by P + s (r, E ) upon each collision.The second factor in Eq. ( 14) is the probability that the adjoint scattering interaction occurs with an atom of type A in the medium: This is a true probability as the sum of the possible numerators totals the denominator: E ).The third factor is the probability for interaction type j to occur amongst the possible scattering interactions: Again the sum of the possible numerators, i.e. partial adjoint macroscopic cross sections for scattering, is the denominator: As for the forward case, the adjoint macroscopic cross sections are related to adjoint microscopic cross sections by the number density of atoms in the collision material: The adjoint microscopic cross sections, σ + j,A (E ), also have units of barns.They are calculated by integration of the differential cross section for the corresponding forward interaction.The integration range covers all incident photon energies that could yield an outgoing photon with energy equal to the colliding pseudo-photon: The final factor in Eq. ( 14) is the probability density function for the possible outcomes of adjoint scattering: The purpose of the adjoint simulation is to calculate the tally response of the corresponding transport forward problem.The response is given by: R = . . .S (r, E , Ω )

Interactions
The PHITS adjoint mode currently supports coherent and incoherent scattering, and implicit photoelectric absorption.This section outlines calculation of the adjoint microscopic cross sections, and the Monte Carlo sampling strategies for the outcome of adjoint collisions.

Forward
Coherent scattering between a photon and an atom leads to a change in the direction of travel of the photon, but, to very good approximation, the energy of the photon is preserved.The differential cross section for the interaction is [9]: where r e = 2.817 × 10 −15 m is the classical electron radius, and μ = Ω • Ω is the cosine of the scattering angle.Eq. ( 22) consists of two parts: Thomson's cross section for scattering off a free electron, and the form factor correction to this, F(Z, x ), which accounts for atomic binding effects.Z is the atomic number of the collision atom, and is the inverse wavelength of the recoiling atom: Here e = 1.602 × 10 −19 J eV −1 , h = 6.626 × 10 −24 J s is Planck's constant, and c = 2.998 × 10 8 m s −1 is the speed of light.
In a PHITS forward mode calculation, coherent scattering interactions are simulated as follows.The photon cross section library (based on JENDL-4.0 [10]) contains tabulated point data for the microscopic cross sections, σ coh (E ), calculated by integrating Eq. ( 22) numerically.The cross sections are stored against logarithmic values of the incoming energy E .The library also contains tables for the form factor at point values of x , and for the integral of the form factor at point values of x 2 .
Upon a coherent scattering interaction, the scattering angle cosine is sampled from f coh (E , μ) using the rejection technique method from Carter and Cashwell [11].A scattering azimuth is selected randomly from [0, 2π].

Adjoint
The adjoint microscopic cross section for pseudo-photon coherent scattering from (E , Ω ) to (E, Ω) is: The steps above use the facts that E = E, μ = μ and x = x for this interaction.The adjoint microscopic cross sections are thus identical to the forward microscopic cross sections at all energies, i.e. σ + coh (E ) = σ coh (E ).Likewise the probability density function for the scattering cosine is the same as the forward case: PHITS thus samples a pseudo-photon coherent scattering event identically to a forward photon coherent scattering event.The code employs the same cross sections and rejection technique for the pdf as per the forward case.

Forward
Incoherent scattering occurs between a photon and a quasifree atomic electron.There is a correlated change in the energy and the direction of propagation of the photon after the interaction.The outgoing photon energy is: where k = E /E e is the incident photon energy relative to the rest energy of an electron, E e = 0.511 MeV.The differential cross section is [9]: The factor S (Z, x ) is the incoherent scattering function, which corrects the Klein-Nishina cross section for atomic binding effects.The second argument is the wavelength of the recoiling target: For a forward incoherent scattering event, the PHITS cross section library contains tables for σ incoh (E ) against ln(E ), obtained by numerical integration of Eq. ( 27).The library also contains S (Z, x ) values against x .
To sample the pdf for the scattering cosine, a candidate is drawn in accordance to the Klein-Nishina component of Eq. ( 27) first.Kahn's rejection method is employed for photons below 1.5 MeV [12] and Koblinger's method for E > 1.5 MeV [13], following Blomquist and Gelbard [14].The candidate is then accepted or rejected against the incoherent scattering function [11].Note for the calculation of S (Z, x ), Monte Carlo codes (e.g.MCNP [15]) commonly approximate x using Eq. ( 23) instead of Eq. (28).

Adjoint
In adjoint mode, pseudo-photons gain (or preserve) energy upon incoherent scattering.The outgoing pseudo-photon energy is: The adjoint microscopic cross section is given by integrating numerically: During numerical integration, the scattered pseudo-photon energy E must be calculated using μ via Eq.( 29) in order to evaluate x.It is necessary to set a minimum scattering cosine for when , where E max is the maximum energy of interest in adjoint photon calculations: (31) We chose E max = 3.0 MeV for the PHITS adjoint mode, as this range covers the main radioactive decay photons.Without setting E max , the integral in Eq. (30) would diverge for E ≥ E e /2.
An example of the adjoint microscopic incoherent scattering cross section is given in Fig. 1 for scattering off an oxygen atom.The adjoint cross section increases from low energies to peak at E max /(1 + 2E max /E e ).At energies above this, the cross section decreases because the pseudo-photon scattering pathways to energies above E max are truncated.
The shape of the adjoint cross section curve and the position of the peak in Fig. 1 necessarily depend on the value chosen for E max .As an example, the curves for germanium under two different choices of E max are given in ref. [5].
The probability density function for the outcome of pseudo-photon incoherent scattering is: The 'adjoint Klein-Nishina component' of this pdf can be sampled using the following rejection method.First a candidate μ is selected uniformly in [μ min , 1].This is accepted if: where r is a uniform random deviate in [0, 1] and m is the maximum value of the numerator in [μ min , 1]: otherwise.
(34) Note the formula for m is conditional on the incoming pseudo-photon energy E .This is because the maximum of the adjoint Klein-Nishina distribution switches from occurring at μ = μ min to μ = 1 around a crossover energy at E = 2.020 MeV (for E max = 3.0 MeV).The crossover energy can be established for E max 3.0 MeV by numerically solving The acceptance efficiency of the adjoint Klein-Nishina distribution rejection routine per candidate scattering cosine is: The variation of e with the incident pseudo-photon energy E is shown in Fig. 2. In forward mode, Kahn's method has a fairly flat acceptance efficiency of ≈ 60 % between 1 keV and 1.5 MeV.The rejection method developed for the adjoint mode generally has good acceptance efficiency, however there is a dip in efficiency between 150 and 650 keV.Each candidate μ sampled from the adjoint Klein-Nishina distribution must then be accepted or rejected in line with the incoherent scattering function.Again μ must be used to evaluate x (approximated with Eq. ( 23), as per the forward case) in order to rejection sample S (Z, x).

Photoelectric effect
PHITS adjoint mode implements an implicit treatment of the photoelectric effect.This marks a departure from forward mode, which offers more detailed physics simulation for the relaxations associated with photoelectric interactions.
In adjoint mode, the weight of a pseudo-photon is multiplied by P + s upon each collision, then either a coherent or incoherent scattering event is forced.As it is possible for P + s be greater than one, pseudo-photon weights can become large and cause high variances in calculation results [5].This risk should be monitored in adjoint calculation results.

Adjoint Cross Section Library
The original PHITS photon cross section library, which is based on JENDL-4.0 [10], was extended to facilitate the adjoint transport function.The adjoint specific library files are encoded in the ACE format.There is one file for each element, with filenames XY000.j90p.Characters X and Y are substituted with the element symbol, e.g.sodium is Na000.j90p.The second character is substituted with an underscore for elements with single character symbols, e.g.O_000.j90p for oxygen.
The first nine blocks of data in the new adjoint library files are identical to the original PHITS photon cross section files (Table 1).The blocks contain energies and microscopic cross sections for forward incoherent, coherent, photoelectric and pair production interactions, all stored logarithmically.They are followed by blocks containing incoherent scattering functions, integrated coherent form factors, coherent form factors, and average heating numbers.
Two new blocks follow to facilitate adjoint incoherent scattering.Block 10 contains the logarithms of energies as per block 1, but with an additional data point added at E max /(1 + 2E max /E e ).The maximum in the microscopic adjoint incoherent scattering cross section occurs at this energy (Fig. 1).Block 11 contains logarithms of the microscopic cross sections for adjoint incoherent scattering, obtained by numerical integration of Eq. (30).PHITS interpolates between the adjoint incoherent scattering cross sections using log-log interpolation, as per the forward case.
Note separate data blocks are unnecessary for adjoint coherent scattering or photoelectric effect cross sections.The adjoint coherent scattering cross sections are equal to the forward cross section (section 4.1.2).As the adjoint photoelectric treatment is implicit, adjoint cross sections are not needed for this interaction.However, the forward photoelectric cross sections need to be read in to calculate Σ t (r, E ) for both adjoint transport and collision operations.

Pseudo-photon transport in PHITS
An adjoint simulation is called by the iadjoint = 1 flag in the [Parameters] section of a PHITS input file.This instructs the code to simulate the pseudo-photon particle type.The user defines a source region in the input file that coincides geometrically to the tally region in the corresponding forward transport problem.Pseudo-photons are generated within the source with energies uniformly Source Tally 5 cm 5 cm 5 cm between the minimum energy of interest for the forward tally and the highest photon energy emitted by the forward radiation source.Each pseudo-photon is tracked through the geometry to a collision site, upon which its weight is adjusted by P + s to account for non-absorption by the photoelectric effect.Then either a coherent or incoherent scattering interaction is forced.Coherent scattering is sampled identically to a forward mode calculation (section 4.1.1).A scattering cosine is sampled first from the Thomson distribution and then rejected against the coherent scattering form factors.A scattering azimuth is selected uniformly in [0, 2π].
Upon incoherent scattering, the pseudo-photon scattering cosine is sampled according to the adjoint Klein-Nishina distribution and rejected against the incoherent scattering function (section 4.2.2).A scattering azimuth is selected uniformly in [0, 2π].

PHITS adjoint tally
To obtain the calculation result the user defines the adjoint tally type in the input file ([T-Adjoint]).This tally type is the parallel of the forward track length tally, PHITS tallies a pseudo-photon if it crosses the adjoint tally region with energy within the photon source energy range.The tally calculates the contribution of a pseudophoton history to the calculation result using a track length estimator: WT l /V.Here W is the pseudo-photon weight, T l [cm] is the track length through the tally region, and V [cm 3 ] is the tally volume.Importantly, in a departure from the usual forward case, it is the energy of the pseudophoton history upon generation in the adjoint source region that determines the energy bin that the history is tallied within for the [T-Adjoint] tally.This process recovers the calculation result of the corresponding forward problem (i.e.Eq. ( 21) is solved by [T-Adjoint]).

Examples
We present two example calculations to demonstrate the PHITS adjoint function.The first is a calculation for a !" photon spectrum demonstrating the accuracy of the adjoint method.The second is an application to a simple large source, small tally transport problem, to gauge the potential performance of the adjoint method for this class of problem.

Spectrum calculation
This calculation was for gamma ray spectrum inside an infinite oxygen medium, density ρ = 1.0 g cm −3 .The source and tally volumes were adjacent cubes, volume 125 cm 3 (Fig. 3).Photons were emitted from the source with uniform probability between 0.95 and 1.0 MeV.The photon flux was tallied in 50 keV wide energy bins.
The calculated spectra are shown in Fig. 4. The adjoint method gives a basically identical result to the forward method at higher energies, to within statistical error.There is however a discernible difference between the fluxes calculated in forward and adjoint mode for the lowest energy bins (0.05-0.1 MeV).This is an artifact of the implicit treatment of the photoelectric effect in adjoint mode, which is an important interaction at these low energies.The forward mode calculation includes relaxation photons from photoelectric interactions, leading to the higher flux results in the lower energy bins.The second test problem consisted of a photon source that was relatively much larger than the tally (Fig. 5).The source and tally were within an infinite scattering medium consisting of oxygen, density ρ = 0.01 g cm −3 .A 1.0 m 3 source emitted photons uniformly in the energy range 0.05-1.0MeV.A 1.0 cm 3 tally was located on the centre line from one of the source faces.The separation between the source and tally was 150 cm.The calculation result was the total flux over the 0.05-1.0MeV energy range at the tally.
The calculation was performed in adjoint mode for 7.5 min run time on a 3.4 GHz four-core desktop computer.For comparison with an accelerated (biased) forward method, a calculation was executed with the PHITS point tally ([T-Point] [8]) for the same CPU run time.Finally the calculation was run in unbiased forward mode with the [T-Track] tally for 436 min, i.e. 58 times longer than the point and adjoint calculations.
The results of the fluxes and tally relative errors for the three calculation modes are given in Table 2.The relative error for the unbiased forward tally was an order of magnitude higher than for the point and adjoint calculation methods.This is despite the much longer run time and larger number of histories simulated.The adjoint and point tallies gave comparable flux results and calculation relative errors.

Conclusion and future developments
A new adjoint function has been added to PHITS for photon transport problems with continuous energy representation.Currently the function can treat coherent and incoherent scattering, and implicit photoelectric absorption.The function is based on three developments to PHITS: a new adjoint cross section library, a pseudo-photon simulation mode, and an adjoint tally function.
The new mode is a candidate for solving transport problems with small detector sizes relative to the source.In a simple test case for this class of transport problem, the adjoint mode offered orders of magnitude speed-up versus an unbiased forward calculation.The adjoint performance was comparable to a biased forward calculation using the point tally.A benefit of the adjoint method over the common point and ring tally forward biasing methods is it allows the tally volume to be any shape.
A number code developments are planned to improve the functionality of the PHITS adjoint mode.The developments are based on the point energy function [5], where pseudo-photons are forced to scatter to discrete energies.This will allow the incorporation of pair production inter-actions [4].Pseudo-photon histories will be scattered to E e to enable this interaction to occur.
The point energy function will also allow calculations for photon sources with discrete energies (e.g.line sources, such as radioactive decay gamma rays) [5].Currently such calculations are not possible in the PHITS adjoint mode, as pseudo-photon histories are unlikely to cross the adjoint tally region with the requisite discrete energies to contribute to the calculation result.With a point energy function, pseudo-photons histories will be split off and forcibly scattered to the discrete energies after each collision.These sub-histories will then contribute to the calculation result if they cross the adjoint tally region.

Figure 1 .
Figure 1.Forward and adjoint incoherent scattering cross sections for oxygen.

Figure 2 .
Figure 2. Acceptance efficiency of rejection methods for sampling the forward and adjoint Klein-Nishina distributions.

Figure 3 .
Figure 3. Geometry of the flux spectrum calculation.
[T-Track] [8].The geometry chosen for the [T-Adjoint] tally ([mesh] input) must correspond to the source geometry in the forward problem.The energy mesh of this tally ([e-type] input) should be the same as for the forward case.Finally the user defines the energy distribution of the photon source in the forward problem via input flags at the end of [T-Adjoint].

Figure 4 .
Figure 4. Photon flux spectrum calculated with forward and adjoint Monte Carlo for the geometry in Fig. 3.

Figure 5 .
Figure 5. Geometry of simple large source, small detector test problem.

Table 1 .
Data blocks in PHITS adjoint ACE cross section files.The block length N is defined in the ACE file header.Note: zero barn cross sections are stored as 0.0 in the tables.

Table 2 .
Tally results, relative errors and CPU run times for flux calculation for the large source, small detector test problem.Forward mode results were calculated with the conventional [T-Track] tally (no acceleration), and with the point tally biasing method ([T-Point]).Adjoint mode results are in the [T-Adjoint] column.