Spherical harmonics, method of

A method for obtaining an approximate solution of a kinetic equation by the decomposition of the phase density of the particles into a finite sum of spherical functions in arguments that define the direction of the velocity of the particles (see [1]). The method is widely used in the solution of problems of neutron physics.

In one-dimensional plane geometry the stationary, integro-differential, kinetic transport equation (given an isotropic scattering of the particles)

$$ \tag{1 } \mu \frac{d \psi ( x, \mu ) }{dx} + \psi ( x, \mu ) = \ \frac{c}{2} \int\limits _ { -1}^ { +1} \psi ( x, \mu ^ \prime ) d \mu ^ \prime $$

is replaced by an approximate system of differential equations for $ \widetilde \psi _ {n} ( x) $— the approximate values of the Fourier coefficients

$$ \tag{2 } \psi _ {n} ( x) = \int\limits _ { - 1} ^ { +1 } \psi ( x, \mu ) P _ {n} ( \mu ) d \mu ,\ \ n = 0 \dots 2N- 1. $$

A system of the form

$$ \tag{3 } n \frac{d \widetilde \psi _ {n-1} ( x) }{dx} + ( 2n+ 1- c \delta _ {n0} ) \widetilde \psi _ {n} ( x) + ( n+ 1) \frac{d \widetilde \psi _ {n+1} ( x) }{dx} = 0 \ \ $$

arises, under the condition

$$ \tag{4 } \widetilde \psi _ {2N} ( x) \equiv 0. $$

Here, $ \psi ( x, \mu ) $ is the phase density of the particles that are scattered in the matter, $ c $ is the average number of secondary particles arising from one act of interaction with the particles of matter, and $ P _ {n} ( \mu ) $ is the Legendre polynomial of degree $ n $( cf. Legendre polynomials). System (3) defines the $ P _ {2N-1} $- approximation of the method of spherical harmonics for equation (1). The approximate value of the phase density is

$$ \tag{5 } \widetilde \psi ( x, \mu ) = \ \sum _ {n=0} ^ { 2N-1} 2n+ \frac{1}{2} \widetilde \psi _ {n} ( x) P _ {n} ( \mu ). $$

For equation (1), the typical boundary conditions take the form

$$ \tag{6 } \left . \begin{array}{ll} \psi ( 0, \mu ) = 0 & \textrm{ for } 0 < \mu \leq 1, \\ \psi ( h, \mu ) = 0 & \textrm{ for } - 1 \leq \mu < 0. \\ \end{array} \right \} $$

Boundary conditions of this type exist in, for example, the problem in neutron physics of the critical conditions of a layer of thickness $ h $ with free surfaces $ x= 0 $ and $ x= h $( boundaries with vacuum). In this problem, a positive solution of (1) and (6), as well as the eigenvalue $ c $, must be found.

In the method of spherical harmonics, it is natural to replace (6) by

$$ \tag{7 } \left . \begin{array}{c} \int\limits _ { 0 } ^ { 1 } P _ {n} ( \mu ) \widetilde \psi ( 0, \mu ) d \mu = 0, \\ \int\limits _ { -1} ^ { 0 } P _ {n} ( \mu ) \widetilde \psi ( h, \mu ) d \mu = 0, \\ n = 0 \dots 2N- 1. \\ \end{array} \right \} $$

However, this approach involves twice as many conditions than are required for a partial solution of (3). In practice, a different set of values has been tested in (7). The best result is obtained by conditions with $ n= 2k+ 1 $, $ k = 0 \dots N- 1 $. For the single-velocity transport equation of a general type, the Vladimirov variational principle leads to a system of equations of the method of spherical harmonics and the given boundary conditions (when choosing the experimental functions in the form of a linear combination of spherical harmonics). In three-dimensional geometry, the boundary conditions can be written in the form

$$ \tag{8 } \left . \int\limits _ {( \Omega n) < 0 } ( \pmb\Omega \mathbf n ) \widetilde \psi ( \mathbf r , bold \Omega ) Y _ {2k,i} ( \pmb\Omega ) d \pmb\Omega \right | _ {r \in \Gamma } = 0, $$

$$ i = 0, \pm 1 \dots \pm 2k; \ k = 0 \dots N- 1. $$

Here $ \mathbf r $ is a vector of space coordinates, $ \pmb\Omega $ is a unit vector of the velocity of the particle with spherical coordinates $ \{ \theta , \phi \} $, $ \mathbf n $ is the unit vector of the exterior normal to the piecewise-smooth surface $ \Gamma $ that bounds the convex domain in space in which the problem is solved,

$$ Y _ {k,i} ( \pmb\Omega ) = P _ {k} ^ {(|} i|) ( \cos \theta ) \sin \ | i | \phi ,\ \ i = - k \dots - 1, $$

$$ Y _ {k,i} ( \pmb\Omega ) = P ^ {(i)} ( \cos \theta ) \cos i \phi ,\ i = 0 \dots k, $$

are the spherical functions, $ P _ {k} ^ {(i)} ( \mu ) $ are the associated Legendre functions of the first kind ( $ P _ {k} ^ {(0)} ( \mu ) = P _ {k} ( \mu ) $ are the Legendre polynomials).

The lowest approximations $ ( P _ {1} , P _ {3} ) $ of the method of spherical harmonics are widely used in solving problems of neutron physics and give good results well away from the boundaries of the domain, sources and strong neutron absorbers. The theory of growth is also constructed in a $ P _ {1} $- approximation. The generalized solution of the method of spherical harmonics converges to the solution of the transport equation when $ N \rightarrow \infty $( see [2]). The rate of the convergence $ \widetilde \psi ( x , \mu ) \rightarrow \psi ( x, \mu ) $ is easily estimated, by comparing the integral equations for $ \widetilde \psi _ {0} ( x) $ and $ \psi _ {0} ( x) $, i.e. by estimating the proximity of their kernels. Equation (1) with the boundary conditions (6) reduces to an integral equation with the kernel

$$ E _ {1} ( | x- s | ) = \int\limits _ { 0 } ^ { 1 } \frac{e ^ {- | x- s | / \mu } } \mu d \mu . $$

System (3), given boundary conditions analogous to (6),

$$ \tag{9 } \left . \begin{array}{ll} \widetilde \psi ( 0, \mu _ {i} ) = 0 & \textrm{ for } 0 < \mu _ {i} < 1 , i = 1 \dots N, \\ \widetilde \psi ( h, \mu _ {i} ) = 0 & \textrm{ for } - 1 < \mu _ {i} < 0, i = 1 \dots - N, \\ \end{array} \right \} $$

where $ \mu _ {i} $ are the roots of $ P _ {2N} ( \mu ) $, reduces to an integral equation with kernel

$$ \widetilde{E} _ {1} ( | x- s | ) = \sum_{i=1}^ { N } a _ {i} \frac{e ^ {- | x- s | / \mu _ {i} } }{\mu _ {i} } , $$

where $ a _ {i} $ are the weights of the Gauss quadrature formula for the system of nodes $ - 1 < \mu _ {i} < 1 $. The singularity of the function $ e ^ {- | x- s | / \mu } / \mu $ leads to a slow convergence for large $ N $:

$$ \int\limits _ { 0 } ^ { h } | E _ {1} ( x) - \widetilde{E} _ {1} ( x) | dx < \frac{\textrm{ const } }{N} ,\ \ | \psi _ {0} ( x) - \widetilde \psi _ {0} ( x) | < \frac{\textrm{ const } }{N} . $$

The approximate eigenvalue converges to an exact eigenvalue at a rate of $ 1/N ^ {2} $.

The boundary conditions (9) arise naturally in the solution of the kinetic equation by the method of discrete ordinates, which consists of replacing (1) by the approximate system

$$ \tag{10 } \mu _ {i} \frac{\partial \widetilde \psi tilde _ {i} ( x) }{\partial x } + \widetilde \psi tilde _ {i} ( x) = \ \frac{c}{2} \sum _ { i=- N }^ { N } a _ {i} \widetilde \psi tilde _ {i} ( x). $$

The method of discrete ordinates in one-dimensional geometry is equivalent to the method of spherical harmonics (see [3]), since the system (10) can be obtained from (3) by means of a linear transformation of the unknown functions:

$$ \widetilde \psi _ {n} ( x) = \ \sum _ { i=- N} ^ { N } a _ {i} \widetilde \psi tilde _ {i} ( x) P _ {n} ( \mu _ {i} ). $$

However, in multi-dimensional problems, the method of spherical harmonics in the lowest approximations is more accurate than the method of discrete ordinates.

References

[1]	G.I. Marchuk, V.I. Lebedev, "Numerical methods in the theory of neutron transport" , Harwood (1986) (Translated from Russian)
[2]	U.M. Sultangazin, "Convergence of the method of spherical harmonics for the non-stationary transport equation" USSR Comp. Math. Math. Phys. , 14 : 1 (1974) pp. 165–176 Zh. Vychisl. Mat. i Mat. Fiz. , 14 : 1 (1974) pp. 166–178
[3]	R.D. Richtmeyer, K.W. Morton, "Difference methods for initial-value problems" , Wiley (Interscience) (1967)

References

[a1]	G.I. Marchuk, "Methods of numerical mathematics" , Springer (1975) (Translated from Russian)
[a2]	B. Davison, J.B. Sykes, "Neutron transport theory" , Clarendon Press (1957)
[a3]	K.M. Case, P.F. Zweifel, "Linear transport theory" , Addison-Wesley (1967)