Multi-species neutron transport equation

The Neutron Transport Equation (NTE) describes the flux of neutrons through inhomogeneous fissile medium. Whilst well treated in the nuclear physics literature (cf. [9, 27]), the NTE has had a somewhat scattered treatment in mathematical literature with a variety of different approaches (cf. [8, 25]). Within a probabilistic framework it has somewhat undeservingly received little attention in recent years; nonetheless, probabilistic treatments can be found see for example [19, 26, 24, 29, 4, 3]. In this article our aim is threefold. First we want to introduce a slightly more general setting for the NTE, which gives a more complete picture of the different species of particle and radioactive fluxes that are involved in fission. Second we consolidate the classical c0-semigroup approach to solving the NTE with the method of stochastic representation which involves expectation semigroups. Third we provide the leading asymptotic of our multi-species NTE, which will turn out to be crucial for further stochastic analysis of the NTE in forthcoming work [6, 5]. The methodology used in this paper harmonises the culture of expectation semigroup analysis from the theory of stochastic processes against c0-semigroup theory from functional analysis. In this respect, our presentation is thus part review of existing theory and part presentation of new research results based on generalisation of existing results.


Introduction
The neutron transport equation (NTE) describes the flux of neutrons across a directional planar cross-section in an inhomogeneous fissile medium (typically measured is number of neutrons per cm 2 per second). As such, flux is described as a function of time, t, Euclidian location, r ∈ R 3 , direction of travel, ∈ S 2 , speed c > 0 (and hence velocity υ = c ), and neutron energy, E ∈ R. It is not uncommon in the physics literature, as indeed we shall do here, to assume that energy is a function of velocity (E = m|υ| 2 /2), thereby reducing the number of variables by one. This allows us to describe the dependency of flux more simply in terms of time and, what we call, the configuration variables (r , υ) ∈ D × V where D ⊆ R 3 is a smooth, open, connected and bounded domain of concern such that ∂ D has zero Lebesgue measure and V is the velocity space, which can now be taken to be V = {v ∈ R 3 : υ min < |v| < υ max }, where 0 < υ min < υ max < ∞.
Before stating the NTE, let us remind the reader of some elementary nuclear physics, which is required to describe the evolution of neutron flux. In the most basic of flux models, there are essentially only four processes at the level of the atomic nuclei which contribute to the evolution of neutron flux.
The first is spontaneous neutron emission from unstable nuclei. This comes from radioactive isotopes whose nuclei are excited. They cause what is known as non-transmutation emissions, in which a neutron is ejected with an escape velocity (neutron emission), or, conversely, what are called transmutation emissions in which the nucleus instantaneously fragments into two or more nuclei (spontaneous fission) with a range of possible masses, emitting one or more neutrons with escape velocities in the process.
The second process pertains to neutron scattering. This is where a neutron travelling with a given velocity passes in close proximity to an atomic nucleus, which, in our model, results in an instantaneous change of velocity.
The third process is neutron-induced fission. This is the classical setting in which a neutron travelling with a given velocity strikes an atomic nucleus sending it into an excited state, from which it instantaneously fragments into two or more nuclei, simultaneously releasing one or more neutrons.
The fourth and final process is neutron capture. In this setting, a neutron travelling with a given velocity strikes an atomic nucleus, but instead of causing nuclear fission, it is absorbed into the nucleus. It can also be the case that neutrons decay into other subatomic particles, and thus disappear from the system. To all intents and purposes, we can treat this as neutron capture.
When modelling the transmission of neutrons in a fissile material, those neutrons which have been released from nuclei are known as prompt neutrons.
With more advanced modelling, one can also take account of the fact that some of the processes described above can also involve other types of nuclear emissions, often in addition to neutrons. These include alpha and beta particles and gamma radiation. Whilst the former two are not sufficiently energetic to cause fission, sufficiently energetic gamma rays are able to induce fission.
Spontaneous fission and neutron-induced fission can also produce what are known as delayed neutrons. These are neutrons released from a fission product (isotope) some time after fission has occurred. In terms of modelling, they are spontaneous neutron emissions which occur at the site of neutron-induced fission but at a moment later in time. Delayed neutrons are only in a delayed state until they are released after which they are considered as prompt neutrons.
We refer to models which take account of the full range of flux profiles as multi-species models.

Neutron Transport Equation
Let us now write down the basic neutron transport equation (prompt neutrons only), which has been widely considered in a variety of physics and engineering literature (cf. [8,28], to name but two classical references), and somewhat more sporadically studied in the mathematical literature. See [6,17,25] for the three most authoritative mathematical texts in more recent times, as well as e.g. [3,4,13,22,26] for some of the rarer examples of the probabilistic treatment of the NTE.
It is normal to assume that all quantities are uniformly bounded away from infinity. It is also usual to assume the additional boundary conditions ⎧ ⎨ ⎩ 0 (r , υ) = g(r , υ) for r ∈ D, υ ∈ V , t (r , υ) = 0 f o rt ≥ 0 and r ∈ ∂ D if υ · n r < 0, (2.2) where n r is the outward facing normal of D at r ∈ ∂ D and g : D × V → [0, ∞) is a bounded, measurable function which we will later assume has some additional properties. Roughly speaking, as the forward equation describes where particles could have evolved from in order to contribute to the current configuration, this boundary condition means that particles from outside the domain with incoming velocity are not taken into account. The second of the above two boundary condition is sometimes written t | ∂ D − = 0, where It is also usual to set Q = 0 when considering a rector with a multiplying medium, as the resulting fission will overwhelm the radioactive source term. The notion of a solution of the form (2.1) turns out to be too strong to expect to make mathematical sense of it. This is predominantly due to the non-diffusive nature of the equation, in particular the non-local nature of the scattering and fission operators as well as regularity issues on the domain D × V in relation to continuity properties of e.g. the operator υ · ∇. It is much more natural to look for solutions that belong to e.g. an appropriate L 2 space. This is, moreover, helpful when looking to understand (2.1) as a backwards equation, rather than a forwards equation.
With some rearrangements, the components of (2.1) separate into transport, scattering and fission. Specifically, such that all operators are defined on D × V and their action is zero otherwise. Let us momentarily consider the operator on the right-hand side of (2.1) as acting on L 2 (D × V ), the space of square integrable functions on D × V , and write for the associated inner product. Note that, for f , g ∈ L 2 (D × V ) such that both υ · ∇ f and υ · ∇g are well defined as distributional derivatives, which are also in the space L 2 (D × V ), with g respecting the second of the boundary conditions in (2.2), we can verify with a simple integration by parts that, for υ ∈ V , providing we insist that f respects the boundary f (r , υ) = 0 for r ∈ ∂ D if υ · n r > 0. Moreover, Fubini's theorem also tells us that, for example, with f , These computations tell us that, for f , g ∈ L 2 (D × V ) such υ · ∇g and υ · ∇ f are well defined in the distributional sense and, moreover, that g(r , υ) = 0 for r ∈ ∂ D if υ · n r < 0, and for f ∈ where now we identify the transport, scattering and fission operators as such that all operators are defined on D × V with zero action otherwise. The reader will immediately note that, although the terms in the sum F are identifiable as the adjoint of the terms in the sum → T + → S + → F, the same can not be said for the individual 'T', 'S' and 'F' operators. That is to say, the way we have grouped the terms does not allow us to say that e.g. ← T is the adjoint operator to → T and so on. The reason for this difference in grouping of terms lies with how one reads the operators in terms of infinitesimal generators as a probabilist. Although this will not make any difference in the analysis of this paper, we keep to this notation for the sake of consistency with further related articles which offer a probabilistic perspective on the backwards NTE; see [5,12,14].
Roughly speaking, ← T, with an appropriately defined domain, is the generator of the rather simple Markov process consisting of a deterministic motion with velocity υ, i.e. transport due to pure advection, with killing on exiting the domain D. Similarly, with an appropriately defined domain, the operator ← S is the generator corresponding to scattering, in which a particle travelling with velocity υ at position r is removed at rate σ s and replaced by a new particle at r with velocity υ chosen with probability π s (r , υ, υ )dυ . Taking advantage of the fact that V π s (r , υ, dυ )dυ = 1 we can also write and also note that it takes the classical form of a difference operator. Finally ← F is the generator action of a fission even in which a particle travelling with velocity υ at position r is removed at rate σ f and replaced by an average number of particles π f (r , υ, υ )dυ moving onwards from r with velocity υ . This leads us to the so called backwards neutron transport equation (which is also known as the adjoint neutron transport equation) given by with additional boundary conditions ⎧ ⎨ ⎩ ψ 0 (r , υ) = g(r , υ) for r ∈ D, υ ∈ V , ψ t (r , υ) = 0 f o rt ≥ 0 and r ∈ ∂ D if υ · n r > 0. (2.7) Similarly to previously, the second of these two conditions is often written ψ t | ∂ D + = 0, where ∂ D + := {(r , υ) ∈ ∂ D × V : υ · n r > 0}. The NTE has played a prominent role in real-world modelling and, for many years, has found a home in commercial software which is used in the nuclear safety industry. In particular, this is most prominent in the modelling and design of environments which are exposed to radioactive material, from nuclear reactor cores and hospital equipment, through to equipment used to irradiate produce that is sold in supermarkets, thereby prolonging its shelflife. More recently, with the notion of human interplanetary space exploration becoming less of a sci-fi fantasy and more of a fast approaching reality, an understanding of how long-lasting and compact nuclear power sources, for e.g. Moon or Mars bases has become increasingly important. Figure 1 illustrates a typical geometrical model of a reactor core rod, cladding and outer shielding. 2 The structural design of such a reactor can easily be stored as virtual environment (i.e. storing the coordinates of the different geometrical domains and the material properties in each domain) with around 150 MB of data, on to which extensive data libraries of numerical values for the respective quantities σ s , σ f , π s , π f can be mapped. (It is an otherwise little known fact that countries which are heavily invested in nuclear power, such as the UK, USA, France, China, etc., are all in possession of such numerical libraries of cross sections, which have been carefully built up over decades.) One of the principal ways in which neutron flux is understood is to look for the leading eigenvalue and associated ground state eigenfunction. Roughly speaking, this means looking for an associated triple of eigenvalue λ ∈ R, non-negative right eigenfunction ϕ : As such, this introduces the notion of fissile stability, in particular in the case that λ = 0. This is naturally the desired scenario 3 for a nuclear reactor.
In the physics literature, it is thus often understood that, to leading order, the NTE (2.6) is solved in the approximate sense ψ t (r , υ) = e λt g,φ ϕ(r , υ) + o(e λt ), t ≥ 0. (2.8) Note that the scenario that λ > 0 is obviously to be avoided in practice as this would correspond to a set-up that could result in exponential growth in fission.
The approximation (2.8) can be seen as a functional version of the Perron-Frobenius Theorem and has given rise to a number of different numerical methods for estimating the value of the eigenvalue λ as well as the eigenfunctions ϕ andφ. One approach pertains to the discretisation of (2.1) followed by the use numerical analytic methods; see [31]. Another pertains to the previously alluded to identification of the solution to the NTE as the linear semigroup of a Markov branching process, which in turn implies Monte Carlo methods involving the simulation of the aforesaid branching process. Such methods are computationally expensive, as branching processes, being tree-like structures, are complex to simulate, e.g. from the point of view of parallelisation. In related papers to this one, we will discuss a new Monte Carlo approach to the NTE based on some of the stochastic analysis we deal with in this article as well as in related work undertaken by the authors of this paper; see [5,12,14].
The aim of this paper is manifold. First and foremost, we aim to reposition the theory of the NTE into a contemporary probabilistic setting. We will do this by explaining a precise relationship between the NTE and a two different families of Markov processes via Feynman-Kac type formulae. Indeed, this article is one of a cluster of forthcoming pieces of work, which take a new and predominantly probabilistic point of view of the NTE; cf. [5,12,14]. Next we want to introduce the notion of the (multi-species) NTE into the literature, which generalises (2.1) by simultaneously modeling the flux of all species of particles and radiation involved in the process of nuclear fission. In doing so we will show that, just as in the classical setting, one may develop the notion of a lead eigenvalue and eigenfunction, which is an important part of describing fissile stability. As such, the current article is part review of existing theory and part presentation of new research results based on generalisation of existing results.
Together with the accompanying papers [5,12,14], we believe that the probabilistic perspective presented here, i.e. coupling the solutions to the NTE with averaging procedures of certain Markov processes, opens up the possibility of many questions that can be considered at depth in the arena of stochastic analysis and Monte Carlo algorithms, which are currently missing from the literature. Indeed, returning to the kind of environments seen in Fig. 1, there are many questions concerning how to analyse and numerically generate the leading eigenfunctions and eigenvalue to a reasonable degree of precision. Such questions might include: What is the connection of the eigendecomposition discussed in this paper and e.g. R-theory or the theory of general Harris recurrence for stochastic processes (cf. [9,23,24]

Organisation of the Paper
In the next section, we give a brief overview of the key mathematical literature for the NTE. (Note we do not stray beyond mathematical literature, as the physics and engineering literature is significantly more expansive.) Thereafter in Sect. 5, we introduce the multi-species NTE (MNTE) and its rigorous formulation, existence, uniqueness and asymptotics in the setting of an abstract Cauchy problem. In particular, we show how the unique solution is identified as a c 0 -semigroup in the appropriate L 2 space. In Sect. 6, we introduce a spatial branching process that is constructed using the cross sections that appear in the NTE to describe its stochastic evolution. Here we introduce its expectation semigroup. In Sect. 7, we provide a second stochastic representation to the expectation semigroup introduced in the previous section via a classical method of the many-to-one formula.
Ideally, we would like to claim that the expectation semigroup discussed in Sects. 6 and 7 agree with the c 0 -semigroup introduced in Sect. 5 (its formal definition appearing just above Theorem 5.2). This is particularly desired as it forms the foundations of how Monte Carlo simulation of the physical process can be used to develop a numerical solution to the MNTE. In Sect. 8, we consolidate the two notions of semigroup and show that there is partial agreement in an appropriate sense. As far as we are aware, this is a point which is currently not clearly discussed in the literature. Finally we end the paper with a proof of one of the main theorems in Sect. 6 which provides the asymptotic behaviour of the solution to the MNTE in terms of the lead eigenfunction. This is a new result in the multi-species setting in the sense that we have allowed for multiple types of prompt emissions (both particles and radioactive emissions) rather than the case of only one type of prompt emission dealt with in [25]; we also allow for multiple types of delayed emissions (that is, emissions that are pre-emptively held in an unstable radioactive isotope product from an earlier fission event). Our proof nonetheless takes inspiration from the classical approach of [6,25], and remains loyal to the techniques there.

Historical Remarks on the Mathematical Treatment of the NTE
Classical texts such as Davison and Sykes [8] were once hailed as a bible of mathematical knowledge during the 1950s post Manhattan Project era when rapid technological advances lead to the construction of the very first nuclear reactors driving commercial power stations. Around this time, there was an understanding of how to treat the NTE in special geometries and also by imposing an isotropic scattering and fission, see for example Lehner [18] and Lehner and Wing [19,20]. It was also understood quite early on that the natural way to cite the NTE is via the linear differential transport equation associated to a suitably defined operator on a Banach space. Moreover, it was understood that in this formulation, a spectral decomposition should play a key role in representing solutions, see e.g. Jörgens [15], Pazy and Rabinowitz [29]. This notion was promoted by the work of R. Dautray and collaborators, who showed how c 0 -semigroups form a natural framework within which one may analyse the existence and uniqueness of solutions to the NTE; see [6,7]. Moreover, a similar approach has also been pioneered by Mokhtar-Kharroubi [25].
The probabilistic interpretation of the NTE was appreciated from the very first treatments of the NTE (see e.g. [8] and references therein, as well as Bell [2]). Indeed, the physical description of nuclear fission, when governed by basic principles, allowing for additional randomness, is nothing more than a branching Markov process. Numerous derivations of the NTE from this perspective can be found in the literature to various degrees of rigour; see e.g. Bell [2], Mori et al. [26], Pazy and Rabinowitz [30], Lewins [21] and Pázsit and Pál. [28].
A more modern treatment of the probabilistic representation through Feynman-Kac expectation semigroups and the connection to the theory of Markov diffusions is found in Dautray et al. [7]. A purely probabilistic can be found in Lapeyre et al. [17]. See also the accompanying papers to this one [5,12,14].
We finish this section by noting that there is a body of literature that pertains to the numerical analysis of the NTE. Recent work in this field, including the notion of uncertainty quantification, can be found in e.g. [22,27,31]. See also references therein.

Multi-species (Backwards) Neutron Transport Equation
In the following discussion, rather than talk about typed particles, we prefer to say typed 'emissions' as the different types correspond to particles, electromagnetic rays (e.g. gamma rays) and isotopes (which are considered to be carriers for delayed emissions).
Let us now introduce an advanced version of the NTE, which takes account of both nontransmutation emissions as well as transmutation emissions, in particular, allowing for the inclusion of all types of emissions, prompt neutrons, delayed neutrons, alpha, beta and gamma emissions etc. An important feature (and arguably a restriction) of our model is that only prompt neutrons can produce delayed emissions.
In order to keep track of the various emission types, we define the type space I := {1, . . . , m} for some m ∈ N, ordered such that Finally, the set of admissible velocities for each of the types i can be embedded within a common space V = {υ ∈ R 3 : υ min ≤ |υ| ≤ υ max }, with 0 < υ min ≤ υ max < ∞). We now consider the flux, ψ t (i, r , υ) of type i emissions through a given region r ∈ D with velocity υ ∈ V at time t ≥ 0. We are interested in the so called multi-species neutron transport equation (MNTE) which takes the form for prompt emissions i = 1, . . . , , whereas, in the case of delayed emissions, i = + 1, . . . , m satisfies which is of a simple form because it describes only how these emissions are held in a suspended state (no advection) before being converted back to prompt emissions. Similarly to before, have the following interpretation: : the rate at which scattering occurs for a type i emission with incoming velocity υ, : the rate at which fission occurs for a type i emission with incoming velocity υ, : the sum of the rates σ i f + σ i s and is known as the total cross section for a type i emission, π i s (r , υ, υ )dυ : the scattering yield at velocity υ from incoming velocity υ for a type i emission, satisfying V π i s (r , υ, υ )dυ = 1, π i, j f (r , υ, υ )dυ : the average type j yield at velocity υ from fission with incoming velocity υ for a type i emission satisfying : the average type j (unstable) isotope yield from a fission event due to a type 1 particle with incoming velocity υ, λ i : the decay rate for a type i isotope.
There are a number of assumptions about the many cross sections that appear in the above equations that will remain in force throughout the remainder of this text. Assumption 5.1 All cross sections are non-negative, measurable and uniformly bounded from above. Moreover, all prompt emissions scatter and hence, without loss of generality, we also assume that for for each i = 1, . . . , , the terms σ i s π i s are uniformly bounded away from the origin on D × V . We need not assume that the cross sections σ i f π i, j f are uniformly bounded away from the origin for 1 ≤ i, j ≤ , with the exception of i = 1, for which it only makes sense that σ 1 f m j is uniformly bounded away from 0 for each j = + 1, . . . , m. Without loss of generality, we can assume that 0 < λ +1 < · · · < λ m .
We also assume similar boundary conditions to the single-type case in the sense that emissions exiting the physical domain D are killed. That is to say For the second condition, we will write Classical literature suggests that one can integrate delayed neutrons into the setting of the NTE by adding an inhomogeneity corresponding to the integral of incoming delayed neutrons from time −∞ to the present; see e.g. [8]. A vectorial representation such as the one above can be found, however, in the work of [25]. There, only one category of prompt emissions are considered with multiple species of delayed neutrons.
As before, let us define the multi-species backward transport, scattering and fission operators as they appear in MNTE (5.1) and (5.2), with zero action otherwise.
It is not often that MNTE is stated as above in (5.1) and (5.2) in existing literature; see e.g. [25] for presentation of the NTE in a similar vectorial format, which allows for only one category of prompt neutrons.
The requirement that all cross sections are uniformly bounded is by far not the weakest assumption we can make (see e.g. Chapter XXI of [6]).
The precise mathematical sense in which we must understand solutions to the coupled system (5.1) and (5.2) needs some discussion before we can proceed. To this end, we shall first introduce some notational conventions.
As alluded to above, we are interested in an vector space of functions, written as the column vector g(·) = (g(1, ·), . . . , g(m, ·)) T , whose entries g(i, , which is easily verified to be itself an L 2 space with inner product given by Generally speaking, for a scalar quantity which is indexed by i, say a(i), when written without the index, we will understand it to be a column vector. Sometimes we will want to put f ∈ m j=1 L 2 (D × V ) on the diagonal of an m × m matrix, in which case we will write diag( f ). For our transport, scattering and fission operators, we will understand The operator ← S can be handled similarly. We are fundamentally interested in a classical solution to the so-called (initial-value) abstract Cauchy problem (ACP) where u t is treated as a column vector belonging to the space m j=1 L 2 (D × V ), for t ≥ 0. Specifically this means that, (u t , t ≥ 0) is continuously differentiable in this space. In other words, there exists aψ t ∈ m The theory of c 0 -semigroups gives us a straightforward approach to describing the unique solution to (5.5). Recall that a c 0 -semigroup also goes by the name of a strongly continuous semigroup and, in the present context, this means a family of time-indexed operators, To see how c 0 -semigroups relate to (5.5), let us define A and standard theory (cf. [11]) tells us that V t [g] ∈ Dom( ← A) for all t ≥ 0, with g as above. Proposition II.6.2 of [11] now gives us the relevance to (5.5). The reader may well have wondered where the second boundary condition in (5.3) has gone in the above formulation. This is a matter of interpretation of ( ← T, Dom( ← T)), and hence the generator ( ← A, Dom( ← A)), as we now discuss. We are interested in the advection semigroup with exponential killing and killing on the boundary of D, where κ D r ,υ := inf{t > 0 : r + υt / ∈ D}.
It is a straightforward exercise, see e.g. Theorem 2 in Chapter XXI of [6], to show that U is a c 0 -semigroup with generator ← T. Its domain satisfies Here, by υ · ∇g ∈ L 2 (D × V ) we mean that υ · ∇g exists in the distributional sense and is integrable in the space  T). To see why, we need only consider that the linear operators of the form where α and π i, j are nonnegative, measurable and uniformly bounded. The proof is a straightforward exercise which uses the Cauchy-Schwarz inequality; see for example Lemma XXI.1 of [6]. It follows that Dom( Note there is no particular necessity to put solutions in an L 2 space, one might equally work with the space m i=1 L p (D×V ), for p ∈ (1, ∞). As the reader might suspect, solutions of the backwards equation in an L p space comes hand in hand with a similarly formulated solution to the forward equation in the conjugate space m i=1 L q (D × V ), where q −1 + p −1 = 1. See for example Chapter XXI of [6] or [25]. The reader will note the exclusion of the L 1 and L ∞ conjugacy. The reason for the exclusion boils down to the cumbersome nature of the advection operator ← T = υ · ∇. Quite simply it is not possible to verify the strong continuity property of the advection semigroup where κ D r ,υ := inf{t > 0 : r + υt / ∈ D}. Hence we cannot give a meaning to υ · ∇ as a c 0 -semigroup on L ∞ (D × V ). This is unfortunate as the latter is the more natural setting for probabilistic interpretation of solutions to the ACP. Having said that, the backwards scattering and fission operators, respectively One of our main results will be to establish the asymptotic (2.8) but now in the current setting. Recall that we have assumed that D ⊆ R 3 is a smooth open pathwise connected bounded domain of concern such that ∂ D has zero Lebesgue measure. and To give a precise value for ε, suppose we enumerate the eigenvalues of ← A in decreasing order by the set {λ (1) , . . . , λ (n) } (noting from earlier that we have at least λ (1) = λ c ). Then λ (n) > −λ +1 and we can take any ε such that ε < λ c − (λ (2)

Remark 5.1
It could be argued that the assumptions in the above theorem rule out the possibility that we may, for example, include alpha or beta emissions emissions in the model for that particular conclusion. Whilst alpha and beta emissions may scatter, they are not energetic enough to cause fission. The irreducibility conditions (5.11) and (5.12) would thus fail. On the other hand, it is also known that when such particles are energetic enough, they can draw gamma radiation or positrons out of nuclei when passing in close proximity. If the latter are sufficiently energetic, then they can induce fission.

Multi-species Neutron Branching Process
Heuristically speaking, (5.5) can be thought of as being closely related to the expectation semigroup of a Markov branching process, or Multi-species nuclear branching process (MNBP) as we shall call it, whose infinitesimal generator is In order to describe the system as Markovian, we will represent it by the atomic measures where A is a Borel subset of D × V and δ is the Dirac measure defined on the same space. Then the system can be described via the m-tuple X t (·) = (X t (1, ·), . . . , X t (m, ·)), t ≥ 0, which evolves as follows.
A emission of type i ∈ {1, . . . , } with configuration (r , υ) moves in a straight line with velocity υ from the point r until one of the following events occur: • The emission leaves the domain, at which point it is killed.
• Independently of all other emissions, a scattering event occurs when a emission comes in close proximity to an atomic nucleus and, accordingly, makes an instantaneous change of velocity. For an emission in the system of type i ∈ {1, . . . , } with initial position and velocity (r , υ), if we write T i s for the random time until the next scattering occurs, then, independently of any other physical event that may affect the emission, • When scattering of an emission of type i ∈ {1, . . . , } occurs at space-velocity (r , υ), the new velocity is selected independently with probability π i s (r , υ, υ )dυ . • Independently of all other emissions, a fission event occurs when an emission smashes into an atomic nucleus. For an emission in the system with initial position and velocity (r , υ), we will write T i f for the random time that the next fission occurs. Then independently of any other physical event that may affect the emission, • When fission occurs, the smashing of the atomic nucleus releases a random number of other prompt emissions of type i = 1, . . . , , say N i, j ≥ 0, which are ejected from the point of impact with randomly distributed, and possibly corollated, velocities, say {υ i, j k : k = 1, . . . , N i, j }. When fission occurs at location r ∈ D from a emission with incoming velocity υ ∈ V , the quantity π i, j f (r , υ, υ )dυ describes the average number of type j prompt emissions released from nuclear fission with outgoing velocity in the infinitesimal neighbourhood of υ . In particular • Note, the possibility that Pr(N i, j = 0) > 0 is possible. If i = j = 1 then this is tantamount to neutron capture or further decomposition into subatomic particles which are not counted. • Further, if the initial emission is a (type 1) neutron, a fission event (occurring at rate σ 1 f ) may result in the production of unstable isotopes (which later release delayed emissions). On this event, an average number, m j (r , υ), of type j ∈ { + 1, . . . , m} isotopes will be produced from a collision at position r from a neutron with incoming velocity υ. The isotopes will inherit the configuration of the incoming nucleus at the time of collision.
An isotope of type i ∈ { + 1, . . . , m} with inherited physical configuration (r , υ) stays in the same place for an exponentially distributed amount of time with rate λ i . At this point, it produces a random number of type j ∈ {1, . . . , } prompt emissions, the average number of which, along with their corresponding velocities, are chosen according to π i, j f (r , υ, υ ), in a similar way to previously described. We note that although unstable isotopes stay in the same spatial position, we will still assign them a velocity as a 'mark'.
In all cases, it is a natural make the following physical assumption which will remain in force throughout. For non-negative and uniformly bounded g where P δ (i,r,v) is law of the process started from a single type i emission with configuration (r , υ) with corresponding expectation operator E δ (i,r,v) . As we have assumed that all cross sections are uniformly bounded, ignoring spatial trajectories of neutrons (in particular those that are killed by leaving the domain D), it is straightforward to compare the growth of (ψ t [g], t ≥ 0) against that of a continuous-time Galton-Watson process with growth rate η{(m × n max ) − 1}, where η = sup 1≤i≤ ,r ∈D,υ∈V σ i f (r , υ) + max +1≤i≤m λ i . The rate of growth η{(m × n max ) − 1} simply assumes that each emission of type i gives rise to at most n max emissions of any other type and at a rate which is uniformly bounded by a uniform upper bound of all possible rates at which fission events occur. Note this rate takes account of the emission count introduced into the system at a fission event and the single emission removed from the system which caused the fission event.
It is also straightforward to stochastically upper bound the process 1, X t , t ≥ 0, by the aforesaid continuous-time Galton Watson process on the same probability space. The latter process branches whenever X does, topping up the number of offspring always to n max , but also it has additional independent branching events at rate (η − 1 (1≤i≤ ) σ i f (r , υ) − 1 ( +1≤i≤m) λ i ) always producing precisely n max offspring of each of the m possible emissions.
If we denote this Galton-Waton process by (Z t , t ≥ 0), then we have both the stochastic bound 1, X t ≤ Z t ≤ Z t+s , for all s, t ≥ 0 and the upper estimate sup 1≤i≤m,r ∈D,υ∈V If we put g in the smaller space m i=1 C + (D × V ), the space of non-negative, continuous and uniformly bounded vector functions on (D × V ), then we also have by a dominated convergence argument, lim t→0 ψ t [g] = g in the pointwise sense. Otherwise the latter convergence is not necessarily clear.
The name 'expectation semigroup' is earned thanks to the behaviour of (ψ t , t ≥ 0) under an application of the Markov branching property. Indeed, associated to the MNBP are the probabilities P μ for atomic measures of the form (1,r i,1 ,υ i,1 ) , . . . , ,r i,m ,υ i,m ) . =: (μ 1 , . . . , μ m ). (6.5) The Markov branching property dictates that, for g ∈ m i=1 L 2 (D × V ) as before and t ≥ 0, Here we are abusing our earlier notation in (5.4) and writing for finite atomic measures μ of the form (6.5), Hence, by conditioning on the configuration of the system at time t ≥ 0, we have, for s ≥ 0, The expectation semigroup property of (ψ t , t ≥ 0) does not imply that it is necessarily a c 0 -semigroup on m i=1 L 2 (D × V ). Recalling our earlier discussion, if we were able to work with (5.5) in the setting of a c 0 -semigroup on m i=1 L ∞ (D × V )), then we would be much closer to being able to match the expectation semigroup (ψ t , t ≥ 0) to the solution (u t , t ≥ 0). But even then, problems would occur with verifying strong continuity at the origin.
Nonetheless, classical literature supports the view that it is the physical processes, i.e. in this setting the MNBP, that provides a stochastic representation of the solution to the backward MNTE. The authors are not aware of a formal proof of this fact. We will nonetheless try to address this point shortly in Sect. 8. In the mean time, let us present an alternative 'mild' form of the MNTE (also called a Duhamel solution in the PDE literature) which the semigroup (ψ t , t ≥ 0) more comfortably solves.

Lemma 6.1 The expectation semigroup
Before proceeding to the proof, let us remark that, in the statement of the theorem, we are not working with , but a pointwise shift operator. The reader will recall from the discussion preceding (5.10) that Proof of Lemma 6.1 First suppose we start with an emission of type i. By splitting the expectation in the definition of ψ t [g] at the first scattering or fission event, and remembering that the time κ D r ,υ defined in (5.8) is deterministic, we have for r ∈ D and υ ∈ V , Now appealing to an analogue of Lemma 1.2, Chapter 4 in [10] (see also the Appendix of [16]), we can transfer the exponential integrals in each of the terms on the right-hand side above to a potential term in the integral so that we end with which agrees with (6.8), for 1 ≤ i ≤ . Following a similar approach, for + 1 ≤ i ≤ m, r ∈ D, υ, ∈ V , we also get (6.10) noting in particular that, for + 1 ≤ i ≤ m, ← S i ≡ 0. Now putting (6.9) and (6.10) together we obtain (6.8).
For uniqueness, suppose that (ψ (i) t , t ≥ 0), i = 1, 2 are two bounded solutions to (6.8). (2) t [g]| and note that, for i = 1, . . . , m, (6.11) for some constants C 1 , C 2 ∈ (0, ∞), where the final inequality follows on account of all cross sections being uniformly bounded. Now defineχ t [g] = sup 1≤i≤m,r ∈D,υ∈V χ t [g](i, r , υ), t ≥ 0. From (6.11) we have that Reversing the order of integration on the right-hand side above and then applying Grönwall's Lemma allows us to conclude that χ t [g] ≡ 0, which shows uniqueness. To describe this more precisely, we need to introduce the notion of a multi-species neutron random walk (MNRW). In the current setting this means a continuous-time typed random walk by (J t , R t , ϒ t ), t ≥ 0, on {1, . . . , m}×(D × V ) with additional cemetery state { †} when it exits the physical domain D or an emission otherwise disappears from the system. The MNRW is described by two fundamental quantities (which are functions of the current particle type, spatial position and velocity). First, a scattering rate α i (r , υ), i ∈ {1, . . . , m}, r ∈ D, υ, υ ∈ V , such that α i (r , υ) = λ i , for i ∈ { +1, . . . , m}. Second, a scattering probability kernel π i, j (r , υ, υ ), i, j ∈ {1, . . . , m}, r ∈ D, υ, υ ∈ V . In the spirit of the description of the MNBP, the MNRW is described as follows.

Multi-species Neutron Random Walk and the Many-to-one Lemma
When the MNRW is of type i ∈ {1, . . . , } with configuration (r , υ), it moves in a straight line with velocity υ from the point r until one of the following events occur: • When the MNRW position moves out of D or e.g. it decomposes into an emission type that is not counted, or is captured in a nucleus, it is instantaneously killed. • A scattering event occurs and, accordingly, the MNRW keeps the same emission type but makes an instantaneous change of velocity. If we write T i s for the random time until the next scattering occurs, then, • When scattering of an emission of type i ∈ {1, . . . , } occurs at space-velocity (r , υ), the new velocity is selected independently with probability π i (r , υ, υ )dυ .
Otherwise, if + 1 ≤ i ≤ m, then the emission remains motionless, i.e. the random walk is dormant, holding its initial position r , but retaining the velocity υ as a mark. After an independent and exponentially distributed random time with rate λ i , the particle transfers it type j ∈ {1, . . . , } and acquires a new velocity υ with probability density π i (r , υ, υ ).
We can associate to the MNRW the infinitesimal generator We thus refer to the process as an ← L-MNRW. With the notion of the MNRW in hand, let us consider the following algebraic manipulations . For i ∈ {1, . . . , }, j ∈ {1, . . . , m}, (r , υ) Note, in particular, that for each fixed 1 ≤ i ≤ m, r ∈ D and υ ∈ V , π i, j (r , υ, υ ) is a probability distribution on {1, . . . , m} × V in the sense that m j=1 V π i, j (r , υ, υ )dυ = 1. Note also that the assumption With simple algebra, we may now identify where, for f ∈ Dom( ← A) (for which it was remarked earlier that it is equal to Dom( ← T)), and ← L is given by (7.2). Heuristically speaking, we have algebraically gathered all of the operators into the infinitesimal generator of an ← L-MNRW and local potential β. This has the attraction of leading us the aforementioned single emission representation of the solution to the MNTE using a single-emission Feynman-Kac representation. Said another way, this means that one would expect that, in the appropriate sense, the solution to the NTE to be represented in the form Here P (i,r ,v) for the law of the ← L-MNRW starting from a single emission with configuration (i, r , υ), and E (i,r ,v) for the corresponding expectation operator.
Appealing to the Markov property for (J , R, ϒ), it is not difficult to show that a semigroup property similar to (6.7) holds. That is to say, for s, , r , υ).
Similarly to the case of (ψ t [g], t ≥ 0), if we put g in the smaller space m i=1 C + (D × V )) then we also have lim t→0 φ t [g] = g in the pointwise sense, but otherwise strong continuity at t = 0 is unclear. Note also that, since all cross sections are uniformly bounded, then so is β (in all of its variables) by a constant, sayβ. Hence, for As with the case of (ψ t [g], t ≥ 0), the notion that (φ t [g], t ≥ 0), solves (5.5) is not a straightforward claim. Nonetheless, as one might expect, these two expectation semigroups are equal and, we can see this by relating back to (6.8).
Indeed, by conditioning the expectation in the definition of φ t [g] on the first scattering event, and then appealing to the Lemma 1.2, Chapter 4 in [10] in a similar manner to what was done in the proof of Lemma 6.1, one easily deduces the below result. In the the spatial branching process literature, this would be called a 'many-to-one' lemma.

Consolidating the ACP with the Expectation Semigroup
We want to understand how the m i=1 L 2 (D × V ) semigroup (V t , t ≥ 0) that represents the unique solution to the Abstract Cauchy Problem (5.5) relates to the expectation semigroups (ψ t , t ≥ 0) and (φ t , t ≥ 0) that offer two different stochastic representations to the mild Eq. (6.8).
We start by noting that if g ∈ m i=1 L + ∞ (D × V ), then, on account of the fact that , it makes makes sense to consider the comparison with (V t [g], t ≥ 0) (defined in (5.6)) for the more restrictive choice g ∈ m i=1 L ∞ (D × V ). The natural setting in which to make the comparison is in the space m i=1 L 2 (D × V ) as, by (6.4), ψ t [g] ∞ < ∞ and the latter implies ψ t [g] 2 < ∞, again thanks to the fact that Vol Before moving to its proof, the reader should take care to note that this does not imply that (V t , t ≥ 0) and (ψ t , t ≥ 0) agree as c 0 -semigroups on m i=1 L 2 (D × V ). In particular, the comparison between the two semigroup operators is only made on m i=1 L 2 (D × V ), and (ψ t , t ≥ 0) was not (and in fact cannot be) shown to demonstrate the strong continuity property on m i=1 L 2 (D × V ). Remark 8.1 If we consider Theorem 8.1 in light of Theorem 5.3, noting that (ψ t [g], t ≥ 0) is a uniformly bounded sequence, it is tempting to want to say that the leading eigenfunction ϕ belongs to m i=1 L ∞ (D × V ). This is not the case necessarily and remains to be proved. In the setting of a single type of emission, this will be demonstrated in the forthcoming paper [14].

Proof of Theorem 8.1 Consider the adjusted ACP with inhomogeneity given by
By taking the difference of two solutions and invoking the uniqueness of the ACP in m i=1 L 2 (D × V ) with initial data g = 0, we note that the solution to (8.1) is unique in However, on the one hand, it is straightforward to verify that On the other hand, taking account of the fact that (V t [g], t ≥ 0) solves (5.5), it is also the case that

solves (8.1). Uniqueness thus tells us that on
where in the second equality we have reversed the direction of integration. In conclusion, where as (ψ t [g], t ≥ 0) solves (6.8) in the pointwise sense, On the other hand, we know that To this end, let us note that, for T > 0, and where in the first inequality we have used Jensen's inequality and Cauchy-Schwarz in the second. Moreover, where the inequality follows as a consequence that, for each υ, the integral D 1 (s<κ D r,υ ) v(i, r + υs, υ) 2 dr integrates over a subdomain of D. Also, we have for the operator ← S (and similarly for where the constant C comes from the fact that σ is uniformly bounded. The final inequality in (8.5) together with Grönwall's Lemma now tells us that ω t 2 = 0, for all t ≤ T . Since T is chosen arbitrarily, it follows that (ψ t [g], t ≥ 0) and The conclusion of this section is that it is not unreasonable to now understand the expectation semigroups (ψ t [g], t ≥ 0) and (φ t [g], t ≥ 0) for non-negative, bounded and measurable g on D × V as the 'solution' to the MNTE in place of (V t [g], t ≥ 0) for the same class of g. Indeed, the two agree both in m i=1 L 2 (D × V ) and hence (dr × dυ)-Lebesgue almost everywhere.
The reader will also note that from the perspective of Monte Carlo simulation, the expectation semigroup (φ t [g], t ≥ 0) carries the potential to be exploited in a way that (ψ t [g], t ≥ 0) cannot. More precisely, where branching trees are difficult to simulate and are not convenient for Monte Carlo computational parallelisation, random walks are. This simple idea is explored in greater detail in the accompanying paper to this one [5].

Asymptotic Behaviour of the MNTE: Proof of Theorem 5.3
In this section we return to the fundamental notion that the solution to the MNTE in the form (5.5) is described by its leading asymptotics for large times. That is to say, we give the proof of Theorem 5.3. Our proof follows closely ideas found in Chapters 4 and 5 of [25].
Recall that the quantities α i , π i, j , β i , i, j = 1, . . . , m were defined in (7.3), (7.4) and (7.5) respectively. They were arranged into the operator These are integral operators, which take the form A similar computation to (8.4) also shows that K i, j g ∈ L 2 (D × V ) when g ∈ L 2 (D × V ). Then from (5.1) and (7.4), taking care to note the use of the indicators for the inclusion of terms for different indices, we can write, for 1 ≤ i ≤ , for g ∈ Dom( ← A), With this notation, write Then the abstract Cauchy problem (5.5) on m j=1 L 2 (D × V ) may now be written in matrix form ∂ ∂t u t = Au t , t ≥ 0.
where A = T + K and which generates the strongly continuous semigroup (U T t , t ≥ 0) given by . In order to prove Theorem 5.3, we consider a different operator that is related to A as follows. Consider the eigenvalue problem (ϕ(1, ·), . . . , ϕ( , ·)) and ϕ • (·) = (ϕ ( + 1, ·), . . . , ϕ(m, ·)) so that ϕ is the concatenation (ϕ • , ϕ • ). Separating this into prompt and delayed initial emissions, it can be written as where I m− is the (m − ) × (m − ) identity matrix. Substituting the second equation into the first, we get and Our strategy is to show that there exists a λ c such that (λI − T) −1 K • (λ) has a leading eigenvalue 1, and that this is equivalent to λ c being an eigenvalue of A. The tool we shall use to do this is the Krein-Rutman theorem, which we recall here for convenience in a format that is appropriate for our use; c.f. [6, p. 286]. Our proof of Theorem 5.3 requires the following intermediary result below. Before stating it, the reader is reminded that the eigenvalues λ +1 , . . . , λ m are arranged so that λ +1 is the smallest. Thus, the condition λ > −λ +1 ensures that K • (λ) is well defined. In particular, (λI m− + ) is invertible. We will use the obvious meaning for I .

Proposition 9.1 Under the assumptions of Theorem
Proof In relation to the Krein-Rutman theorem stated above, our Banach space is X = m i=1 L 2 (D × V ) and the corresponding cone is It is clear that this cone is convex, and since every L 2 function can be written as the difference of its positive and negative parts, C satisfies the assumptions of the theorem. We now break the rest of the proof into a number of steps which are stated with a proof immediately afterwards.
Step 1 First we claim that (λI − T) −1 K • (λ) is a compact operator. Fix 1 ≤ i, j ≤ m. By Fubini's Theorem we have that r → K i, j f (r , υ) is measurable for g ∈ L 2 (D × V ). The operators K i, j are also integral operators and therefore are continuous on L 2 (V ) and compact. The assumed piecewise continuity of the cross sections σ i s π i s and σ i f π i, j f and the boundedness of the domain V is sufficient to ensure that r → K i, j · (r , ·) is continuous under the operator norm on L 2 (V ) and hence {K i, j · (r , ·) : r ∈ D} forms a relatively compact set in the space of linear operators on L 2 (V ). With these properties, the mapping r → K i, j · (r , ·), for r ∈ D, is said to be regular. One similarly (but more easily) shows that r → M i, j · (r , ·) is regular for r ∈ D as operators on L 2 (V ). By linearity, this implies that, for 1 ≤ i, j ≤ , the mapping r → K • (λ) i, j is also regular. Hence, by [25,Theorem 4 We use de Pagter's Theorem, cf. [25,Theorem 5.7], which says that the spectral radius of an irreducible operator is strictly positive; that is to say r (λI − T) −1 K • (λ) > 0. In turn the Krein-Rutman theorem 9.1 states that r (λI − T) −1 K • (λ) is thus an eigenvalue for the operator (λI − T) −1 K • (λ) with a corresponding non-negative eigenfunction ϕ • λ .
Proof of Theorem 5.3 (i) In looking for a non-negative eigenfunction of ← A with real eigenvalue, our earlier discussion tells us we must equivalently look for a solution to (9.5) and hence (9.7). This is equivalent to finding a real value λ c such that r (λ c I −T) −1 K • (λ c ) = 1. We again achieve this goal in steps.
Step 3 In this penultimate step, we show that we have found a non-negative function of ← A, with eigenvalue λ c . We have the existence of a λ c > −λ +1 such that r ((λI − T) −1 K • (λ)) = 1. That is to say, thanks to Proposition 9.1, we have found ϕ • = ϕ • λ c which solves (9.7), which in turn, thanks to (9.6) gives us that ϕ • = (λI m− + ) −1 K • ϕ • λ c so that with the concatenation we have the eigensolution which is equivalent to ← Aϕ = λ c ϕ.
Step 4 For the final step we need to show that λ c is the leading real eigenvalue of  where σ ( A) is the spectrum of the operator A or equivalently of ← A. Moreover we need to show that it is simple and isolated.
We first note that since we have shown that λ c ∈ σ ( A), in particular that the spectrum is non-empty, it follows from [ Note the Theorem from which the above proposition is derived in [1,p. 359,Theorem 22] requires as a sufficient condition that (λI − T ) −1 K is compact, where I is an m × m identity matrix. This fact easily follows from the conclusion in Step 1 of the proof of Proposition 9.1.
Finally we can complete the proof of Theorem 5.3.
Suppose we enumerate the eigenvalues in σ ( A) ∩ {λ : Re(λ) > −λ +1 } in decreasing order by the set {λ (1) , . . . , λ (n) } (noting from earlier that we have at least λ (1) = λ c and λ (n) > −λ +1 ). Then, from [6, p. 265], for g ∈ Dom( ← A), we have as t → ∞, where k are projectors in Dom( ← A). We are really only interested in the projection onto the eigenfunction that we know exists in the real part of the spectrum. The projector 1 can be written in the form whereφ is the left-eigenfunction with eigenvalue λ c , which is guaranteed to exist by examining the preceding arguments for Note that since, according to Proposition 9.2, λ c is isolated, there exists a ε > 0 such that λ (2) ∨ (−λ +1 ) < λ c − ε, where we understand λ (2) = −∞ if n = 1. The statement of part (ii) of Theorem 5.3 now follows.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Glossary
For convenience, at the request of the referee, we include a short glossary of the shorthand terminology.