The radiative transport equation with heterogeneous cross-sections

TL;DR

This paper analyzes the radiative transport equation with heterogeneous cross-sections, providing explicit bounds on the integral operator, convergence rates, and solution estimates, with applications to nuclear reactor safety.

Contribution

It offers explicit estimates for the integral operator and solution of the RTE with heterogeneous coefficients, and applies these results to eigenvalue problems in reactor safety.

Findings

Bound on the integral operator norm explicit in cross-sections

Convergence rate estimate for source iteration

Eigenvalues in reactor safety are real and positive

Abstract

We consider the classical integral equation reformulation of the radiative transport equation (RTE) in a heterogeneous medium, assuming isotropic scattering. We prove an estimate for the norm of the integral operator in this formulation which is explicit in the (variable) coefficients of the problem (also known as the cross-sections). This result uses only elementary properties of the transport operator and some classical functional analysis. As a corollary, we obtain a bound on the convergence rate of source iteration (a classical stationary iterative method for solving the RTE). We also obtain an estimate for the solution of the RTE which is explicit in its dependence on the cross-sections. The latter can be used to estimate the solution in certain Bochner norms when the cross-sections are random fields. Finally we use our results to give an elementary proof that the generalised eigenvalue problem arising in nuclear reactor safety has only real and positive eigenvalues.

Equations

$$\Omega \cdot \nabla \psi (\mathbf{r},\Omega )+\sigma (\mathbf{r})\psi (\mathbf{r},\Omega )=\frac{{\sigma }_S(\mathbf{r})}{4\pi }{\int }_{{\mathbb{S}}^2}\psi (\mathbf{r},{\Omega }^{\prime })\text{d}{\Omega }^{\prime }+Q(\mathbf{r}),$$

$$\sigma (\mathbf{r})={\sigma }_S(\mathbf{r})+{\sigma }_A(\mathbf{r}),\forall \mathbf{r}\in V,$$

$$0<({\sigma }_S{)}_{\min }\le {\sigma }_S(\mathbf{r})\le ({\sigma }_S{)}_{\max },0<({\sigma }_A{)}_{\min }\le {\sigma }_A(\mathbf{r})\le ({\sigma }_A{)}_{\max },$$

$$0<(\sigma {)}_{\min }\le \sigma (\mathbf{r})\le (\sigma {)}_{\max },$$

$$\psi (\mathbf{r},\Omega )=0,\text{when}\mathbf{n}(\mathbf{r})\cdot \Omega <0,\mathbf{r}\in \partial V,$$

$$\mathcal{T}\psi (\mathbf{r},\Omega )=\Omega \cdot \nabla \psi (\mathbf{r},\Omega )+\sigma (\mathbf{r})\psi (\mathbf{r},\Omega )\text{and}\mathcal{P}(\cdot )=\frac{1}{4\pi }{\int }_{{\mathbb{S}}^2}(\cdot )\text{d}\Omega ,$$

$$\mathcal{T}\psi (\mathbf{r},\Omega )={\sigma }_S(\mathbf{r})\varphi (\mathbf{r})+Q(\mathbf{r}),$$

$$\mathcal{T}{\psi }^{i+1}(\mathbf{r},\Omega )={\sigma }_S{\varphi }^i(\mathbf{r})+Q(\mathbf{r}),\text{subject to}\text{eqref}eq:3DnonzeroBCs\text{and}{\varphi }^{i+1}:=\mathcal{P}{\psi }^{i+1}.$$

$$\mathcal{T}\psi (\mathbf{r},\Omega )=g(\mathbf{r}),\mathbf{r}\in V,\Omega \in {\mathbb{S}}^2,$$

$$\psi (\mathbf{r},\Omega )={\int }_0^{d(\mathbf{r},\Omega )}\exp (-\tau (\mathbf{r},\mathbf{r}-s\Omega ))g(\mathbf{r}-s\Omega )\mathrm{d}s,$$

$$d(\mathbf{r},\Omega )=\inf \{s>0:\mathbf{r}-s\Omega \notin V\}.$$

$$\varphi (\mathbf{r}):=(\mathcal{K}g)(\mathbf{r}),$$

$$(\mathcal{K}v)(\mathbf{r}):={\int }_{V}k(\mathbf{r},{\mathbf{r}}^{\prime })v({\mathbf{r}}^{\prime })\text{d}{\mathbf{r}}^{\prime },\text{with kernel}k(\mathbf{r},{\mathbf{r}}^{\prime }):=\frac{\exp \left(-\tau (\mathbf{r},{\mathbf{r}}^{\prime })\right)}{4\pi {\Vert \mathbf{r}-{\mathbf{r}}^{\prime }\Vert }_2^2}.$$

$$\displaystyle\phi^{i+1}\$$

$$\displaystyle\text{and}\quad e^{i+1}\$$

$$\varphi -\mathcal{K}{\sigma }_S\varphi =\mathcal{K}Q,$$

$${\Vert \mathcal{K}{\sigma }^{\ast }\Vert }_{L^2(V,\sigma )}\le {\Vert \frac{{\sigma }^{\ast }}{\sigma }\Vert }_{L^{\infty }(V)}.$$

$$\mathcal{L}:={\sigma }^{1/2}\mathcal{K}{\sigma }^{1/2}.$$

$$\displaystyle\psi({\bf r},\Omega)g({\bf r})\$$

$$=\frac{1}{2}\nabla \cdot (\Omega {\psi }^2(\mathbf{r},\Omega ))+\sigma (\mathbf{r}){\psi }^2(\mathbf{r},\Omega ).$$

$$\displaystyle\int_{V}\psi({\bf r},\Omega)g({\bf r})\mathrm{d}r\$$

$$\ge {\int }_V\sigma (\mathbf{r}){\psi }^2(\mathbf{r},\Omega ),$$

$$\displaystyle\langle g,\mathcal{K}g\rangle\ =\ \int_{V}\phi({\bf r})g({\bf r})\mathrm{d}r\$$

$$⟨g,\mathcal{L}g⟩\le \Vert g{\Vert }_{L^2(V)}^2,\text{for all}g\in L^2(V).$$

$$\mathcal{T}\psi ={\sigma }^{1/2}g,$$

$$\mathcal{P}(\Omega \cdot \nabla \psi )+\sigma \varphi ={\sigma }^{1/2}g,$$

$$\displaystyle\sigma^{1/2}\mathcal{K}\sigma^{1/2}g\$$

$$\displaystyle\langle g,\mathcal{L}g\rangle\$$

$$\frac{1}{4\pi }{\int }_{{\mathbb{S}}^2}{\int }_V{\sigma }^{-1/2}(\mathbf{r})g(\mathbf{r})\Omega \cdot \nabla \psi (\mathbf{r},\Omega )\mathrm{d}r\mathrm{d}\Omega .$$

$$\frac{1}{4\pi }{\int }_{{\mathbb{S}}^2}{\int }_V\left(\psi (\mathbf{r},\Omega )\Omega \cdot \nabla \psi (\mathbf{r},\Omega )\mathrm{d}\mathbf{r}\mathrm{d}\Omega +{\sigma }^{-1}(\Omega \cdot \nabla \psi (\mathbf{r},\Omega ){)}^2\right)\mathrm{d}\mathbf{r}\mathrm{d}\Omega$$

$$\text{mbox}\ge \frac{1}{8\pi }{\int }_{{\mathbb{S}}^2}{\int }_V\nabla \cdot ({\psi }^2(\mathbf{r},\Omega )\Omega )\mathrm{d}\mathbf{r}\mathrm{d}\Omega =\frac{1}{8\pi }{\int }_{{\mathbb{S}}^2}{\int }_{\partial V}{\psi }^2(\mathbf{r},\Omega )\Omega \cdot \mathbf{n}(\mathbf{r})\mathrm{d}s\mathrm{d}\Omega ,$$

$$\Vert \mathcal{L}g{\Vert }_{L^2(V)}^2=⟨g,{\mathcal{L}}^{2}g⟩\le \Vert g{\Vert }_{L^2(V)}^2,\text{for all}g\in L^2(V).$$

Full text

1 Introduction

In this note we present some elementary estimates for the steady-state mono-energetic Radiative (Boltzmann or Neutron) Transport Equation (RTE) with heterogeneous cross-sections. Although these are relatively straightforward to prove, and there is a huge literature on this topic, we were unable to find proofs of these precise results in the literature and so it seems useful to record them here. We emphasise that there are many excellent classical references for the general field discussed here - for example BelGla:70 ; LeMi:84 ; DaLi:12 . The estimates given here provide the tools needed to estimate the solution of the RTE explicitly in terms of the data (the so-called ‘cross-sections’). They also allow an explicit estimate for the rate of convergence of source iteration in terms of the scattering ratio. These estimates have recently proved essential for the rigorous analysis of uncertainty quantification techniques for the RTE GrPaSc:18 ; GrPaSc:18a ; Pa:18 . Most of the estimates presented here appeared in the University of Bath PhD theses of Fynn Scheben and Jack Blake Sch:11 ; Bl:16 ; GrSc:11 . An application to problems with random data appears in the more recent University of Bath thesis of Matthew Parkinson Pa:18 - see also GrPaSc:18 ; GrPaSc:18a .

For ${\bf r}\in V\subset\mathbb{R}^{3}$, where $V$ is a bounded convex spatial domain and $\Omega\in{\mathbb{S}^{2}}$, the unit sphere in 3D, and assuming isotropic scattering, the RTE source problem takes the form

[TABLE]

where $\nabla$ denotes the gradient with respect to ${\bf r}$, $Q$ is the source and ${\sigma}$, ${\sigma_{S}}$ are, respectively, the total and scattering cross sections that satisfy

[TABLE]

where ${\sigma_{A}}$ is the absorption cross section. All cross-sections are assumed to be pointwise bounded above and below on $V$ by strictly positive constants, i.e., for all   ${\bf r}\in V$,

[TABLE]

and thus

[TABLE]

which ensures that $\|{\sigma_{S}}/\sigma\|_{L^{\infty}(V)}<1$.

In (1), the angular flux $\psi$ is to be found, subject to given boundary conditions. Here we only consider the vacuum boundary condition on the inflow boundary:

[TABLE]

where $\mathbf{n}$ denoted the ourward normal from $V$.

Introducing the transport and averaging operators:

[TABLE]

(1) can be rewritten as

[TABLE]

where $\phi=\mathcal{P}\psi$ is called the scalar flux. This is to be solved, subject to (3),

Throughout, $L^{2}(V)$ will denote the space of square integrable functions on $V$, with inner product denoted $\langle\cdot,\cdot\rangle$. For any uniformly positive and bounded weight function $w$ on $V$, we will use $L^{2}(V,w)$ to denote the space of functions $v$, for which $\|v\|_{L^{2}(V,w)}^{2}:=\int_{V}|v({\bf r})|^{2}w({\bf r})\mathrm{d}r<\infty$. When $w\equiv 1$, this reduces to the standard $L^{2}$ norm $\|v\|_{L^{2}(V)}^{2}=\langle v,v\rangle$.

Before studying (1), (3) it is useful to first consider the ‘pure transport problem’ $\mathcal{T}\psi=g$, subject to boundary condition (3) and with source $g\in L^{2}(V)$. If the scalar flux $\phi$ were known, then (5) would allow computation of $\psi({\bf r},\Omega)$ for all ${\bf r}$ and any fixed $\Omega$ by solving a single transport equation with vacuum boundary condition. Such ‘transport sweeps’ are ‘easy’ operations both in theory (the solution of the transport problem can be written down using characteristics - see Lemma 1), and in numerical practice (e.g., when discontinuous Galerkin methods are used to discretise the transport operator, the resulting linear system can usually be solved by a single sweep through the elements). This motivates the use of the ‘source iteration’ for solving (5), which computes a sequence of approximations to $\phi$, starting with an initial guess $\phi^{0}$ and iterating by solving :

[TABLE]

The fact that this iteration is always well-defined and is equivalent to a fixed point iteration for a certain self-adjoint weakly singular integral operator is established in the following lemma.

Lemma 1

Let $g\in L^{2}(V)$ and consider the pure transport problem: Solve

[TABLE]

for $\psi({\bf r},\Omega)$, subject to boundary condition (3). This problem has a unique solution given by

[TABLE]

where \tau({\bf r},{\bf r}^{\prime})=\displaystyle{\int_{l({\bf r},{\bf r}^{\prime})}{\sigma}({\bf z})\;\textrm{d}l({\bf z})}\ (the integral of $\sigma$ along the line $l({\bf r},{\bf r}^{\prime})$ joining ${\bf r}$ and ${\bf r}^{\prime}$) and

[TABLE]

Moreover, the correponding scalar flux $\phi:=\mathcal{P}\psi$ can be expressed as

[TABLE]

where $\mathcal{K}$ is the integral operator defined by

[TABLE]

Proof

It is easy to show (using the method of characteristics) that the formula (7) provides the unique solution to the pure transport problem. This is well-known in the neutron transport literature, e.g. BelGla:70 ; PrLa:10 and $\tau$ is called the optical path length. Once (7) is established, the formula (8) is obtained by applying $\mathcal{P}$ to each side of (7) and rewriting the result using spherical polar coordinates.

Remark 1

An interesting observation is that, although the original transport problem is far from being self-adjoint, the integral operator $\mathcal{K}$ is self-adjoint with a positive kernel (and in fact is a positive definite operator, as we see below). Since $\mathcal{K}$ has a weakly singular kernel and the domain of integration is a bounded Euclidean domain, it is to be expected that the solution $\phi$ will have (weak) boundary singularities. This property has been analysed in the classical literature

• see, e.g., PiSc:83 and the references therein.

Returning to source iteration (6), we have the following simple corollary:

Corollary 1

The iterates $\phi^{i}$ and the errors $e^{i}:=\phi-\phi^{i}$ satisfy the equations

[TABLE]

Hence iteration (6) will converge if and only if there is a norm in which the operator $\mathcal{K}{\sigma_{S}}$ is a contraction. (Here we emphasise that $\mathcal{K}{\sigma_{S}}$ denotes the composition of the operator of multiplication by ${\sigma_{S}}$ with the integral operator $\mathcal{K}$.) In Theorem 2.1 we prove that this is the case in a certain weighted $L^{2}$ norm induced by the total cross-section ${\sigma}$. Then Corollary 2 provides the result on the convergence of source iteration.

In fact Theorem 2.1 has several other ramifications. Combining Lemma 1 with (5) we obtain that the scalar $\phi$ satisfies the second kind weakly singular integral equation

[TABLE]

and Theorem 2.1 (and the Banach lemma) then readily tells us that this equation has a unique solution and provides a bound on its norm explicit in the cross-sections (Corollary 3).

One reason for providing these results in this paper is that their proofs are hard to locate in the literature. Another reason is that they have direct relevance to the modern theory of unertainty quantification for the transport equation. When the cross-sections $\sigma_{S}$ and $\sigma$ are random fields, then both the error estimates for numerical methods for computing $\phi$ and also the rates of convergence of iterative methods for computing realisations of the scalar flux $\phi$ depend expilcitly on the cross-sections through the theorems presented here. This dependence is used explicitly in recent work on UQ for transport problems GrPaSc:18 ; Pa:18 ; GrPaSc:18a .

It is known that source iteration converges when solving the neutron transport equation with constant cross sections - see, e.g., (Sch:11, , Chapter 4). Ashby et. al (AshBroDorHin:95, , Section 4) prove a similar result with spatially dependent cross sections for a special discrete case. This work motivated us to consider a general proof in the heterogeneous case. The results here are for the underlying operator before discretization. Extension to general discretizations is a complicated question.

We will present the theory for the full 3D case where ${\bf r}\in V\subset\mathbb{R}^{3}$ and $\Omega\in\mathbb{S}^{2}$ the unit sphere in 3D, but the result also applies to the 2D reduction where ${\bf r}\in V\subset\mathbb{R}^{2}$ and $\Omega\in\mathbb{S}^{1}$ and to the case when space and angle are one-dimensional (the so called slab geometry case). Details of the proof in this case are given in Bl:16 .

2 The main result

Our main goal in this section will be to prove the following theorem.

Theorem 2.1

For any function $\sigma^{*}$ satisfying $0<\sigma^{*}_{\mathrm{min}}\leq\sigma^{*}({\bf r})\leq\sigma^{*}_{\mathrm{max}}$ for all ${\bf r}\in\overline{\Omega}$,

[TABLE]

(The left hand side of the inequality in (12) denotes the operator norm of $\mathcal{K}\sigma^{*}$ on the space $L^{2}(V)$, equipped with the weighted norm $\|\cdot\|_{L^{2}(V,\sigma)}$.)

The proof depends on several lemmas. In these it is useful to introduce the operator

[TABLE]

Lemma 2

The operators $\mathcal{K}$ and $\mathcal{L}$ are compact, self-adjoint and positive definite on $L^{2}(V)$.

Proof

First, $\mathcal{K}$ is compact on $L^{2}(V)$ because it is a weakly singular operator of potential type, see (KanAki:82, , p.332). To see the positive definiteness, let $g$ be an arbitrary function in $L^{2}(V)$ and let $\psi$ be the solution of $\mathcal{T}\psi=g$, subject to vacuum boundary conditions (3). Then

[TABLE]

Hence, integrating over $V$ and using the divergence theorem, we obtain

[TABLE]

where we used the vacuum boundary condition to get the final inequality in (13). Now introducing $\phi=\mathcal{P}\psi$, and recalling (8), we also have $\mathcal{P}\psi=\phi=\mathcal{K}g$. Applying $\mathcal{P}$ to each side of (13), we then obtain

[TABLE]

This proves the postitive definiteness of $\mathcal{K}$. Since $\mathcal{L}$ is a simple left and right scaling of $\mathcal{K}$ with the positive-valued function $\sigma^{1/2}$, the proof for $\mathcal{L}$ follows directly.

Our next result concerns an upper bound on $\mathcal{L}$.

Lemma 3

[TABLE]

Proof

In a variation of the proof of Lemma 2, let $g\in L^{2}(V)$, but this time let $\psi$ be the solution of

[TABLE]

subject to vacuum boundary conditions. Then set $\phi=\mathcal{P}\psi$, implying that $\mathcal{P}\psi=\phi=\mathcal{K}(\sigma^{1/2}g)$. This time, applying $\mathcal{P}$ directly to (14) and recalling that $g$ and $\sigma$ are both independent of $\Omega$, we get

[TABLE]

and so

[TABLE]

Multiplying each side of this relation by $g$ and integrating over $V$, we get

[TABLE]

Examining the second term on the right-hand side of (15), we see that this may be written

[TABLE]

Multiplying (14) by $\sigma^{-1}$, we obtain the formula $\sigma^{-1/2}g=\psi+\sigma^{-1}\Omega.\nabla\psi$, Using this in (16) and then the divergence theorem again, we see that (16) is

[TABLE]

the inner integral on the right-hand side being over the surface $\partial V$. This is non-negative because of the vacuum boundary conditions. Hence (16) is non-negative and combining this with (15), we obtain the result.

Lemma 4

[TABLE]

Proof

By (RieNag:55, , Chapter 104) the positive-definite self-adjoint operator $\mathcal{L}$ possesses a unique positive-definite self-adjoint square root, $\mathcal{L}^{1/2}$. Take any $g\in L^{2}(V)$. Then, using the self-adjointness of $\mathcal{L}^{1/2}$ and Lemma 3 (twice), we obtain

[TABLE]

as required.

Using the above results we are now in a position to prove the main theorem.

*Proof of Theorem 2.1. *    For any $g\in L^{2}(V,\sigma)$, we have $g\in L^{2}(V)$ and we can write

[TABLE]

Using Lemma 4, we then have

[TABLE]

and the result follows.

Remark 2

Although we have given the proof here only in the 3D case, the same result holds for classical 2D and 1D model problems. For example in the 1D “slab geometry” case formulated on the unit interval, the transport equation is

[TABLE]

where $x\in(0,1)$ and $\mu\in(-1,1)$. The counterpart of the intergral operator $\mathcal{K}$ is

[TABLE]

with $\tau$ denoting the optical path and $E_{1}$ the exponential integral. The counterpart of Theorem 2.1 for this case is proved using almost identical arguments to those given above (see, e.g., Bl:16 ).

2.1 Some applications of Theorem 2.1

Convergence of source iteration

From Theorem 2.1 and Corollary 1 we immediately have the following result on the convergence of source iteration.

Corollary 2

Under the definitions above we have:

[TABLE]

Since, by (2), $\left\|\dfrac{{\sigma_{S}}}{{\sigma}}\right\|_{L^{\infty}(V)}<\ 1$, we have $e^{i}\rightarrow 0$ as $i\rightarrow\infty$.

Remark 3

A stronger estimate than (17) can be obtained in the case of constant cross-sections. In Bl:16 it was shown that on a spatial domain $V$ with diameter $d$

[TABLE]

So on small domains source iteration can still converge rapidly, even if the scattering ratio is close to $1$.

Data-explicit estimates for the RTE

The next corollary gives data-explicit estimates for the pure transport problem and for the RTE.

Corollary 3

(i) Consider the pure transport problem $\mathcal{T}\psi=g$ with vacuum boundary conditions, as in Lemma 1, and let $\phi$ be the corresponding scalar flux. Then

[TABLE]

(ii) Consider the RTE (1) with vacuum boundary conditions (3) and let $\phi$ be the corresponding scalar flux. Then

[TABLE]

Proof

(i) By Theorem 2.1, we have $\phi=\mathcal{K}g$, so

[TABLE]

(using Lemma 4), and this yields the result.

(ii) By Theorem 2.1, the operator $\mathcal{K}{\sigma_{S}}$ is a contraction on $L^{2}(V,\sigma)$ with norm bounded by $\|{\sigma_{S}}/\sigma\|_{L^{\infty}(V)}<1$. Hence, by (11) and the Banach Lemma, we can write $\phi=(I-\mathcal{K}\sigma_{S})^{-1}\mathcal{K}Q$, with

[TABLE]

Writing $\sigma^{1/2}\mathcal{K}Q=\mathcal{L}(\sigma^{-1/2}Q)$ and proceeding as in part (i), we obtain (ii).

Remark 4

(i) The bounds in Corollary 3 provide the mechanism for estimating the flux $\phi$ in appropriate Bochner norms when the data $\sigma,\sigma_{A},\sigma_{S}$ are random fields, i.e., when we wish to quantify how uncertainty in data propagates to uncertainty in the fluxes or in the criticality (see the next subsection). Integrability in probability space of the right-hand sides in each of the estimates (i) or (ii) above immediately implies the same integrability properties for the resulting flux. In Pa:18 ; GrPaSc:18 ; GrPaSc:18a this is worked out in detail and the theory of multilevel Monte Carlo methods for computing quantities of interest is presented for the RTE in one and two-dimensional models.

Spectral properties of the RTE

In the study of nuclear reactor stability, one is concerned with the eigenvalues $\lambda$ of the generalised eigenproblem:

[TABLE]

with vacuum boundary condition. Here (in this simplified model problem), $\sigma_{F}$ is the fission cross-section which is also assumed bounded above and below on $V$ by positive constants, and now

[TABLE]

In fact one is concerned with the fundamental eigenvalue of (18), the smallest in absolute value. The reactor is stable and efficient provided the fundamental eigenvalue is close to $1$. In this case the neutrons produced by fission balance the neutrons lost by scattering and absorption.

It is a not completely obvious fact that the spectrum of the problem (18) is in fact discrete, real and bounded below by a positive number. This fact can be obtained from the elementary properties which we have derived above.

Corollary 4

The eigenvalues of problem (18) are real and positive.

Proof

Let $(\lambda,\psi)$ be an eigenpair of (18). Then, as in (11),

[TABLE]

Multiplying through by ${\sigma_{S}}^{1/2}$ and setting $v=\sigma_{S}^{1/2}\phi$, we have

[TABLE]

where $\mathcal{L}_{\sigma_{S}}:=\sigma_{S}^{1/2}\mathcal{K}\sigma_{S}^{1/2}$. Now, since $\mathcal{L}_{\sigma_{S}}v=(\sigma_{S}/\sigma)^{1/2}\mathcal{L}((\sigma_{S}/\sigma)^{1/2}v)$, we can use Lemma 4 to obtain

[TABLE]

Since $\|{\sigma_{S}}/\sigma\|_{L^{\infty}(V)}<1$, $\mathcal{L}_{{\sigma_{S}}}$ is a contraction on $L^{2}(V)$, and $I-\mathcal{L}_{\sigma_{S}}$ is an invertible operator. Thus $\lambda$ cannot vanish in (19) (since if $\lambda=0$, then $v=0$, which implies $\phi=0$ and hence by (18), $\psi=0$). Hence (19) is equivalent to

[TABLE]

where $\mathcal{M}=(I-\mathcal{L}_{\sigma_{S}})^{-1}\mathcal{L}_{\sigma_{S}}$. It is easy to see that $\mathcal{M}$ is self-adjoint and compact and has eigenmalues $(1-\mu)^{-1}\mu$ where $\mu$ denotes an eigenvalue of $\mathcal{L}_{\sigma_{S}}$. Since $\mu$ is always positive and less than $1$ it follows that the eigenvalues of $\mathcal{M}$ are all positive and so $\mathcal{M}$ is positive definite. Then setting $w=(\sigma_{F}/\sigma_{S})^{1/2}v=\sigma_{F}^{1/2}\phi$, we have

[TABLE]

where

[TABLE]

Since $\mathcal{N}$ is also self-adjoint and positive definite, the result follows.

Remark 5

(i) A more sophisticated argument based on the Krein-Rutman theorem can be used to show that the fundamental eigenvalue is simple with a positive eigenfunction.

(ii) The eigenvalue problem is discussed in detail in the fundamental reference DaLi:12 . (iii) The quantity $1/\lambda$ is called “$k-$effective” in the nuclear engineering literature.

Acknowledgement  We thank Professor Paul Smith (Wood plc., Poundbury, Dorset, UK) for many helful discussions over many years’ collaboration, and for supporting the PhD theses of Fynn Scheben, Jack Blake and Matt Parkinson, whose work is partially reported here.