APPLICATION OF THE PHOTOMETRIC THEORY OF THE RADIANCE FIELD IN THE PROBLEMS OF ELECTRON SCATTERING

The physical model of the radiance field is similar in some aspects to the elementary particle transport theory under the assumptions of the classical mechanics. Disregarding the differences in the used nomenclatures, it can be shown that the transport equations for the radiance field are identical to those for the particle flux density. Since the end of the 19th century, both theories have been developing in parallel, thereby enriching each other. In other words, a breakthrough, which has been made in one theory, readily contributes to the significant progress in another one. Nowadays the accuracy achieved in the experiments with particles is close to the limit, which allows validating the relationships derived within the light scattering theory. Besides, the experiments with particles are free from uncertainties in the scattering medium, which are typical for atmospheric remote sensing applications. In this paper, a new algorithm is described, which is derived by analogies between these theories. It is applied for calculating the electron flux elastically scattered by plane-parallel layers of a solid with the strongly forward peaked phase functions. The calculations are compared against the experimental angular distributions of electrons, which are elastically reflected by the two-layer solid samples.


Introduction
The formulation of the transport theory for elementary particles is similar to the laws of the light beam propagation considered in the ray approximation framework.Moreover, these laws are valid for all classical particles, as long as they can be localized in space.The ray approximation can be regarded as a case of a more general quantum physical concept [1], in which a photon is considered as a small particle moving along a trajectory; the latter is referred to as a 'ray.'The density of the photon flux is associated with the radiance ( , ) L rl at the given point r in the direction l .Under this setup, the radiance field is given by the radiative transfer equation (RTE): here ε ,  are the extinction and scattering coefficients, respectively, while ˆ( , ) x  ll is the single scattering phase function.Initially, the RTE was introduced in [2] for a medium without scattering.In this case, for optically thin media, the RTE is reduced to the Bouguer law, i.e. exponential brightness attenuation along a ray.
The transfer equation for particles looks similar to Equation (1), namely, ˆˆˆˆˆˆ( , ) ( , ) ( ) ( , ) ( , ) ( , ) , 4  rl is the particle flux density at the given point r along direction l , σ el is the elastic scattering cross-section, while in  is the inelastic scattering cross-section.Equations ( 1) and ( 2) can be used to determine the flux density of photons and electrons reflected by multi-layered inhomogeneous media with the underlying surface.Essentially, they are the basis for the forward modelling in the Earth remote sensing retrieval codes and in the electron spectroscopy processing algorithms (in particular, describing the process of electron lithography, determining the backscatter factor in the X-ray spectral analysis as well as analysing the X-ray photoelectron spectroscopy (XPS) and reflected electron spectroscopy (RES) data).
The description methodologies of both electron and optical scattering are based on the transport equations.Therefore, the methods designed for solving optical problems can be used in the problems related to atomic particle scattering in solids [3][4][5].
By using the electron spectroscopy framework, the problem of quantitative composition analysis of samples can be readily solved.However, such an analysis of the experimental data typically requires a model of multiple elastic and inelastic scattering of electrons in solids.
Retrieval of the component composition of multi-layered samples is an ill-posed problem.It incorporates the analysis of the energy spectra of emitted electrons and extracting the information of the interest.In this regard, there are additional requirements and constraints on the scattering signal modelling algorithms and their robustness.
The state-of-the-art methods of the XPS and the elastic peak electron spectroscopy (EPES) analysis are based on rather simplistic models, in which the process of multiple elastic scattering is neglected.One of the representatives of such models is the straight-line approximation (SLA) [6,7].The main advantage of SLA-like models consists in their simplicity.However, SLA leads to a biased estimate of the scattered electron signal [8].Essentially, the error of the SLA approach is because the elastic scattering cross section el  in the situations, which are relevant for XPS, RES, and EPES, exceeds the inelastic scattering cross section, i.e.
el in   .Several techniques were developed in order to compensate the SLA bias.However, they are ad hoc and do not take rigorously into account all the factors that cause the methodological errors.Unlike in remote sensing applications, in the electron scattering spectroscopy, there are several independent methods for layer-by-layer component analyses of samples.In its turn, it is possible to experimentally validate techniques which describe the reflection processes from multi-layered inhomogeneous media.Besides, it is possible to control the composition of the sample under study, even at the stage of its preparation.In this regard, layer-by-layer composition retrieval based on the analysis of electron spectra can be performed by using techniques provided by the radiative transfer theory.In this paper, it will be shown that the experimental data on angular distributions of elastically reflected electrons [9][10][11] is beneficial for validation of the radiation codes and models used in the Earth remote sensing operational algorithms [12][13][14][15].The similar ideas were behind the verification study [16].The invariant imbedding method was proposed by Ambartsumian in the 40s of the last century to describe the processes of radiation transfer in the atmospheres of stars and planets [17,18].Essentially, this method converts the RTE for the radiance into the equations for the reflection and transmission coefficients of a slab.The method was further developed in the works of Chandrasekhar, Sobolev, and others [19,20], primarily, for spherical and Rayleigh single scattering phase functions [17][18][19][20].In these cases, an iterative procedure appeared to be efficient.As the first approximation, a single scattering solution was used [18].Typically, four iterations were enough to achieve the convergence.However, strongly forward peaked phase functions with dominant small-angle scattering are of great practical importance.In particular, these are the cases relevant for electron scattering in solids and photon scattering in a turbid medium.In this paper, it will be shown that in the case of dominant small-angle scattering it is possible to simplify the solution procedure since the nonlinear Chandrasekhar equations can be linearized.
The forward peaked phase functions were considered by Gaudsmith and Saunderson [21,22] to solve the electron transport equation by using the small-angle approximation and the spherical harmonics method.Let us consider an infinite medium with the sources of light (or particles) placed in the center (the medium boundary is placed at 0 z  ); the source satisfies the following condition: In [23] Scott developed a technique for solving transport equations by using the small-angle approximation.Dashen [3] applied the Ambartsumian equation, which solves the boundary problem of reflection from a semi-infinite medium, to describe the electron scattering process.Successful attempts to solve the linearized Ambartsumian and Chandrasekhar equations in the small-angle approximation were made in [4,24].Throughout the paper, the small-angle approximation is associated with the following small parameter In this paper, the small-angle approximation is used to derive analytical expressions for radiation reflection processes both from a semi-infinite layer and from layers of finite thickness.The small-angle solutions of the Chandrasekhar equation for the transmission function will be given.It will be shown that the small-angle approximation can be used only for Chandrasekhar equations (written for the transmission function [19]), but not for the alternative formulation from [25].An iterative procedure will be constructed that solves the problem of reflection from multilayer samples with an underlying surface at the bottom.
The main advantage of approximate analytical solutions consists in the high computational speed, which is a prerequisite for solving inverse problems by the fitting procedure.Note that the latter is the most robust approach to deal with ill-posed problems of mathematical physics, namely, remote sensing retrieval and quantitative composition analysis using electron spectroscopy [26,27].
The paper will present methods for the numerical solution of the Chandrasekhar equations.Currently, Monte-Carlo (MC) simulations [28,29] are mainly used to describe the energy and angular spectra of electrons.However, MC codes require much computational time: e.g., on a standard laptop, this time ranges up to several minutes; whereas the same simulations but based on the numerical solution of the Chandrasekhar equations, presented in this paper, are performed within fractions of a second.Note, that the idea of using numerical methods to interpret the spectra of electron spectroscopy is inherited from the radiative transfer theory.
The approbation of the approximate small-angle methods developed in this work is performed by using a comparison with exact numerical solutions as well as with the experimental data.

Chandrasekhar equations for reflection and transmission functions, linearization procedure, solution by using the small-angle approximation
The equation for the reflection function derived by Chandrasekhar [19] for a layer of finite thickness reads as follows: where τ is the dimensionless layer thickness, m is the azimuthal expansion index,  is the reflection function, x is the single scattering phase function, / ( ) with  0 =arccos  0 , 0  and  =arccos  ,  being the polar and azimuth angles of probing and viewing angles, respectively; the polar angles are defined with respect to the axis, which is perpendicular to the sample surface and looks towards the surface.Similar to Equation (4), the equation for the transmission function T reads as follows: ) .

Numerical solution for reflection and transmission functions
We note that the methods discussed in this section were first developed for solving the problems of radiative transfer and significantly enriched the theory of electron transfer as well.To get the matrix form of equation ( 1), the continuous dependency on  should be replaced with a discrete set of N values i   , while the integrals should be represented through the quadrature formulas.Then the reflection function m S turns into a matrix of dimension NN  , where i s are the weights of the quadrature method, i   are the quadrature nodes for the cosines of the incidence/observation angles, diag ( ) / ii ws   .Bearing that in mind, equation ( 4) takes the following form:  .The index + in the phase function is related to the process in which the propagation direction after the scattering event is preserved, while the index -indicates that the upward propagation changes into downward one and vice versa.Equation ( 7) is referred to as the differential algebraic Riccati equation [30,31].It can be solved numerically by several numerical techniques [31][32][33].In this paper, we use the BDF (Backward Differential Formula) method [34].The situation is somewhat simplified since the matrix m x  is symmetric.
Having performed a similar discretization for particles elastically scattered inside a layer, we obtain the following matrix equation for the transmission function: Let us find the solutions of the equations obtained, considering the nonlinear terms, and the solutions of the linearized matrix equations based on the BDF method.For elastic scattering we use the Henyey-Greenstein phase function, which is well-known in optics.The numbers in Figures 1 and 2 show the value of the asymmetry parameter g , which determines the degree of elongation of the Henyey-Greenstein phase function (2 1) ( ), ( 12 ) where P l are the Legendre polynomials.
The solution of equation ( 7) (and consequently, equation ( 4)) is shown in Figure 1 with solid lines, while the dashed line depicts the solution to equation ( 7) with non-linear term neglected (i.e., the second term on the right-hand side).The plots in Fig. 1 reveal the increasing of the computation error of the linearized equations.That is quite an expected result since the phase functions with 0,5 g  do not have a dominant small-angle scattering part.However, even in case of violation of condition (3), the linearized equations provide a result with an error less than 10%.
In Figure 2, the solid lines represent the solution of equation ( 8), the dashed line shows the solution of equation ( 8) with neglected terms containing   0 m  .In other words, this is a solution to equation (6), in which the last two terms on the right-hand side are neglected.

Analytical solutions of linearized equations for the reflection and transmission functions
The analysis performed on the basis of numerical solutions indicates that for strongly forward peaked phase functions (see condition (3)), the processes of reflection and transmission through the layer can be described by using the following equations: ) and The method of spherical harmonics is based on the representation of functions through the Legendre polynomial series.However, to be able to apply the method of spherical harmonics (which is a standard procedure for analytic solution of equations ( 10) and ( 11)), based on the orthogonality property of the Legendre polynomials within the domain [-1,1], it is necessary to perform the analytic extension of the integrands to the domain (-1, 1).For equation (10) this step is trivial since in the domain (0, -1) the integrand goes to zero due to condition (3); the analytical continuation for equation ( 11) is not straightforward.There are several approaches to how to deal with this problem, (see, e.g.[2]).The most efficient implementation is given in [35,36].The solution of equation (11), in which the integration limits are extended into the range [-1, 1] in the integral term can be found by using the idea of Goudsmit and Saunderson: namely, in the multiplier 1/  in the first term of the left-hand side of (11) and in the exponent in the first term of the right-hand side of (11) (with respect to 0 /  ), the values 0 ,  are considered to be constants.This approximation means replacing the real path by the projective one.A significant error occurs if the transport path of a photon or electron   The error in this case does not exceed 5%, as long as /1 tr el ll .Taking into account the assumptions made after the substitution of the Legendre polynomial expansions into equation (11), we obtain a system of separable differential equations; applying the boundary condition ( ) 1 lm T   , the solution for the transmission function can be derived: Equation ( 10) is solved by the method of iterations, which provides the analytical extension in a systematic way.A detailed description of this procedure is given in [35,36].The final result is


where   is the integral exponent.
Figures 1 and 2 illustrate the comparison between the exact numerical solution and the smallangle approximation model.

Reflection from multilayer structures
We consider reflection from multilayer structures using the angular distributions of electrons, elastically reflected from solids, as an example.In the literature, there are experimental data on the angular distributions of electrons reflected by homogeneous solids as well as by multilayer samples [9][10][37][38][39][40][41][42].
Consider a two-layer sample (see Fig. 3).In accordance with the described scheme, the reflection function for a two-layer sample can be represented as or, exploiting the one-speed approximation and the small angle approximation, , , , , , 1/ 1/ , , , , .
The computations for three-layer systems can be performed by using the relation, which is similar to equation ( 16), namely, , lm R  is computed by using equation ( 16) with the following index perturbation: 1 2; 2 3  .
Figure 4 shows the angular distributions of electrons, elastically reflected by the two-layer systems (the layer of Be on the Au substrate), computed by using the exact numerical solutions and the small-angle solution (16).

Figure 1
Figure 1 The reflection functions for the semi-infinite medium.The normal angle of incidence.The calculations are performed for the Henyey-Greenstein phase function.Plot (a) shows the influence of the non-linear term on the reflection function.Plot (b) illustrates the increase of error of the small-angle approximation with the phase function smoothness.The solid line is the numerical solution (MDOM), the dashed line is the small-angle approximation.The normal angle of incidence.The single scattering albedo is 0.67.Numbers with arrows are the asymmetry parameters of the Henyey-Greenstein phase function.

Figure 2 .
Figure 2. The transmission functions for the layers of thickness 0.1 tr l .Plot (a) shows the influence of the non-linear term.The solid line is the result with the non-linear term, while the dashed line shows the result without non-linear term taken into account.Plot (b) shows the increasing of the small-angle approximation error with the smoothness of the phase function.The solid line is the numerical solution (MDOM), the dashed line is the small-angle approximation.The angle of incidence is 45 °, the single scattering albedo is 0.54.Numbers with arrows are the asymmetry parameters of the Henyey-Greenstein phase function.

Figure 3 .
Figure 3.The two-layer model of reflection.The curved line corresponds to the reflection function, while the straight line is associated with the transmission function.