Homogenization in the energy variable for a neutron transport model

Abstract

This article addresses the homogenization of linear Boltzmann equation when the optical parameters are highly heterogeneous in the energy variable. We employ the method of two-scale convergence to arrive at the homogenization result. In doing so, we show the induction of a memory effect in the homogenization limit. We also provide some numerical experiments to illustrate our results. In the Appendix, we treat a related case of the harmonic oscillator.

Keywords

Linear Boltzmann equation self-shielding two-scale convergence memory effects

1. Introduction

Neutron transport plays a key role in a number of scientific and engineering areas. It is, for example, relevant in understanding certain atmospheric processes and it also plays a major role in the field of nuclear engineering – the safety of nuclear reactors and shielding. The evolution of neutrons in all these problems can be modelled by a linear Boltzmann equation whose coefficients, usually called optical parameters or microscopic cross sections, present rapid oscillations related to the neutron energy. In this work, we derive two-scale homogenization results related to this phenomenon which reveal the presence of a memory effect in the limit equation.

In what follows, Ω is a bounded domain of $R^{d}$ with $C^{1}$ boundary and $Y : = {(0, 1)}^{d}$ is the unit cube in $R^{d}$ . Denoting by $f = f (t, x, v)$ the population density of neutrons which are located at $x \in Ω$ (position) at time $t \in R_{+}$ and travelling with velocity $v \in V \subset R^{d}$ , the neutron gas can be described by the linear Boltzmann equation (sometimes called the neutron transport equation) $\begin{matrix} (1) & \partial_{t} f + v \cdot \nabla_{x} f + σ (x, v) f - \int_{V} κ (x, v \cdot v^{'}) f (t, x, v^{'}) d v^{'} = 0, \end{matrix}$ where the non-negative functions σ and κ denote the total cross-section of the background material and the scattering kernel respectively. In the following, we will sometimes refer to the pair $(σ, κ)$ as the optical parameters. Note that the velocity variable v can be expressed via the couple $(ω, E)$ , where $ω = v / | v |$ is the trajectory angle of the neutron and $E = m | v |^{2} / 2$ is the kinetic energy, m being the mass of the neutron. Assuming that E ranges in $[E_{min}, E_{max}]$ , the same equation can equivalently be written in terms of the neutron flux $\begin{matrix} φ (t, x, ω, E) = φ (t, x, v) : = | v | f (t, x, v), \end{matrix}$ which satisfies the equation $\begin{array}{l} \sqrt{\frac{m}{2 E}} \partial_{t} φ + ω \cdot \nabla_{x} φ + σ (x, ω, E) φ \\ (2) & - \int_{E_{\min}}^{E_{\max}} \int_{| ω^{'} | = 1} κ (x, ω \cdot ω^{'}, E, E^{'}) φ (x, ω^{'}, E^{'}) d ω^{'} d E^{'} = 0, \end{array}$ where the optical parameters of the linear Boltzmann equation are appropriately redefined. The above evolution equation shall be supplemented by suitable initial data, i.e., $φ (0, x, ω, E) = φ_{in} (x, ω, E)$ , and boundary data on the incoming phase-space boundary. For simplicity, we shall consider absorption-type boundary data, i.e., $\begin{array}{l} φ (t, x, ω, E) = 0 \forall t > 0 and for (x, ω) \in Γ_{-} : = {(x, ω) \in \partial Ω \times S^{d - 1} : n (x) \cdot ω < 0}, \end{array}$ where $n (x)$ denotes the unit exterior normal to the boundary $\partial Ω$ at the point x.

The numerical computation of solutions to (2) is challenging for the following reasons:

the problem is high-dimensional: its solution depends on $(2 d + 1)$ variables and the integral term couples the angular and energy variables and gives rise to dense matrices,

the solutions have low regularity,

the optical parameters can present high oscillations in space (due to the spatial heterogeneity of the materials) and in energy.

The first two items have been the main focus of numerous works on the numerical analysis of equation (1) and (2). We refer to [1,11,13] for general references on the topic and to a recent work [5] where modern operator compression and adaptive techniques have been applied to efficiently attack the high-dimensional aspects of the problem.

The last item has motivated a considerably huge amount of literature in the theory of homogenization (see [7] and references therein). However, to the best of our knowledge, the existing mathematical theory addresses only high oscillations in the spatial variable and no rigorous results seem to address high oscillations in the energy variable. This point has been treated thus far only in the engineering community where the problem is known as energy self-shielding or resonant absorption, which will be further commented on towards the end of this introduction. In this context, the present contribution is to bridge the gap between theory and practice by arriving at some rigorous homogenization results. This involves certain modelling assumptions on the multiscale behavior of the optical parameters which are in agreement with physical observations. Real experiments reveal strong oscillations in σ as a function of E when the neutrons interact with relevant materials like, for example, Uranium 238 (see Fig. 1). A similar behaviour is also observed for the scattering kernel κ. These facts motivate us to study the multi-scale linear Boltzmann equation $\begin{array}{l} \sqrt{\frac{m}{2 E}} \partial_{t} φ^{ε} + ω \cdot \nabla_{x} φ^{ε} + σ^{ε} (x, ω, E) φ^{ε} \\ (3) & - \int_{E_{\min}}^{E_{\max}} \int_{| ω^{'} | = 1} κ^{ε} (x, ω \cdot ω^{'}, E, E^{'}) φ^{ε} (x, ω^{'}, E^{'}) d ω^{'} d E^{'} = 0, \end{array}$ where $0 < ε ≪ 1$ is a small parameter and $\begin{array}{l} (4) & σ^{ε} (x, ω, E) = σ (x, ω, E, \frac{E}{ε}); κ^{ε} (x, ω \cdot ω^{'}, E, E^{'}) = κ (x, ω \cdot ω^{'}, E, E^{'}, \frac{E^{'}}{ε}), \end{array}$ where $σ (x, ω, E, y)$ and $κ (x, ω \cdot ω^{'}, E, E^{'}, y^{'})$ are both assumed to be periodic in the y and $y^{'}$ variables respectively. The equation is complemented with zero incoming flux condition on the phase-space boundary and an initial condition $φ_{in}$ which we assume to be in $L^{2} (Ω \times S^{d - 1} \times (E_{min}, E_{max}))$ .

In addition to the above hypotheses, we also assume that there exists $α > 0$ such that for all $ε > 0$ , $\begin{array}{l} (5) & σ^{ε} (x, ω, E) - {\bar{κ}}^{ε} (x, ω, E) ⩾ α and σ^{ε} (x, ω, E) - {\tilde{κ}}^{ε} (x, ω, E) ⩾ α, \end{array}$ where $\begin{matrix} (6) & \{\begin{matrix} {\bar{κ}}^{ε} (x, ω, E) : = \int_{E_{min}}^{E_{max}} \int_{S^{d - 1}} κ^{ε} (x, ω \cdot ω^{'}, E, E^{'}) d ω^{'} d E^{'} \\ {\tilde{κ}}^{ε} (x, ω, E) : = \int_{E_{min}}^{E_{max}} \int_{S^{d - 1}} κ^{ε} (x, ω \cdot ω^{'}, E^{'}, E) d ω^{'} d E^{'} . \end{matrix} \end{matrix}$ From a physical point of view, these assumptions mean that we place ourselves in the so-called subcritical regime where absorption phenomena dominate scattering.

Fig. 1.

Total cross-section σ of Uranium 238 as a function of the energy according to the JEFF 3.1 library [12]. Note the highly oscillatory interval for $E \in [1 eV; 10^{4} eV]$ .

From the homogenization viewpoint, transport dominated equations such as (3) are particularly challenging since the structure of the partial differential equation becomes more complex after taking the homogenization limit. This is due to the memory effects induced in the limit that make the dynamics be no longer defined by a semigroup [3,15,20]. This, in turn, entails difficulties in the numerical solution of the homogenized equation since the memory effects dramatically increase the computational complexity in terms of the number of degrees of freedom to be used in order to retrieve a certain target accuracy.

Our main contribution is the homogenization result given in Section 3, where we derive a homogenized equation for the neutron transport problem when the optical parameters oscillate periodically in the energy variable. The result is derived employing the theory of two-scale convergence. For technical reasons that we will explain later on in more detail, it is more convenient to work with scattering kernels $κ^{ε}$ exhibiting separation in the E and $E^{'}$ variables as follows: $\begin{matrix} κ^{ε} (x, ω \cdot ω^{'}, E, E^{'}) : = κ_{1} (x, ω \cdot ω^{'}, E) κ_{2} (x, ω \cdot ω^{'}, E^{'}, \frac{E^{'}}{ε}) \end{matrix}$ with $κ_{2} (x, ω \cdot ω^{'}, E^{'}, y^{'})$ being periodic in the $y^{'}$ variable. A result without this assumptions on $κ^{ε}$ requires further investigation as it is not apparent whether one can derive closed form homogenized equations in the latter case.

The homogenized equation is integro-differential and presents a memory term. To derive the results, we base our strategy on the method of characteristics, and hence we first derive a homogenization result for the associated ordinary differential equation (Section 2). In there, we also show that this result is in agreement with previous works on memory effects by Tartar [19], [20, chapter 35]. An interesting result in its own right is that our technique gives an explicit expression of the memory kernel that, in the situation studied by Tartar, is equal to the implicit expression given in [19], [20, chapter 35]. Finally we consider in Section 4 some numerical experiments on a simpler model, mainly focused on the convergence rate towards the homogenized solution as the period of the oscillating terms tends to zero (i.e., rapid oscillations). In the Appendix, we treat the related case of a harmonic oscillator.

We conclude this introduction by comparing our approach to some standard methods from the nuclear engineering community for treating self-shielding phenomena. To the best of the authors’ knowledge, the most widespread technique is a two-stage method originially proposed in [14] by M. Livolant and F. Jeanpierre (we refer to [17, Chapters 8 and 15] for an introductory overview). It consists in finding first the averaged optical parameters which are then plugged into a multigroup version of equation (1) to compute the behavior of the flux on large and geometrically complex domains such as nuclear reactors. The pre-computation of the averaged parameters is done on a cell problem involving a much simpler spatial geometry and simplified physics. It is nevertheless carefully designed with elaborate physical considerations in a way to keep as much consistency as possible with respect to the original problem. Note that, while this approach implicitly assumes that the homogenized equation is of the same nature as the original Boltzmann problem (1), our starting point is fundamentally different in the sense that we do not postulate any final form of the limit equation. Our goal is precisely to discover its form from the only assumption that the optical parameters oscillate in energy. As a result, our methodology and conclusions are different from the ones discussed in [14] and do not involve a pre-computation on a cell problem. Another approach, also based on averaging the optical parameters, is the so-called multi-band method (see [4]), where, like in the previous method, the limit equation is assumed to be a Boltzmann equation. Finally, a more recent approach based on averaging arguments taken from results of homogenization of pure transport equations has recently been proposed in [8]. The initial problem there is a time-independent Boltzmann source problem with no oscillations in the scattering kernel.

2. Homogenization of an ordinary differential equation

It turns out that the homogenization of the multi-scale linear Boltzmann equation (3) is closely related to the homogenization of an evolution equation studied by Luc Tartar in the 1980’s. Hence, we begin this section by recalling his example and by giving an alternate proof for the phenomenon of memory effect by homogenization, as was demonstrated by him in [19]. For the unknown $u^{ε}$ , consider the differential equation $\begin{matrix} (7) & \partial_{t} u^{ε} + σ (\frac{x}{ε}) u^{ε} = 0; u^{ε} (0, x) = u_{in} (x) . \end{matrix}$ Even though Tartar treats a more general setting, Theorem 1 recalls the result from [19] adapted to when σ has only periodic oscillations. In the sequel, we use the following notation for the Laplace transform (in the time variable) of a function: $\begin{matrix} \hat{f} (p) : = \int_{0}^{\infty} e^{- p s} f (s) d s for p > 0 . \end{matrix}$ In the section, for any $v \in L^{1} (Y)$ , $\begin{matrix} ⟨ v ⟩ : = \int_{Y} v (y) d y \end{matrix}$ denotes the average of v in Y.

Theorem 1 (Tartar, [19]).

Let the coefficient $σ (\cdot)$ in ( 7 ) be a strictly positive, bounded and purely periodic coefficient of period Y. Then the $L^{\infty}$ weak ∗ limit $u_{hom} (t, x)$ of the solution family $u^{ε}$ satisfies the following integro-differential equation $\begin{matrix} (8) & \{\begin{matrix} \partial_{t} u_{hom} (t, x) + ⟨ σ ⟩ u_{hom} (t, x) - \int_{0}^{t} M (t - s) u_{hom} (s, x) d s = 0 \\ u_{hom} (0, x) = u_{in} (x) \end{matrix} \end{matrix}$ where the memory kernel $M (τ)$ is given in terms of its Laplace transform as follows $\begin{array}{l} (9) & \hat{M} (p) = p + ⟨ σ ⟩ - B (p) = \int_{Y} (p + σ (y) - B (p)) d y \forall p > 0, \end{array}$ with the constant $B (p)$ taking the value $\begin{array}{l} B (p) : = {(\int_{Y} \frac{d y}{p + σ (y)})}^{- 1} . \end{array}$

Note that the memory kernel is expressed implicitly by means of the Laplace transform (see formula (9)). This is due to the fact that the proof given in [19] goes via taking the Laplace transform in the t variable. Contrary to this, our alternative approach yields directly an explicit expression without necessitating the use of Laplace transform in the proof. This is because we place ourselves in the periodic homogenization setting, which will in turn aid us in homogenizing the neutron transport equation (3).

Our starting point is the following evolution equation (we consider oscillatory sources as this will resemble more our final targeted equation (3)): $\begin{matrix} (10) & \{\begin{matrix} \partial_{t} u^{ε} (t, x) + σ^{ε} (x) u^{ε} (t, x) = f^{ε} (t, x) & for (t, x) \in (0, T) \times Ω, \\ u^{ε} (0, x) = u_{in}^{ε} (x) & for x \in Ω, \end{matrix} \end{matrix}$ where $T < \infty$ is arbitrary. The coefficients and data in (10) are of the form $\begin{matrix} σ^{ε} (x) : = σ (x, \frac{x}{ε}), f^{ε} (t, x) : = f (t, x, \frac{x}{ε}), u_{in}^{ε} (x) : = u_{in} (x, \frac{x}{ε}), \end{matrix}$ with the locally periodic (in space) coefficient $\begin{matrix} σ (x, y) \in L^{\infty} (Ω; C_{per} (Y)) \end{matrix}$ and data $\begin{matrix} f (t, x, y) \in L^{\infty} ((0, T) \times Ω; C_{per} (Y)), u_{in} (x, y) \in L^{2} (Ω; C_{per} (Y)) . \end{matrix}$ Here we have used the standard notation $C_{per} (Y)$ to denote continuous functions on $R^{d}$ which are Y-periodic. Furthermore we assume that there exists a positive constant $σ_{\min}$ such that $\begin{matrix} σ (x, y) ⩾ σ_{min} \forall (x, y) \in Ω \times Y . \end{matrix}$

Since our analysis relies on the notion of two-scale convergence, before stating our homogenization result, we first recall some basic results and introduce some notation (for proofs and further discussions, readers are encouraged to refer to the works cited here). Two-scale convergence was also used by J.-S. Jiang [9] to treat the case with $u_{in}^{ε} = 0$ and $f^{ε} (x) = f (x)$ for all $ε > 0$ . In that work, the resulting memory kernel was again given implicitly, this time in terms of a certain Volterra equation.

Two-scale convergence was introduced by Gabriel Nguetseng [16] and further developed by Grégoire Allaire [2]. We start by recalling its definition.

Definition 2.
A family of functions $v^{ε} (x) \subset L^{2} (Ω)$ is said to two-scale converge to a limit $v^{0} (x, y) \in L^{2} (Ω \times Y)$ if, for any smooth test function $ψ (x, y)$ , Y-periodic in the y variable, we have $\begin{array}{l} lim_{ε \to 0} \int_{Ω} v^{ε} (x) ψ (x, \frac{x}{ε}) d x = \int_{Ω} \int_{Y} v^{0} (x, y) ψ (x, y) d x d y . \end{array}$

A distinctive result in the two-scale convergence theory is the following result of compactness which says that the above notion of convergence is not void.
Theorem 3 ([2,16]).

Suppose a family $v^{ε} (x) \subset L^{2} (Ω)$ is uniformly bounded, i.e., $\begin{matrix} {‖ v^{ε} ‖}_{L^{2} (Ω)} ⩽ C \end{matrix}$ with constant C being independent of ε. Then, we can extract a sub-sequence (still denoted $v^{ε}$ ) such that $v^{ε}$ two-scale converges to some limit $v^{0} (x, y) \in L^{2} (Ω \times Y)$ .

We additionally have the following useful property of the two-scale limit.

Proposition 4 ([2,16]).

Let $v^{ε}$ be a sequence of functions in $L^{2} (Ω)$ which two-scale converges to a limit $v^{0} \in L^{2} (Ω \times Y)$ . Then $v^{ε} (x)$ converges to $⟨ v ⟩ (x) = \int_{Y} v^{0} (x, y) d y$ weakly in $L^{2} (Ω)$ , i.e., $\begin{matrix} lim_{ε \to 0} \int_{Ω} v^{ε} (x) φ (x) d x = \int_{Ω} φ (x) \int_{Y} v^{0} (x, y) d y d x for all φ \in L^{2} (Ω) . \end{matrix}$

Note that the notion of two-scale convergence is a weak-type convergence as it is given in terms of test functions (see Definition 2).

For any given $g \in L^{\infty} (Y)$ , the linear operator $\begin{matrix} L_{g} v : = g v - ⟨ g v ⟩ \forall v \in L_{per}^{2} (Y) \end{matrix}$ is bounded in $L_{per}^{2} (Y)$ as $\begin{array}{l} ‖ L_{g} h ‖_{L_{per}^{2} (Y)}^{2} = \int_{Y} {| g (y) h (y) - ⟨ g h ⟩ |}^{2} d y = \int_{Y} {| g (y) h (y) |}^{2} d y - {⟨ g h ⟩}^{2} \end{array}$ and we have by Cauchy–Schwarz $\begin{array}{l} | ⟨ g h ⟩ | = | \int_{Y} g (y) h (y) d y | ⩽ {(\int_{Y} {| g (y) h (y) |}^{2} d y)}^{\frac{1}{2}} . \end{array}$ As a consequence, $L_{g} : L_{per}^{2} (Y) \mapsto L_{per}^{2} (Y)$ is the infinitesimal generator of a uniformly continuous semigroup given by $\begin{matrix} e^{t L_{g}} = \sum_{n = 0}^{\infty} \frac{t^{n}}{n!} L_{g}^{n} . \end{matrix}$

We next present the main homogenization result for the evolution (10).

Theorem 5.
Let $u^{ε} (t, x)$ be the solution to the evolution problem ( 10 ). Then $\begin{matrix} u^{ε} ⇀ u_{hom} weakly in L^{2} ((0, T) \times Ω) \end{matrix}$ and $u_{hom} (t, x)$ solves the following integro-differential equation $\begin{matrix} (11) & \{\begin{matrix} \partial_{t} u_{hom} (t, x) + ⟨ σ ⟩ (x) u_{hom} (t, x) - \int_{0}^{t} K (t - s, x) u_{hom} (s, x) d s = S (t, x) \\ u_{hom} (0, x) = ⟨ u_{in} ⟩ (x) \end{matrix} \end{matrix}$ where the memory kernel is given by $\begin{matrix} (12) & K (τ, x) : = \int_{Y} σ (x, y) e^{- τ L_{σ}} L_{1} σ (x, y) d y \end{matrix}$ and the source term is given by $\begin{array}{rcl} S (t, x) & : = & ⟨ f ⟩ (t, x) - \int_{0}^{t} \int_{Y} σ (x, y) e^{- (t - s) L_{σ}} L_{1} f (s, x, y) d y d s \\ (13) & - \int_{Y} σ (x, y) e^{- t L_{σ}} L_{1} u_{in} (x, y) d y . \end{array}$
Remark 6.
The expression (12) should be understood as the action of the semigroup $e^{- t L_{σ}}$ on the quantity $σ (y) - ⟨ σ ⟩$ followed by the multiplication of the periodic function σ. The thus obtained result is then integrated over the periodic variable y in the unit cell Y. Similar interpretation can be given for the terms involved in the expression (13) for the source $S (t, x)$ .
Proof of Theorem 5.
To employ compactness results of the two-scale convergence theory, we first derive uniform (w.r.t ε) estimates on the solution $u^{ε} (t, x)$ . Owing to the simplicity of the evolution equation (10), we can explicitly solve it and get $\begin{array}{l} u^{ε} (t, x) = u_{in} (x, \frac{x}{ε}) e^{- σ^{ε} (x) t} + \int_{0}^{t} e^{- σ^{ε} (x) (t - s)} f (s, x, \frac{x}{ε}) d s . \end{array}$ The regularity properties of the initial condition $u_{in}$ and of the source term f, together with the fact that $σ^{ε} ⩾ 0$ , imply that $\begin{matrix} {‖ u^{ε} ‖}_{L^{\infty} ((0, T); L^{2} (Ω))} ⩽ C < \infty, \end{matrix}$ where the constant is C independent of the parameter ε. Therefore, by Theorem 3, there exists a subsequence $u^{ε}$ which two-scale converges to a function $u^{0} \in L^{2} ((0, T) \times Ω \times Y)$ .

Next, Passing to the limit as $ε \to 0$ in the sense of two-scale in the above expression, we obtain $\begin{array}{l} u^{0} (t, x, y) = u_{in} (x, y) e^{- σ (x, y) t} + \int_{0}^{t} e^{- σ (x, y) (t - s)} f (s, x, y) d s \end{array}$ which essentially says that the limit $u^{0}$ solves the following two-scale evolution: $\begin{matrix} (14) & \{\begin{matrix} \partial_{t} u^{0} (t, x, y) + σ (x, y) u^{0} (t, x, y) = f (t, x, y) & for (t, x, y) \in (0, T) \times Ω \times Y, \\ u^{0} (0, x, y) = u_{in} (x, y) & for (x, y) \in Ω \times Y . \end{matrix} \end{matrix}$ By Proposition 4, the sequence $u^{ε}$ converges weakly in $L^{2} ((0, T) \times Ω)$ to $\begin{matrix} u_{hom} (t, x) : = ⟨ u^{0} ⟩ (t, x) \end{matrix}$ and we can then decompose the two-scale limit into a homogeneous part and a remainder which is of zero mean over the periodic cell, i.e., $\begin{array}{l} (15) & u^{0} (t, x, y) = u_{hom} (t, x) + r (t, x, y) where ⟨ r ⟩ = 0 . \end{array}$ Substituting (15) into the evolution equation for the two-scale limit (14), we obtain $\begin{matrix} (16) & \partial_{t} u_{hom} + σ (x, y) u_{hom} + \partial_{t} r + σ (x, y) r = f (t, x, y) . \end{matrix}$ Integrating the above equation over the periodicity cell Y yields $\begin{matrix} (17) & \partial_{t} u_{hom} + ⟨ σ ⟩ (x) u_{hom} = ⟨ f ⟩ (t, x) - ⟨ σ (x, \cdot) r (t, x, \cdot) ⟩ \end{matrix}$ as the reminder r is of zero average in the y variable. Using (17) in (16) we obtain an equation for the remainder term: $\begin{array}{l} \partial_{t} r + σ (x, y) r - \int_{Y} σ (x, y) r (t, x, y) d y \\ (18) & = (⟨ σ ⟩ (x) - σ (x, y)) u_{hom} + f (t, x, y) - ⟨ f ⟩ (t, x) . \end{array}$ As a result, equations (17) and (18) yield the following coupled system for $u_{hom} (t, x)$ and $r (t, x, y)$ : $\begin{matrix} (19) & \{\begin{matrix} \partial_{t} u_{hom} + ⟨ σ ⟩ (x) u_{hom} = ⟨ f ⟩ (t, x) - ⟨ σ (x, \cdot) r (t, x, \cdot) ⟩ \\ \partial_{t} r + L_{σ} r = - u_{hom} L_{1} σ + L_{1} f \\ u_{hom} (0, x) = ⟨ u_{in} (x) ⟩ \\ r (0, x, y) = L_{1} u_{in} . \end{matrix} \end{matrix}$ In the above couple system, we can solve for the remainder term $r (t, x, y)$ in terms of the homogenized limit $u_{hom}$ as follows: $\begin{array}{rcl} r (t, x, y) & = & e^{- t L_{σ}} L_{1} u_{in} (x, y) + \int_{0}^{t} e^{- (t - s) L_{σ}} L_{1} f (s, x, y) d s \\ (20) & - \int_{0}^{t} e^{- (t - s) L_{σ}} L_{1} σ (x, y) u_{hom} (s, x) d s . \end{array}$ Substituting this expression for the remainder in the evolution for $u_{hom} (t, x)$ in (19) yields $\begin{array}{rcl} \partial_{t} u_{hom} + ⟨ σ ⟩ (x) u_{hom} & = & ⟨ f ⟩ (t, x) - \int_{Y} σ (x, y) v (t, x, y) d y - \int_{Y} σ (x, y) w (t, x, y) d y \\ = & ⟨ f ⟩ (t, x) + \int_{0}^{t} \int_{Y} σ (x, y) e^{- (t - s) L_{σ}} L_{1} σ (x, y) u_{hom} (s, x) d y d s \\ - \int_{0}^{t} \int_{Y} σ (x, y) e^{- (t - s) L_{σ}} L_{1} f (s, x, y) d y d s \\ - \int_{Y} σ (x, y) e^{- t L_{σ}} L_{1} u_{in} (x, y) d y, \end{array}$ from which the evolution equation (11) for $u_{hom} (t, x)$ easily follows. □
Remark 7.
The result in Theorem 5 can easily be extended to a system of ordinary differential equations of the form $\begin{matrix} (21) & \{\begin{matrix} \partial_{t} u^{ε} (t, x) + Σ^{ε} (x) u^{ε} (t, x) = f^{ε} (t, x) \\ u^{ε} (0, x) = u_{in} (x, \frac{x}{ε}) \end{matrix} \end{matrix}$ where the unknown $u^{ε} \in R^{n}$ . In this case, $\begin{matrix} Σ^{ε} (x) : = Σ (x, \frac{x}{ε}) = {(Σ_{i, j} (x, \frac{x}{ε}))}_{1 ⩽ i, j ⩽ n}, Σ_{i, j} \in L^{\infty} (Ω; C_{per} (Y)) \end{matrix}$ and we assume that $Σ^{ε}$ is diagonalizable in the sense that there exists $P \in R^{n \times n}$ invertible and $D^{ε} \in R^{n \times n}$ diagonal such that $Σ^{ε} = P D^{ε} P^{- 1}$ . The right hand side is now $\begin{matrix} f^{ε} (t, x) : = f (t, x, \frac{x}{ε}) = {(f_{i} (t, x, \frac{x}{ε}))}_{1 ⩽ i ⩽ n}, f_{i}^{ε} \in L^{\infty} ((0, T) \times Ω; C_{per} (Y)) . \end{matrix}$ Since $w^{ε} = P^{- 1} u^{ε} = {(w_{i}^{ε})}_{1 ⩽ i ⩽ n}$ satisfy for $1 ⩽ i ⩽ n$ , $\begin{matrix} (22) & \{\begin{matrix} \partial_{t} w_{i}^{ε} (t, x) + D_{i, i} (x \frac{x}{ε}) w_{i}^{ε} (t, x) = {(P^{- 1} f)}_{i} (t, x, \frac{x}{ε}) \\ w_{i}^{ε} (0, x) = {(P^{- 1} u_{in})}_{i} (x, \frac{x}{ε}), \end{matrix} \end{matrix}$ we can now apply Theorem 5 to problem (22) and derive the homogenized equation and corresponding solution $w_{i, hom}$ . We derive the final result by undoing the change of variables, namely by computing $u_{hom} = P w_{hom}$ .

Note that our extension required that Σ be diagonalizable in $R$ . At the end of the paper (see Appendix x), we present an extension of our ODE setting (10) to the case where σ and f lie in the complex plane which allows us to work with any matrix Σ.

We next show that expression (12) of our memory kernel is consistent with the one given by Tartar in (9). For this, note that when $σ (x, y) = σ (y)$ , i.e., when the coefficient is purely periodic, the memory kernel $K$ takes the form $\begin{array}{l} (23) & K (τ) : = \int_{Y} σ (y) e^{- τ L_{σ}} (σ - ⟨ σ ⟩) (y) d y . \end{array}$ The following result shows that, in this setting, we have $\hat{K} = \hat{M}$ , which means that our expression coincides with that of Tartar. Theorem 8.
For any $p > 0$ , $\begin{matrix} \hat{K} (p) = \hat{M} (p) . \end{matrix}$
Proof.
The Laplace transform of our expression (23) for the memory kernel is $\begin{array}{l} \hat{K} (p) = \int_{0}^{\infty} e^{- p t} K (t) d t = \int_{0}^{\infty} \int_{Y} σ (y) e^{- p t} e^{- t L_{σ}} (σ - ⟨ σ ⟩) (y) d y d t \end{array}$ and since the Laplace transform of a semigroup yields the corresponding resolvent, we have $\begin{array}{l} (24) & \hat{K} (p) = \int_{Y} σ (y) {[p + L_{σ}]}^{- 1} (σ - ⟨ σ ⟩) (y) d y . \end{array}$ We now prove that $\begin{matrix} (25) & {[p + L_{σ}]}^{- 1} (σ - ⟨ σ ⟩) (y) = \frac{σ (y) - ⟨ σ ⟩}{p + σ} - \frac{B (p)}{p + σ} \int_{Y} \frac{σ (y) - ⟨ σ ⟩}{p + σ (y)} d y . \end{matrix}$ For this, we consider the equation $\begin{array}{l} [p + L_{σ}] g (y) = f (y), y \in Y \end{array}$ for a given $p > 0$ and a given measurable function f of zero mean. Averaging the equation in the y variable yields that it is necessary to have the solution $g (y)$ to be of zero average as well. Hence, for zero average functions, we observe that $\begin{array}{rcl} [p + L_{σ}] g (y) & = & p + σ (y) g (y) - \int_{Y} σ (y) g (y) d y \\ = & p + σ (y) g (y) - \int_{Y} (σ (y) + p) g (y) d y = L_{p + σ} g (y) . \end{array}$ A simple inspection reveals that a general solution to $\begin{matrix} L_{p + σ} g (y) = f (y) \end{matrix}$ is given by $\begin{array}{l} g (y) = \frac{f (y)}{p + σ (y)} + \frac{C}{p + σ (y)} \end{array}$ where the constant C needs to be chosen such that $g (y)$ is of zero average. Hence $\begin{matrix} C = - B (p) \int_{Y} \frac{f (y)}{p + σ (y)} d y with B (p) : = {(\int_{Y} \frac{d y}{p + σ (y)})}^{- 1} \end{matrix}$ and we have found the explicit expression for the solution $g (y)$ .

As a result, by taking $f (y) = σ - ⟨ σ ⟩$ , we obtain the solution to be $\begin{array}{l} {[p + L_{σ}]}^{- 1} (σ - ⟨ σ ⟩) (y) = g (y) = \frac{σ (y) - ⟨ σ ⟩}{p + σ} - \frac{B (p)}{p + σ} \int_{Y} \frac{σ (y) - ⟨ σ ⟩}{p + σ (y)} d y \end{array}$ as anticipated in (25). Using the above observation, the expression (24) yields $\begin{matrix} (26) & \begin{matrix} \hat{K} (p) & = \int_{Y} σ (y) (\frac{σ (y) - ⟨ σ ⟩}{p + σ} - \frac{B (p)}{p + σ} \int_{Y} \frac{σ (y) - ⟨ σ ⟩}{p + σ (y)} d y) d y \\ = \int_{Y} (σ (y) + p) (\frac{σ (y) - ⟨ σ ⟩}{p + σ} - \frac{B (p)}{p + σ} \int_{Y} \frac{σ (y) - ⟨ σ ⟩}{p + σ (y)} d y) d y \end{matrix} \end{matrix}$ as the function $g (y)$ is of zero mean. Observe that the following quotient is of zero mean $\begin{matrix} \frac{p + σ (y) - B (p)}{p + σ (y)}, \end{matrix}$ which implies the following identity $\begin{array}{l} (27) & \int_{Y} \frac{σ (y)}{p + σ (y)} d y = 1 - \frac{p}{B (p)} . \end{array}$ Using it in (26) yields $\begin{array}{l} \hat{K} (p) = ⟨ σ ⟩ - B (p) \int_{Y} \frac{σ (y)}{p + σ (y)} d y = p + ⟨ σ ⟩ - B (p) = \hat{M} (p), \end{array}$ thus proving the equivalence. □

3. Homogenization in energy for the linear Boltzmann equation

We consider the homogenization of the multi-scale linear Boltzmann model (3), where the optical parameters $σ^{ε}$ and $κ^{ε}$ present oscillations in the energy variable that are modeled in accordance to experimental observations (see the Introduction). The first step is the proof of some a priori bounds, which will then be used to prove our main homogenization result. The proofs for the a priori bounds involve classical techniques and similar results can be found in [6, Chapter 21].

In this section, $⟨ \cdot ⟩$ denotes integral over the interval $(0, 1)$ , i.e. averaging over the periodic cell in the energy variable: $\begin{matrix} ⟨ v ⟩ : = \int_{0}^{1} v (y) d y for all v \in L^{1} (0, 1) . \end{matrix}$

3.1. A priori bounds of the solution

In the following, we will sometimes use the shorthand notation $\begin{matrix} V = S^{d - 1} \times (E_{min}, E_{max}) \end{matrix}$ We first prove a positivity property of the Boltzmann operator $\begin{array}{l} Q^{ε} f : = σ^{ε} f - \int_{V} κ^{ε} (x, ω \cdot ω^{'}, E, E^{'}) f (x, ω^{'}, E^{'}) d ω^{'} d E^{'}, \forall f \in L^{2} (Ω \times V) . \end{array}$

Proposition 9.
If $(σ^{ε}, κ^{ε})$ satisfy assumptions ( 4 ) and ( 5 ), then for all $ε > 0$ and all $f \in L^{2} (Ω \times V)$ , $\begin{matrix} (28) & {(Q^{ε} f, f)}_{L^{2} (Ω \times V)} ⩾ α ‖ f ‖_{L^{2} (Ω \times V)}^{2} . \end{matrix}$
Proof.
It follows from the Cauchy–Schwarz inequality and the definition of ${\bar{κ}}^{ε}$ and ${\tilde{κ}}^{ε}$ in (6) that $\begin{array}{l} \int_{V} \int_{V} f (x, ω, E) f (x, ω^{'}, E^{'}) κ^{ε} (x, ω \cdot ω^{'}, E, E^{'}) d ω d E d ω^{'} d E^{'} \\ ⩽ {(\int_{V} \int_{V} {| f (x, ω, E) |}^{2} κ^{ε} (x, ω \cdot ω^{'}, E, E^{'}) d ω d E d ω^{'} d E^{'})}^{1 / 2} \\ \times {(\int_{V} \int_{V} {| f (x, ω^{'}, E^{'}) |}^{2} κ^{ε} (x, ω \cdot ω^{'}, E, E^{'}) d ω d E d ω^{'} d E^{'})}^{1 / 2} \\ = {(\int_{V} {| f (x, ω, E) |}^{2} {\bar{κ}}^{ε} (x, ω, E) d ω d E)}^{1 / 2} {(\int_{V} {| f (x, ω, E) |}^{2} {\tilde{κ}}^{ε} (x, ω, E) d ω d E)}^{1 / 2} . \end{array}$ We then have the following, thanks to Young’s inequality and our assumption (5): $\begin{array}{l} {(Q^{ε} f, f)}_{L^{2} (Ω \times V)} \\ ⩾ \int_{Ω} d x \int_{V} {| f (x, ω, E) |}^{2} σ^{ε} (x, ω, E) d ω d E \\ - \int_{Ω} {(\int_{V} {| f (x, ω, E) |}^{2} {\bar{κ}}^{ε} (x, ω, E) d ω d E)}^{1 / 2} {(\int_{V} {| f (x, ω, E) |}^{2} {\tilde{κ}}^{ε} (x, ω, E) d ω d E)}^{1 / 2} d x \\ ⩾ α ‖ f ‖_{L^{2} (Ω \times V)}^{2} . \end{array}$ □

Using Proposition 9 in the energy arguments, one can obtain the apriori bounds in the following lemma. Lemma 10.
If $(σ^{ε}, κ^{ε})$ satisfy assumptions ( 4 ) and ( 5 ), then there exists $C > 0$ such that for all $ε > 0$ , the solution $φ^{ε}$ to ( 3 ) satisfies $\begin{array}{l} (29) & {‖ φ^{ε} ‖}_{L^{\infty} ((0, T); L^{2} (Ω \times V))} ⩽ C and {‖ φ^{ε} ‖}_{L^{2} ((0, T) \times Ω \times V)} ⩽ C . \end{array}$

3.2. Homogenization result

Following a similar pathway as in our initial ODE model (see Section 2), we derive the following result for the scaled Equation (3), which is proved, without loss of generality, by taking the mass of the neutron $m = 2$ . Since we focus on the rapid oscillations in energy, in what follows we assume no space dependence in the optical parameters.

Theorem 11.
Let $φ^{ε} = φ^{ε} (t, x, ω, E)$ be the solution of the equation $\begin{matrix} (30) & \{\begin{matrix} \partial_{t} φ^{ε} + \sqrt{E} ω \cdot \nabla_{x} φ^{ε} + σ^{ε} (ω, E) φ^{ε} \\ - \int_{E_{\min}}^{E_{\max}} \int_{| ω^{'} | = 1} κ^{ε} (ω \cdot ω^{'}, E, E^{'}) φ^{ε} (ω^{'}, E^{'}) d ω^{'} d E^{'} = 0 \\ φ^{ε} (0, x, ω, E) = φ_{in}^{ε} (x, ω, E), \\ φ^{ε} (t, x, ω, E) = 0 \forall t > 0 \\ and for (x, ω) \in Γ_{-} : = {(x, ω) \in \partial Ω \times S^{d - 1} : n (x) \cdot ω < 0}, \end{matrix} \end{matrix}$ where the coefficients and the data are of the form $\begin{array}{l} σ^{ε} (ω, E) : = \sqrt{E} σ (ω, E, \frac{E}{ε}) with σ (ω, E, y) \in L^{\infty} (V; C_{per} (0, 1)) \\ φ_{in}^{ε} (x, ω, E) : = φ_{in} (x, ω, E, \frac{E}{ε}) with φ_{in} (x, ω, E, y) \in L^{2} (Ω \times V; C_{per} (0, 1)) \\ κ^{ε} (ω \cdot ω^{'}, E, E^{'}) : = \sqrt{E} κ_{1} (ω \cdot ω^{'}, E) κ_{2} (ω \cdot ω^{'}, E^{'}, \frac{E^{'}}{ε}) \\ with κ_{1} (η, E) \in L^{\infty} ([- 1, 1] \times [E_{min}, E_{max}]) \\ and κ_{2} (η, E^{'}, y^{'}) \in L^{\infty} ([- 1, 1] \times [E_{min}, E_{max}]; C_{per} (0, 1)) . \end{array}$ Then, $\begin{matrix} φ^{ε} ⇀ φ_{hom} weakly in L^{2} ((0, T) \times Ω \times V) \end{matrix}$ and $φ_{hom}$ satisfies the following partial integro-differential equation $\begin{array}{l} \partial_{t} φ_{hom} + \sqrt{E} ω \cdot \nabla_{x} φ_{hom} + \sqrt{E} ⟨ σ ⟩ φ_{hom} \\ - \int_{E_{\min}}^{E_{\max}} \int_{S^{d - 1}} \sqrt{E} κ_{1} (ω \cdot ω^{'}, E) \int_{0}^{1} κ_{2} (ω \cdot ω^{'}, E^{'}, y^{'}) φ_{hom} d y^{'} d ω^{'} d E^{'} \\ = \int_{E_{\min}}^{E_{\max}} \int_{S^{d - 1}} \sqrt{E} κ_{1} (ω \cdot ω^{'}, E) \int_{0}^{1} κ_{2} (ω \cdot ω^{'}, E^{'}, y^{'}) \\ \times [e^{- t \sqrt{E^{'}} L_{σ}} L_{1} φ_{in} - \int_{0}^{t} e^{- (t - s) \sqrt{E^{'}} L_{σ}} \sqrt{E^{'}} L_{1} σ (ω^{'}, E^{'}, y^{'}) φ_{hom} d s] d y^{'} d ω^{'} d E^{'} \\ - \sqrt{E} \int_{0}^{1} σ (ω, E, y) [e^{- t \sqrt{E} L_{σ}} L_{1} φ_{in} - \int_{0}^{t} e^{- (t - s) \sqrt{E} L_{σ}} \sqrt{E} L_{1} σ (ω, E, y) φ_{hom} d s] d y, \end{array}$ with initial condition $\begin{matrix} φ_{hom} (0, x, ω, E) = ⟨ φ_{in} (x, ω, E, \cdot) ⟩ \end{matrix}$ and zero absorption condition at the in-flux phase-space boundary.
Proof.
Following the same strategy as in Theorem 5, thanks to the method of characteristics, we can write the explicit solution to (30) as $\begin{matrix} φ^{ε} (t, x, ω, E) \\ = φ_{in}^{ε} (x - \sqrt{E} ω t, ω, E) e^{- t σ^{ε} (ω, E)} \\ + \int_{0}^{t} \int_{V} κ^{ε} (ω \cdot ω^{'}, E, E^{'}) e^{- (t - s) σ^{ε} (ω, E)} φ^{ε} (s, x - (t - s) \sqrt{E} ω, ω^{'}, E^{'}) d ω^{'} d E^{'} d s . \end{matrix}$ Making the change of variables $x - \sqrt{E} ω t = r$ , we get $\begin{matrix} φ^{ε} (t, r + \sqrt{E} ω t, ω, E) \\ = φ_{in}^{ε} (r, ω, E) e^{- t σ^{ε} (ω, E)} \\ + \int_{0}^{t} \int_{V} κ^{ε} (ω \cdot ω^{'}, E, E^{'}) e^{- (t - s) σ^{ε} (ω, E)} φ^{ε} (s, r + s \sqrt{E} ω, ω^{'}, E^{'}) d ω^{'} d E^{'} d s . \end{matrix}$ Defining $\begin{matrix} φ^{ε} (t, r + \sqrt{E} ω t, ω, E) = : ψ^{ε} (t, r, ω, E) \end{matrix}$ we have the following expression for the unknown $ψ^{ε}$ : $\begin{matrix} (31) & \begin{matrix} ψ^{ε} (t, r, ω, E) = & φ_{in}^{ε} (r, ω, E) e^{- t σ^{ε} (ω, E)} \\ + \int_{0}^{t} \int_{V} κ^{ε} (ω \cdot ω^{'}, E, E^{'}) e^{- (t - s) σ^{ε} (ω, E)} ψ^{ε} (s, r, ω^{'}, E^{'}) d ω^{'} d E^{'} d s . \end{matrix} \end{matrix}$ We can pass to the limit in the expression (31) as $ε \to 0$ under the assumptions on the scattering kernel which have been listed in the statement of the theorem.

As the term of interest in the impending asymptotic analysis is the integral term in (31), we shall detail the asymptotic procedure for that term alone. Multiplying the said integral term by a test function $g (E, \frac{E}{ε})$ and integrating over the E variable and exploiting the structure of the scattering kernel $κ^{ε}$ , we have $\begin{array}{l} \int_{E_{\min}}^{E_{\max}} \int_{0}^{t} \int_{E_{\min}}^{E_{\max}} \int_{S^{d - 1}} κ^{ε} (ω \cdot ω^{'}, E, E^{'}) e^{- (t - s) σ^{ε} (ω, E)} ψ^{ε} (s, r, ω^{'}, E^{'}) g (E, \frac{E}{ε}) d ω^{'} d E^{'} d s d E \\ (32) & = \int_{E_{\min}}^{E_{\max}} \int_{0}^{t} \int_{S^{d - 1}} \sqrt{E} κ_{1} (ω \cdot ω^{'}, E) e^{- (t - s) σ^{ε} (ω, E)} w^{ε} (s, r, ω, ω^{'}) g (E, \frac{E}{ε}) d ω^{'} d s d E, \end{array}$ where we have used the notation $\begin{matrix} w^{ε} (s, r, ω, ω^{'}) : = \int_{E_{\min}}^{E_{\max}} κ_{2} (ω \cdot ω^{'}, E^{'}, \frac{E^{'}}{ε}) ψ^{ε} (s, r, ω^{'}, E^{'}) d E^{'} . \end{matrix}$ For every $(s, r, ω, ω^{'}) \in (0, T) \times Ω \times S^{d - 1} \times S^{d - 1}$ , $\begin{matrix} w^{ε} (s, r, ω, ω^{'}) \to w^{0} (s, r, ω, ω^{'}) = \int_{E_{\min}}^{E_{\max}} \int_{0}^{1} κ_{2} (ω \cdot ω^{'}, E^{'}, y^{'}) ψ^{0} (s, r, ω^{'}, E^{'}, y^{'}) d y^{'} d E^{'} \end{matrix}$ pointwise as $ε \to 0$ , where $ψ^{0}$ is the two-scale limit of $ψ^{ε}$ . Therefore we have two-scale convergence of (32), namely $\begin{array}{l} \int_{E_{\min}}^{E_{\max}} \int_{0}^{t} \int_{S^{d - 1}} \sqrt{E} κ_{1} (ω \cdot ω^{'}, E) e^{- (t - s) σ^{ε} (ω, E)} w^{ε} (s, y, ω, ω^{'}) g (E, \frac{E}{ε}) d ω^{'} d s d E \\ ⇀ \int_{E_{\min}}^{E_{\max}} \int_{0}^{t} \int_{S^{d - 1}} \int_{0}^{1} \sqrt{E} κ_{1} (ω \cdot ω^{'}, E) e^{- (t - s) \sqrt{E} σ (ω, E, y)} w^{0} (s, r, ω, ω^{'}) g (E, y) d y d ω^{'} d s d E . \end{array}$ As a result of this calculation, passing to the limit in (31) in the sense of two-scale yields $\begin{array}{l} ψ^{0} (t, r, ω, E, y) \\ = φ_{in} (r, ω, E, y) e^{- t \sqrt{E} σ (ω, E, y)} \\ + \int_{0}^{t} \int_{E_{\min}}^{E_{\max}} \int_{S^{d - 1}} \int_{0}^{1} \sqrt{E} κ_{1} (ω \cdot ω^{'}, E) e^{- (t - s) \sqrt{E} σ (ω, E, y)} \\ \times κ_{2} (ω \cdot ω^{'}, E^{'}, y^{'}) ψ^{0} (s, r, ω^{'}, E^{'}, y^{'}) d y^{'} d ω^{'} d E^{'} d s \end{array}$ which implies that $ψ^{0}$ solves the two-scale evolution $\begin{matrix} \{\begin{matrix} \partial_{t} ψ^{0} (t, r, ω, E, y) + \sqrt{E} σ (ω, E, y) ψ^{0} (t, r, ω, E, y) \\ = \int_{E_{\min}}^{E_{\max}} \int_{S^{d - 1}} \sqrt{E} κ_{1} (ω \cdot ω^{'}, E) \int_{0}^{1} κ_{2} (ω \cdot ω^{'}, E^{'}, y^{'}) ψ^{0} (s, r, ω^{'}, E^{'}, y^{'}) d y^{'} d ω^{'} d E^{'}, \\ ψ^{0} (0, r, ω, E, y) = φ_{in} (r, ω, E, y) . \end{matrix} \end{matrix}$ Proceeding in the same way as in the proof of Theorem 5, we decompose the two-scale limit $ψ^{0}$ into a homogeneous part and a remainder, which is of zero mean over the periodic cell, i.e., $\begin{array}{l} ψ^{0} (t, r, ω, E, y) = ψ_{hom} (t, r, ω, E) + ρ (t, r, ω, E, y) \end{array}$ where $\begin{matrix} ψ_{hom} (t, r, ω, E) : = ⟨ ψ^{0} ⟩ (t, r, ω, E) and ⟨ ρ ⟩ = 0 . \end{matrix}$ We obtain the following coupled system for $ψ_{hom}$ and ρ: $\begin{matrix} (33) & \{\begin{matrix} \partial_{t} ψ_{hom} + \sqrt{E} ⟨ σ ⟩ ψ_{hom} - S (ψ_{hom}) = S (ρ) - \sqrt{E} ⟨ σ ρ ⟩ \\ \partial_{t} ρ + \sqrt{E} (σ ρ - ⟨ σ ρ ⟩) = \sqrt{E} (⟨ σ ⟩ - σ) ψ_{hom} \\ ψ_{hom} (0, r, ω, E) = ⟨ φ_{in} (r, ω, E, \cdot) ⟩ \\ ρ (0, r, ω, E, y) = L_{1} φ_{in} (r, ω, E, y), \end{matrix} \end{matrix}$ where $\begin{matrix} S (ξ) : = \int_{E_{\min}}^{E_{\max}} \int_{S^{d - 1}} \sqrt{E} κ_{1} (ω \cdot ω^{'}, E) \int_{0}^{1} κ_{2} (ω \cdot ω^{'}, E^{'}, y^{'}) ξ (s, r, ω^{'}, E^{'}, y^{'}) d y^{'} d ω^{'} d E^{'} \end{matrix}$ for any $ξ \in L^{2} ((0, T) \times Ω \times V; C_{per} (0, 1))$ . The coupled system written above has the same structure as system (19). Hence, following the same lines as in Section 2, we can deduce $\begin{array}{rcl} ρ (t, r, ω, E, y) & = & e^{- t \sqrt{E} L_{σ}} L_{1} φ_{in} (r, ω, E, y) \\ - \int_{0}^{t} e^{- (t - s) \sqrt{E} L_{σ}} \sqrt{E} L_{1} σ (ω, E, y) ψ_{hom} (s, r, ω, E) d s \end{array}$ and, finally, substituting this expression in the evolution equation for $ψ_{hom}$ , we obtain $\begin{matrix} \{\begin{matrix} \partial_{t} ψ_{hom} + \sqrt{E} ⟨ σ ⟩ ψ_{hom} - S (ψ_{hom}) = S (ρ) - \sqrt{E} ⟨ σ ρ ⟩ \\ ψ_{hom} (0, r, ω, E) = ⟨ φ_{in} (r, ω, E, \cdot) ⟩, \end{matrix} \end{matrix}$ that is $\begin{array}{l} \partial_{t} ψ_{hom} + \sqrt{E} ⟨ σ ⟩ ψ_{hom} - \int_{E_{\min}}^{E_{\max}} \int_{S^{d - 1}} \sqrt{E} κ_{1} (ω \cdot ω^{'}, E) \int_{0}^{1} κ_{2} (ω \cdot ω^{'}, E^{'}, y^{'}) ψ_{hom} d y^{'} d ω^{'} d E^{'} \\ = \int_{E_{\min}}^{E_{\max}} \int_{S^{d - 1}} \sqrt{E} κ_{1} (ω \cdot ω^{'}, E) \int_{0}^{1} κ_{2} (ω \cdot ω^{'}, E^{'}, y^{'}) \\ \times [e^{- t \sqrt{E^{'}} L_{σ}} L_{1} φ_{in} - \int_{0}^{t} e^{- (t - s) \sqrt{E^{'}} L_{σ}} \sqrt{E^{'}} L_{1} σ (ω^{'}, E^{'}, y^{'}) ψ_{hom} d s] d y^{'} d ω^{'} d E^{'} \\ - \sqrt{E} \int_{0}^{1} σ (ω, E, y) [e^{- t \sqrt{E} L_{σ}} L_{1} φ_{in} - \int_{0}^{t} e^{- (t - s) \sqrt{E} L_{σ}} \sqrt{E} L_{1} σ (ω, E, y) ψ_{hom} d s] d y, \end{array}$ with initial condition $\begin{matrix} ψ_{hom} (0, r, ω, E) = ⟨ φ_{in} (r, ω, E, \cdot) ⟩ . \end{matrix}$ By going back to the original variables, we finally obtain the thesis of the theorem. □
Remark 12.
Some comments are in order with regards to the assumption on the optical parameters $σ^{ε}$ and $κ^{ε}$ made in Theorem 11.
The assumption of separability $\begin{matrix} κ^{ε} (ω \cdot ω^{'}, E, E^{'}) : = \sqrt{E} κ_{1} (ω \cdot ω^{'}, E) κ_{2} (ω \cdot ω^{'}, E^{'}, \frac{E^{'}}{ε}) \end{matrix}$ has simplified our computations in the proof. It also lead to a relatively simpler homogenized model.

In the above separable structure, we can further allow the factor $κ_{1}$ to oscillate in the E-variable. The proof of Theorem 11 can be reworked in this case, at the price of arriving at a more complex two-scale system than (33). The memory structure remains the same but with additional terms.

It is apparent from the proof of Theorem 11 that the energy oscillations in σ, not those in the scattering kernel, resulted in the memory effects in the homogenized limit.

4. Numerical tests

We present some results on admittedly simple cases, aiming primarily at studying the convergence rate towards the homogenized solution as the period of the oscillating terms tends to zero.

4.1. Collision term κ inside the integral

We consider the problem $\begin{matrix} (34) & \{\begin{matrix} \partial_{t} φ^{ε} (t, E) + σ (\frac{E}{ε}) φ^{ε} (t, E) = \int_{E_{min}}^{E_{max}} κ (\frac{E^{'}}{ε}) φ^{ε} (t, E^{'}) d E^{'} \\ φ^{ε} (0, E) = φ_{in} (E) \end{matrix} \end{matrix}$ for $(t, E) \in [0, T] \times (E_{min}, E_{max})$ . This equation is a simple instance of the linear Boltzmann equation (30) with neither space nor angular dependence. In the following, we set $(E_{min}, E_{max}) = (0, 1)$ and the final time to $T = 10$ . By Theorem 11, we know that $φ^{ε}$ weakly converges in $L^{2} ((0, T) \times (E_{min}, E_{max}))$ to $φ_{hom}$ under some assumptions on the optical parameters. We study here the convergence rate for different functions σ, κ and initial value $φ_{in}$ . For this, we introduce the family of orthogonal Legendre polynomials in $L^{2} (E_{min}, E_{max})$ which we denote by ${ℓ_{k}}_{k ⩾ 0}$ and define the modes $\begin{matrix} (35) & m_{k}^{ε} (t) : = {(φ^{ε} (t, \cdot), ℓ_{k} (\cdot))}_{L^{2} ([E_{min}, E_{max}])}, k ⩾ 0, \end{matrix}$ of the solution for $t \in [0, T]$ . Likewise, $\begin{matrix} (36) & m_{k}^{hom} (t) : = {(φ_{hom} (t, \cdot), ℓ_{k} (\cdot))}_{L^{2} ([E_{min}, E_{max}])}, k ⩾ 0, \end{matrix}$ are the modes of the homogenized solution $φ_{hom}$ . We know that $φ^{ε} ⇀ φ_{hom}$ weakly in $L^{2} ([E_{min}, E_{max}])$ if $\begin{matrix} (37) & e_{k}^{ε} : = max_{t \in [0, T]} | m_{k}^{ε} (t) - m_{k}^{hom} (t) | \to 0, as ε \to 0 \forall k ⩾ 0 . \end{matrix}$ In the following, we use this sufficient condition to give numerical evidence of weak convergence. Note that, in practice, we cannot inspect all the values $k ⩾ 0$ so we present results for the first eight modes.

First example. We take $\begin{matrix} σ (y) = 2 + \frac{1}{2} sin (2 π y), κ (y^{'}) = 1 + \frac{1}{2} sin (2 π y^{'}), φ_{in} (y) = 1 + sin (2 π y) . \end{matrix}$ We compute the numerical solution $φ^{ε}$ with an explicit Runge Kutta method of order four with time step $d t = T / 50$ . We take a uniform mesh of size h for the discretization of the oscillatory variable. We refine h depending on the value of ε to be able to resolve the oscillations. Our choice is to take $h = ε / 100$ . This discretization might be an overkill and one could probably achieve a good accuracy with less degrees of freedom.

We next give numerical evidence that $φ^{ε}$ weakly converges to $φ_{hom}$ in $L^{2} ((0, T) \times (E_{min}, E_{max}))$ through the plot of the error $e_{k}^{ε}$ in Fig. 2. The plot shows that $e_{k}^{ε}$ tends to 0 as $ε \to 0$ for $k = 0, \dots, 7$ at a rate slightly lower than ε. Next, we illustrate some kind of “strong” two-scale convergence by showing convergence in the norm: $\begin{matrix} (38) & | {‖ φ^{ε} ‖}_{L^{2} ((0, T) \times (E_{min}, E_{max}))} - {‖ φ^{0} ‖}_{L^{2} ((0, T) \times (E_{min}, E_{max}) \times (0, 1))} | \to 0 as ε \to 0 . \end{matrix}$ Figure 2 shows that this is indeed the case. We refer the readers to [2, Theorem 1.8 on p. 1488] for more details on strong two-scale convergence.

Fig. 2.

Function κ inside the integral. First example.

Second example. Taking the first example as a starting point, we study the impact on the convergence rate of the initial value and consider $\begin{matrix} σ (y) = 2 + \frac{1}{2} sin (2 π y), k (y^{'}) = 1 + \frac{1}{2} sin (2 π y^{'}), φ_{in} (y) = 1 + 1_{{y ⩽ \frac{1}{2}}} (y) . \end{matrix}$ Here and in the following $1_{A}$ denotes the indicator function of a set $A \subset [0, 1]$ , i.e., $1_{A} (y) = 1$ if $y \in A$ and 0 otherwise. As a result, in our example $φ_{in} (y)$ has a jump at $y = \frac{1}{2}$ . Figure 3 shows $e_{k}^{ε}$ as a function of ε. We observe that the discontinuity in the initial condition does not affect the overall convergence rate. Here again, we observe that (38) holds.

Fig. 3.

Function κ inside the integral. Second example.

Third example. Taking once again the first example as a starting point, we now investigate the impact on the convergence rate that oscillatory and discontinuous functions σ and k can have. We consider $\begin{matrix} σ (y) = 2 + \frac{1}{2} 1_{{sin (2 π y) ⩾ 0}} (y), κ (y^{'}) = 1 + \frac{1}{2} 1_{{sin (2 π y^{'}) ⩾ 0}} (y^{'}), φ_{in} (y) = 1 + sin (2 π y) \end{matrix}$ and Fig. 4 shows the convergence of $e_{k}^{ε}$ and of the norm difference. We see that weak convergence occurs at a similar rate as before, indicating that the regularity of the oscillatory functions might have a mild impact on the convergence.

Fig. 4.

Function κ inside the integral. Third example.

4.2. Collision term κ outside the integral

It is interesting to compare the converge rate of problem (34) with $\begin{matrix} (39) & \{\begin{matrix} \partial_{t} φ^{ε} (t, E) + σ (\frac{E}{ε}) φ^{ε} (t, E) = κ (\frac{E}{ε}) \int_{E_{min}}^{E_{max}} φ^{ε} (t, E^{'}) d E^{'} \\ φ^{ε} (0, E) = φ_{in} (E) \end{matrix} \end{matrix}$ where the collision kernel κ is outside of the integral. Figure 5 gives numerical evidence of weak and strong convergence of (39) for the same three choices of σ, κ and $φ_{in}$ that were studied in the previous section. Note that, due to the fact that κ is now outside of the integral, the convergence rate is higher and occurs at a rate of about $ε^{2}$ . Note also that in this case the regularity of the initial condition and the oscillatory coefficients significantly affect the convergence rate. For all cases, we also observe strong convergence (the plots are very similar to the ones of the previous section – so we omit them).

Fig. 5.

Function κ outside the integral: convergence rate in ε of the error $e_{k}^{ε}$ .

Footnotes

Acknowledgements

The authors would like to thank the referees for helpful discussion and for their valuable comments. This work was supported by the ANR project Kimega (ANR-14-ACHN-0030-01) and by the Italian Ministry of Education, University and Research (Dipartimenti di Eccellenza program 2018-2022, Dipartimento di Matematica ‘F. Casorati’, Università degli Studi di Pavia). We also received funding from a PEPS grant from the INSMI (Institut National des Sciences Mathématiques et de leurs Interactions). H. Hutridurga and O. Mula started to work on this paper while enjoying the kind hospitality of the Hausdorff Institute for Mathematics in Bonn during the trimester program on multiscale problems in 2017. H. Hutridurga also acknowledges the support of the EPSRC programme grant “Mathematical fundamentals of Metamaterials for multiscale Physics and Mechanics” (EP/L024926/1). The authors thank Vishal Vasan (ICTS) for posing the question on harmonic oscillator addressed here in the Appendix.

Harmonic oscillator

Taking inspiration from Tartar’s work on memory effects [19], we analyse the asymptotic behaviour of a differential equation for the unknown $u^{ε} : [0, \infty) \times Ω \mapsto C$ $\begin{matrix} (40) & \{\begin{matrix} \partial_{t} u^{ε} (t, x) + i b^{ε} (x) u^{ε} (t, x) = 0 & for t > 0, x \in Ω, \\ u^{ε} (0, x) = u_{in} (x) & for x \in Ω . \end{matrix} \end{matrix}$ The coefficient family ${b^{ε} (x)}$ is assumed to be a family of bounded real-valued function. Taking $u^{ε} = w^{ε} + i v^{ε}$ , we obtain a system of equations for the real and imaginary parts: $\begin{matrix} (41) & \{\begin{matrix} \partial_{t} w^{ε} (t, x) - b^{ε} (x) v^{ε} (t, x) = 0 & for t > 0, x \in Ω, \\ \partial_{t} v^{ε} (t, x) + b^{ε} (x) w^{ε} (t, x) = 0 & for t > 0, x \in Ω, \\ (w^{ε}, v^{ε}) (0, x) = (w_{in}, v_{in}) (x) & for x \in Ω, \end{matrix} \end{matrix}$ where we have used that $u_{in} (x) = w_{in} (x) + i v_{in} (x) for x \in Ω$ . Define the vector function $U^{ε} : [0, \infty) \times Ω \mapsto R^{2}$ $\begin{matrix} U^{ε} (t, x) : = (\begin{matrix} w^{ε} \\ v^{ε} \end{matrix}) (t, x) . \end{matrix}$ Consider a skew-symmetric matrix $\begin{matrix} A : = (\begin{matrix} 0 & 1 \\ - 1 & 0 \end{matrix}) . \end{matrix}$ Then the system of equations (41) can be rewritten as $\begin{matrix} (42) & \{\begin{matrix} \partial_{t} U^{ε} (t, x) = b^{ε} (x) A U^{ε} (t, x) & for t > 0, x \in Ω, \\ U^{ε} (0, x) = U_{in} (x) : = (\begin{matrix} w_{in} \\ v_{in} \end{matrix}) (x) & for x \in Ω . \end{matrix} \end{matrix}$ Note that the system (42) is nothing but the harmonic oscillator equation with heterogeneous coefficients. Note further that memory effects have been shown to appear in the homogenization limit for the second order ordinary differential equations (see the work of J.-S. Jiang, K.-H. Kuo and C.-K. Lin [10]).

Our proof of the above theorem, similar to the analysis of Luc Tartar in [19], goes via the Laplace transform. Furthermore, as we are making no structural assumptions (such as periodicity) on the coefficient family $b^{ε} (x)$ , we need to use the notion of Young measures. Hence, before we give the proof of Theorem 13, we recall a cornerstone result of the theory of compensated compactness which gives the existence of certain probability measures called the young mesaures [18].

References

Agoshkov, Boundary Value Problems for Transport Equations. Modeling and Simulation in Science, Engineering and Technology, Birkhäuser, Boston, 1998.

Allaire, Homogenization and two-scale convergence, SIAM Journal on Mathematical Analysis 23(6) (1992), 1482–1518. doi:10.1137/0523084.

Bernard,

Golse and

Salvarani, Homogenization of transport problems and semigroups, Math Method Appl Sci 33(10) (2010), 1228–1234. doi:10.1002/mma.1319.

D.E.

Cullen and

G.C.

Pomraning, The multiband method in radiative transfer calculations, Journal of Quantitative Spectroscopy and Radiative Transfer 24(2) (1980), 97–117. doi:10.1016/0022-4073(80)90024-2.

Dahmen,

Gruber and

Mula, An adaptive nested source term iteration for radiative transfer equations. 2018, ArXiv e-prints https://arxiv.org/abs/1810.07035.

Dautray and

J.-L.

Lions, Mathematical Analysis and Numerical Methods for Science and Technology: Volume 6 Evolution Problems II, Springer Science & Business Media, 2012. doi:10.1007/978-3-642-58004-8.

Dumas and

Golse, Homogenization of transport equations, SIAM J Appl Math 60(4) (2000), 1447–1470. doi:10.1137/S0036139997332087.

Haut,

Ahrens,

Jonko,

Lowrie and

Till, A new multigroup method for cross-sections that vary rapidly in energy, Journal of Quantitative Spectroscopy and Radiative Transfer 187 (2017), 461–471. doi:10.1016/j.jqsrt.2016.10.019.

J.-S.

Jiang, Two-scale homogenization and memory effects of a first order differential equation, Taiwanese Journal of Mathematics (2006), 963–975. doi:10.11650/twjm/1500403887.

10.

J.-S.

Jiang,

K.-H.

Kuo and

C.-K.

Lin, On the homogenization of second order differential equations, Taiwanese Journal of Mathematics (2005), 215–236. doi:10.11650/twjm/1500407797.

11.

Kanschat, Solution of radiative transfer problems with finite elements, in: Numerical Methods in Multidimensional Radiative Transfer,

Kanschat,

Meinköhn,

Rannacher and

Wehrse, eds, Springer, Berlin Heidelberg, 2009, pp. 49–98. doi:10.1007/978-3-540-85369-5_5.

12.

Koning,

Forrest,

Kellett,

Mills,

Henriksson,

Rugama,

Bersillon,

Bouland,

Courcelle,

Duijvestijn et al., The JEFF-3.1 Nuclear Data Library-JEFF Report 21. Technical report, Organisation for Economic Co-operation and Development, 2006. https://www.oecd-nea.org/dbdata/nds_jefreports/jefreport-21/jeff21.pdf.

13.

Lewis and

Miller, Computational Methods of Neutron Transport, J. Wiley & Sons, Inc., New York, 1984.

14.

Livolant and

Jeanpierre, Autoprotection des résonnances dans les réacteurs nucléaires, Technical Report 4533, CEA, 1974. https://inis.iaea.org/collection/NCLCollectionStore/_Public/05/153/5153104.pdf.

15.

Mathiaud and

Salvarani, A numerical strategy for radiative transfer problems with higly oscillating opacities, Applied Mathematics and Computation 221 (2013), 249–256. doi:10.1016/j.amc.2013.06.057.

16.

Nguetseng, A general convergence result for a functional related to the theory of homogenization, SIAM Journal on Mathematical Analysis 20(3) (1989), 608–623. doi:10.1137/0520043.

17.

Reuss, Neutron Physics, EDP Sciences, 2008.

18.

Tartar, Compensated compactness and applications to partial differential equations, in: Nonlinear Analysis and Mechanics IV,

Knops, ed., Pitman Res. Notes Math. Ser., Vol. 39, Longman, Harlow, UK, 1979, pp. 136–211, https://www.researchgate.net/publication/239059426_Compensated_compactness_and_applications_to_partial_differential_equations.

19.

Tartar, Nonlocal effects induced by homogenization, in: Partial Differential Equations and the Calculus of Variations, Vol. II, Progr. Nonlinear Differential Equations Appl., Vol. 2, Birkhäuser Boston, Boston, MA, 1989, pp. 925–938. doi:10.1007/978-1-4615-9831-2_19.

20.

Tartar, An Introduction to Navier–Stokes Equation and Oceanography, Lecture Notes of the Unione Matematica Italiana., Vol. 1, Springer-Verlag, Berlin Heidelberg, 2006. doi:10.1007/3-540-36545-1.