# Numerical analysis of the neutron multigroup SP\_N equations

Erell Jamelot, François Madiot

<figure><img src="image_1.png" /></figure>

## To cite this version:

Erell Jamelot, François Madiot. Numerical analysis of the neutron multigroup SP\_N equations. Comptes Rendus. Mathématique, 2021, 359 (5), pp.533-545. 10.5802/crmath.189. cea-04490604

HAL Id: cea-04490604

https://cea.hal.science/cea-04490604

Submitted on 4 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Numerical analysis of the neutron multigroup
$SP_N$ equations

Erell Jamelot *¹ and François Madiot †²

¹ DEN-Service de Thermohydraulique et de Mécanique des Fluides-STMF, LMSF, CEA, Université Paris-Saclay, F-91191, Gif-sur-Yvette, France.

² DEN-Service d'Etudes des Réacteurs et de Mathématiques Appliquées (SERMA), CEA, Université Paris-Saclay, F-91191, Gif-sur-Yvette, France.

July 20, 2020

Abstract

The multigroup neutron $SP_N$ equations, which are an approximation of the neutron transport equation, are used to model nuclear reactor cores. In their steady state, these equations can be written as a source problem or an eigenvalue problem. We study the resolution of those two problems with an $H^1$-conforming finite element method and a Discontinuous Galerkin method, namely the Symmetric Interior Penalty Galerkin method.

# 1 Introduction

The neutron transport equation describes the neutron flux density in a reactor core. It depends on 7 variables: 3 for the space, 2 for the motion direction, 1 for the energy (or the speed), and 1 for the time.

The energy variable is discretized using the multigroup theory [9]. In this method, the entire range of neutron energies is divided into *G* intervals, called energy groups. In each energy group, the neutron flux density is lumped and all parameters are averaged. We denote by $\mathcal{I}_G := \{1, \dots, G\}$, the set of energy group indices.

Concerning the motion direction, the $P_N$ transport equations are obtained by developing the neutron flux on the spherical harmonics from order 0 to order

<footnote>*erell.jamelot@cea.fr</footnote>

<footnote>†francois.madiot@cea.fr</footnote>
Numerical analysis of the neutron multigroup
$SP_N$ equations

Erell Jamelot *¹ and François Madiot †²

¹ DEN-Service de Thermohydraulique et de Mécanique des Fluides-STMF, LMSF, CEA, Université Paris-Saclay, F-91191, Gif-sur-Yvette, France.

² DEN-Service d'Etudes des Réacteurs et de Mathématiques Appliquées (SERMA), CEA, Université Paris-Saclay, F-91191, Gif-sur-Yvette, France.

July 20, 2020

Abstract

The multigroup neutron $SP_N$ equations, which are an approximation of the neutron transport equation, are used to model nuclear reactor cores. In their steady state, these equations can be written as a source problem or an eigenvalue problem. We study the resolution of those two problems with an $H^1$-conforming finite element method and a Discontinuous Galerkin method, namely the Symmetric Interior Penalty Galerkin method.

# 1 Introduction

The neutron transport equation describes the neutron flux density in a reactor core. It depends on 7 variables: 3 for the space, 2 for the motion direction, 1 for the energy (or the speed), and 1 for the time.

The energy variable is discretized using the multigroup theory [9]. In this method, the entire range of neutron energies is divided into *G* intervals, called energy groups. In each energy group, the neutron flux density is lumped and all parameters are averaged. We denote by $\mathcal{I}_G := \{1, \dots, G\}$, the set of energy group indices.

Concerning the motion direction, the $P_N$ transport equations are obtained by developing the neutron flux on the spherical harmonics from order 0 to order

<footnote>*erell.jamelot@cea.fr</footnote>

<footnote>†francois.madiot@cea.fr</footnote>
N. This approach is very time-consuming. The simplified $P_N(SP_N)$ transport theory [11] was developed to address this issue. The two fundamental hypotheses to obtain the $SP_N$ equations are that locally, the angular flux has a planar symmetry; and that the axis system evolves slowly. The neutron flux and the scattering cross sections are then developed on the Legendre polynomials. From a mathematical point of view, $SP_N$ equations correspond to tensorized $1D P_N$ transport equations, so that some couplings are missing. Consequently, the $SP_N$ equations do not converge to transport equations. Nevertheless, they are commonly used by physicists since their resolution is cheap in terms of computational cost. The order $N$ is odd, and the number of $SP_N$ odd (resp. even) moments is $\tilde{N} := \frac{N+1}{2}$. We will denote by $\mathcal{I}_e$ (resp. $\mathcal{I}_o$) the subset of even (resp. odd) integers of the integer set $\{0, \dots, N\}$.

Finally, the (motion direction and energy) discretization of the neutron flux is such that there are $\tilde{N} \times G$ even and odd moments of the neutron flux.

We will denote by $\phi = ((\phi_m^g)_{m \in \mathcal{I}_e})_{g \in \mathcal{I}_G} \in \mathbb{R}^{\tilde{N} \times G}$ the set of functions containing, for all energy group $g$, the even moments of the neutron flux.

Likewise, we will denote by $\mathbf{p} = (( (p_{x,m}^g)_{m \in \mathcal{I}_o})_{g \in \mathcal{I}_G})_{x=1}^d \in (\mathbb{R}^{\tilde{N} \times G})^d$ the set of functions containing the odd moments of the neutron flux.

Note that while modelling the core of a pressurized water reactor, the number of groups if such that $2 \le G \lesssim 30$, physicists usually choose $N = 1$ or 3, more rarely $N = 5$.

## 2 Setting of the model

The reactor core is modelled by a bounded, connected and open subset $\mathcal{R}$ of $\mathbb{R}^d$, $d = 1, 2, 3$, having a Lipschitz boundary which is piecewise regular. The coefficients are piecewise regular, so that we split $\mathcal{R}$ into $\tilde{N}$ open disjoint parts $(\mathcal{R}_i)_{\tilde{N}=1}$ with Lipschitz, piecewise regular boundaries: $\bar{\mathcal{R}} = \cup_{i=1}^{\tilde{N}} \bar{\mathcal{R}}_i$. For this reason, we will use the following space of piecewise regular functions:

$$ \mathcal{P}W^{1,\infty}(\mathcal{R}) = \{ D \in L^{\infty}(\mathcal{R}) \mid \vec{\nabla} D|_{\mathcal{R}_i} \in (L^{\infty}(\mathcal{R}_i))^d, i = 1, \dots, \tilde{N} \}. $$

For a set of functions $\psi = (\phi_m^g)_{m,g} \in \mathbb{R}^{\tilde{N} \times G}$, we make the following abuse of notation: $\vec{\nabla} \psi = ((\partial_x \psi_m^g)_{m,g})_{x=1}^d \in (\mathbb{R}^{\tilde{N} \times G})^d$.

For a set of vector valued functions $\mathbf{q} = (( (q_{x,m}^g)_{m,g})_{x=1}^d) \in (\mathbb{R}^{\tilde{N} \times G})^d$, we make the following abuse of notation:

$$ \mathrm{div} \, \mathbf{q} = (\mathrm{div} \, ((q_{x,m}^g)_{x=1}^d))_{m,g}, \quad \mathbf{q} \cdot \mathbf{p} = \left( \sum_{x=1}^d q_{x,m}^g p_{x,m}^g \right)_{m,g} \in \mathbb{R}^{\tilde{N} \times G}. $$

Let us use these notations: for $E \subset \mathbb{R}^d$, $L(E) = L^2(E)$; $L := L^2(\mathcal{R})$; $V := H_0^1(\mathcal{R})$; $V' := H^{-1}(\mathcal{R})$ its dual and $Q := H(\mathrm{div}, \mathcal{R})$. For $W = L(E)$, $L$, $V$ or
N. This approach is very time-consuming. The simplified $P_N(SP_N)$ transport theory [11] was developed to address this issue. The two fundamental hypotheses to obtain the $SP_N$ equations are that locally, the angular flux has a planar symmetry; and that the axis system evolves slowly. The neutron flux and the scattering cross sections are then developed on the Legendre polynomials. From a mathematical point of view, $SP_N$ equations correspond to tensorized $1D P_N$ transport equations, so that some couplings are missing. Consequently, the $SP_N$ equations do not converge to transport equations. Nevertheless, they are commonly used by physicists since their resolution is cheap in terms of computational cost. The order N is odd, and the number of $SP_N$ odd (resp. even) moments is $\tilde{N} := \frac{N+1}{2}$. We will denote by $\mathcal{I}_e$ (resp. $\mathcal{I}_o$) the subset of even (resp. odd) integers of the integer set $\{0, \dots, N\}$.

Finally, the (motion direction and energy) discretization of the neutron flux is such that there are $\tilde{N} \times G$ even and odd moments of the neutron flux.

We will denote by $\phi = ((\phi_m^g)_{m \in \mathcal{I}_e})_{g \in \mathcal{I}_G} \in \mathbb{R}^{\tilde{N} \times G}$ the set of functions containing, for all energy group $g$, the even moments of the neutron flux.

Likewise, we will denote by $\mathbf{p} = (( (p_{x,m}^g)_{m \in \mathcal{I}_o})_{g \in \mathcal{I}_G})_{x=1}^d \in (\mathbb{R}^{\tilde{N} \times G})^d$ the set of functions containing the odd moments of the neutron flux.

Note that while modelling the core of a pressurized water reactor, the number of groups if such that $2 \le G \lesssim 30$, physicists usually choose $N = 1$ or 3, more rarely $N = 5$.

## 2 Setting of the model

The reactor core is modelled by a bounded, connected and open subset $\mathcal{R}$ of $\mathbb{R}^d$, $d = 1, 2, 3$, having a Lipschitz boundary which is piecewise regular. The coefficients are piecewise regular, so that we split $\mathcal{R}$ into $\tilde{N}$ open disjoint parts $(\mathcal{R}_i)_{\tilde{N}=1}$ with Lipschitz, piecewise regular boundaries: $\bar{\mathcal{R}} = \cup_{i=1}^{\tilde{N}} \bar{\mathcal{R}}_i$. For this reason, we will use the following space of piecewise regular functions:

$$ \mathcal{P}W^{1,\infty}(\mathcal{R}) = \{ D \in L^\infty(\mathcal{R}) \mid \vec{\nabla} D|_{\mathcal{R}_i} \in (L^\infty(\mathcal{R}_i))^d, i = 1, \dots, \tilde{N} \}. $$

For a set of functions $\psi = (\phi_m^g)_{m,g} \in \mathbb{R}^{\tilde{N} \times G}$, we make the following abuse of notation: $\vec{\nabla} \psi = ((\partial_x \psi_m^g)_{m,g})_{x=1}^d \in (\mathbb{R}^{\tilde{N} \times G})^d$.

For a set of vector valued functions $\mathbf{q} = (( (q_{x,m}^g)_{m,g})_{x=1}^d) \in (\mathbb{R}^{\tilde{N} \times G})^d$, we make the following abuse of notation:

$$ \mathrm{div} \, \mathbf{q} = (\mathrm{div} \, ((q_{x,m}^g)_{x=1}^d))_{m,g}, \quad \mathbf{q} \cdot \mathbf{p} = \left( \sum_{x=1}^d q_{x,m}^g p_{x,m}^g \right)_{m,g} \in \mathbb{R}^{\tilde{N} \times G}. $$

Let us use these notations: for $E \subset \mathbb{R}^d$, $L(E) = L^2(E)$; $L := L^2(\mathcal{R})$; $V := H_0^1(\mathcal{R})$; $V' := H^{-1}(\mathcal{R})$ its dual and $Q := H(\mathrm{div}, \mathcal{R})$. For $W = L(E)$, $L$, $V$ or
$Q$ we define the product space $\underline{W} := W^{\hat{N}\times G}$ endowed with the following scalar product and associated norm:

$$ (\mathbf{u}, \mathbf{v})_{\underline{W}} = \sum_{g \in \mathcal{I}_G} \sum_{m \in \mathcal{I}_{e,o}} (\mathbf{u}_m^g, \mathbf{v}_m^g)_{\underline{W}}, \quad \|\mathbf{u}\|_{\underline{W}}^2 = \sum_{g \in \mathcal{I}_G} \sum_{m \in \mathcal{I}_{e,o}} \|\mathbf{u}_m^g\|_{\underline{W}}^2. \quad (1) $$

We also set $\underline{V}' := (V')^{\hat{N}\times G}$, $\underline{L}(E) = (\underline{L}(E))^d$ and $\underline{L}^p(\cdot) = (L^p(\cdot))^{\hat{N}\times G}$.

Let $\underline{q} \in (\mathbb{R}^{\hat{N}\times G})^d$ and $\underline{M} \in (\mathbb{R}^{\hat{N}\times \hat{N}})^{G\times G}$. We set $\underline{q}_x = (q_{x,m}^g)_{m,g}$ and we use the notation $\underline{M}\underline{q} = (\underline{M}\underline{q}_x)_{x=1}^d$.

Given a source term $S_f \in \underline{L}$, the multigroup SPN equations with zero-flux boundary conditions<fnref n="1" /> read as coupled diffusion-like equations set in a mixed formulation:

$$ \text{Solve in } (\phi, \mathbf{p}) \in \underline{V} \times \underline{Q} \mid \left\{ \begin{array}{l} \mathbb{T}_o \mathbf{p} + \vec{\nabla}(\mathbb{H}\phi) = 0, \\ t\mathbb{H}\text{div } \mathbf{p} + \mathbb{T}_e \phi = S_f. \end{array} \right. \quad (2) $$

When $S_f$ depends on $\phi$, the steady state multigroup SPN equations read as the following generalized eigenproblem:

$$ \text{Solve in } (\lambda, \phi, \mathbf{p}) \in \mathbb{R}^* \times \underline{V} \times \underline{Q} \mid \left\{ \begin{array}{l} \mathbb{T}_o \mathbf{p} + \vec{\nabla}(\mathbb{H}\phi) = 0, \\ t\mathbb{H}\text{div } \mathbf{p} + \mathbb{T}_e \phi = \lambda^{-1}\mathbb{M}_f\phi. \end{array} \right. \quad (3) $$

The physical solution to Problem (3) corresponds to the eigenfunction associated to the smallest eigenvalue, which in addition is simple [7]. In neutronics, the multiplication factor $k_{eff} = \max_{\lambda} \lambda$ characterizes the physical state of the core reactor: if $k_{eff} = 1$: the nuclear chain reaction is self-sustaining; if $k_{eff} > 1$: the chain reaction is diverging; if $k_{eff} < 1$: the chain reaction vanishes.

The matrices $\mathbb{H}$, $\mathbb{T}_e$, $\mathbb{T}_o$, $\mathbb{M}_f \in (\mathbb{R}^{\hat{N}\times \hat{N}})^{G\times G}$ are such that $\forall (g, g') \in \mathcal{I}_G \times \mathcal{I}_G$ ($\delta_{.,.}$ is the Kronecker symbol):

* $(\mathbb{H})_{g,g'} = \delta_{g,g'} \hat{\mathbb{H}} \in \mathbb{R}^{\hat{N}\times \hat{N}}$, with $\forall (i, j) \in \{1, \dots, \hat{N}\}^2$, $\hat{\mathbb{H}}_{i,j} = \delta_{i,j} + \delta_{i,j-1}$.

* $(\mathbb{T}_e)_{g,g} := \mathbb{T}_e^g \in \mathbb{R}^{\hat{N}\times \hat{N}}$ denotes the even removal matrix, such that:

$$ \mathbb{T}_e^g = \text{diag} (t_0 \sigma_{r,0}^g, t_2 \sigma_{r,2}^g, \dots), $$

$(\mathbb{T}_o)_{g,g} := \mathbb{T}_o^g \in \mathbb{R}^{\hat{N}\times \hat{N}}$ denotes the odd removal matrix, such that:

$$ \mathbb{T}_o^g = \text{diag} (t_1 \sigma_{r,1}^g, t_3 \sigma_{r,3}^g, \dots), $$

where $\forall m \in \mathcal{I}_{e,o}$, $\sigma_{r,m}^g := \sigma_t^g - \sigma_{s,m}^{g \to g}$, and $\forall m > 0$, $t_m > 0$.

The coefficient $\sigma_t^g$ is the macroscopic total cross section of energy group $g$, and the coefficients $\sigma_{s,m}^{g \to g}$ denote the $P_N$ moments of the macroscopic self scattering cross sections from energy group $g$ to itself.

<footnote n="1">¹ie: for $1 \le g \le G$, $m \in \mathcal{I}_e$, $(\phi_m^g)|_{\partial \mathcal{R}} = 0$.</footnote>
$Q$ we define the product space $\underline{W} := W^{\hat{N}\times G}$ endowed with the following scalar product and associated norm:

$$ (\mathbf{u}, \mathbf{v})_{\underline{W}} = \sum_{g \in \mathcal{I}_G} \sum_{m \in \mathcal{I}_{e,o}} (\mathbf{u}_m^g, \mathbf{v}_m^g)_{\underline{W}}, \quad \|\mathbf{u}\|_{\underline{W}}^2 = \sum_{g \in \mathcal{I}_G} \sum_{m \in \mathcal{I}_{e,o}} \|\mathbf{u}_m^g\|_{\underline{W}}^2. \quad (1) $$

We also set $\underline{V}' := (V')^{\hat{N}\times G}$, $\underline{L}(E) = (\underline{L}(E))^d$ and $\underline{L}^p(\cdot) = (L^p(\cdot))^{\hat{N}\times G}$.

Let $\underline{q} \in (\mathbb{R}^{\hat{N}\times G})^d$ and $\underline{M} \in (\mathbb{R}^{\hat{N}\times \hat{N}})^{G\times G}$. We set $\underline{q}_x = (q_{x,m}^g)_{m,g}$ and we use the notation $\underline{M}\underline{q} = (\underline{M}\underline{q}_x)_{x=1}^d$.

Given a source term $S_f \in \underline{L}$, the multigroup SPN equations with zero-flux boundary conditions<fnref n="1" /> read as coupled diffusion-like equations set in a mixed formulation:

$$ \text{Solve in } (\phi, \mathbf{p}) \in \underline{V} \times \underline{Q} \mid \left\{ \begin{array}{l} \mathbb{T}_o \mathbf{p} + \vec{\nabla} (\mathbb{H} \phi) = 0, \\ t \mathbb{H} \text{div } \mathbf{p} + \mathbb{T}_e \phi = S_f. \end{array} \right. \quad (2) $$

When $S_f$ depends on $\phi$, the steady state multigroup SPN equations read as the following generalized eigenproblem:

$$ \text{Solve in } (\lambda, \phi, \mathbf{p}) \in \mathbb{R}^* \times \underline{V} \times \underline{Q} \mid \left\{ \begin{array}{l} \mathbb{T}_o \mathbf{p} + \vec{\nabla} (\mathbb{H} \phi) = 0, \\ t \mathbb{H} \text{div } \mathbf{p} + \mathbb{T}_e \phi = \lambda^{-1} \mathbb{M}_f \phi. \end{array} \right. \quad (3) $$

The physical solution to Problem (3) corresponds to the eigenfunction associated to the smallest eigenvalue, which in addition is simple [7]. In neutronics, the multiplication factor $k_{eff} = \max_{\lambda} \lambda$ characterizes the physical state of the core reactor: if $k_{eff} = 1$: the nuclear chain reaction is self-sustaining; if $k_{eff} > 1$: the chain reaction is diverging; if $k_{eff} < 1$: the chain reaction vanishes.

The matrices $\mathbb{H}$, $\mathbb{T}_e$, $\mathbb{T}_o$, $\mathbb{M}_f \in (\mathbb{R}^{\hat{N}\times \hat{N}})^{G\times G}$ are such that $\forall (g, g') \in \mathcal{I}_G \times \mathcal{I}_G$ ($\delta_{.,.}$ is the Kronecker symbol):

* $(\mathbb{H})_{g,g'} = \delta_{g,g'} \hat{\mathbb{H}} \in \mathbb{R}^{\hat{N}\times \hat{N}}$, with $\forall (i, j) \in \{1, \dots, \hat{N}\}^2$, $\hat{\mathbb{H}}_{i,j} = \delta_{i,j} + \delta_{i,j-1}$.

* $(\mathbb{T}_e)_{g,g} := \mathbb{T}_e^g \in \mathbb{R}^{\hat{N}\times \hat{N}}$ denotes the even removal matrix, such that:

$$ \mathbb{T}_e^g = \text{diag} (t_0 \sigma_{r,0}^g, t_2 \sigma_{r,2}^g, \dots), $$

$(\mathbb{T}_o)_{g,g} := \mathbb{T}_o^g \in \mathbb{R}^{\hat{N}\times \hat{N}}$ denotes the odd removal matrix, such that:

$$ \mathbb{T}_o^g = \text{diag} (t_1 \sigma_{r,1}^g, t_3 \sigma_{r,3}^g, \dots), $$

where $\forall m \in \mathcal{I}_{e,o}, \sigma_{r,m}^g := \sigma_t^g - \sigma_{s,m}^{g \to g}$, and $\forall m > 0, t_m > 0$.

The coefficient $\sigma_t^g$ is the macroscopic total cross section of energy group $g$, and the coefficients $\sigma_{s,m}^{g \to g}$ denote the $P_N$ moments of the macroscopic self scattering cross sections from energy group $g$ to itself.

<footnote n="1">¹ie: for $1 \le g \le G$, $m \in \mathcal{I}_e$, $(\phi_m^g)|_{\partial \mathcal{R}} = 0$.</footnote>
• For $g' \neq g$:

($\mathbb{T}_e)_{g,g'} := -\mathbb{S}_e^{g' \to g} \in \mathbb{R}^{\hat{N} \times \hat{N}}$ denotes the even scattering matrix, such that:

$$ \mathbb{S}_e^{g' \to g} = \mathrm{diag} \left( t_0 \sigma_{s,0}^{g' \to g}, t_2 \sigma_{s,2}^{g' \to g}, \dots \right), $$

($\mathbb{T}_o)_{g,g'} := -\mathbb{S}_o^{g' \to g} \in \mathbb{R}^{\hat{N} \times \hat{N}}$ denotes the odd scattering matrix, such that:

$$ \mathbb{S}_o^{g' \to g} = \mathrm{diag} \left( t_1 \sigma_{s,1}^{g' \to g}, t_3 \sigma_{s,3}^{g' \to g}, \dots \right), $$

where $\sigma_{s,m}^{g' \to g}$ are the $P_N$ moments of the macroscopic scattering cross sections from energy group $g'$ to energy group $g$.

• ($\mathbb{M}_f)_{g,g'} := \chi^g \mathbb{M}_f^{g'} \in \mathbb{R}^{\hat{N} \times \hat{N}}$ is such that $\mathbb{M}_f^{g'} \phi^{g'} = {}^t(\nu \sigma_f^{g'} \phi_0^{g'}, 0, \dots)$ where the coefficient $\nu \sigma_f^{g'}$ is the product of the number of neutrons emitted per fission times the macroscopic fission cross section; and the coefficient $\chi_g$ is the fission spectrum of energy group $g$.

The coefficients of the matrices $\mathbb{T}_{e,o}$, $\mathbb{M}_f$ are supposed to be such that:

$$ \left\{ \begin{array}{ll} (0) & \forall g, g' \in \mathcal{I}_G, \forall m \in \mathcal{I}_{e,o}: \\ & (\sigma_{r,m}^g, \sigma_{s,m}^{g' \to g}, \nu \sigma_f^g) \in \mathcal{P}W^{1,\infty}(\mathbb{R}) \times L^\infty(\mathbb{R}) \times L^\infty(\mathbb{R}). \\ (i) & \exists (\sigma_{r,(e,o)})_*, (\sigma^{r,(e,o)})^* > 0 \mid \forall g \in \mathcal{I}_G, \forall m \in \mathcal{I}_{e,o}: \\ & (\sigma_{r,(e,o)})_* \le t_m \sigma_{r,m}^g \le (\sigma^{r,(e,o)})^* \text{ a.e. in } \mathbb{R}. \\ (ii) & \exists (\nu \sigma_f)^* > 0 \mid \forall g \in \mathcal{I}_G, 0 \le \nu \sigma_f^g \le (\nu \sigma_f)^* \text{ a.e. in } \mathbb{R} \text{ and } \exists g' \mid \nu \sigma_f^{g'} \ne 0. \\ (iii) & \exists 0 < \varepsilon < \frac{1}{G-1} \mid \forall m \in \mathcal{I}_{e,o}, \forall g, g' \in \mathcal{I}_G, g' \ne g, \\ & |\sigma_{s,m}^{g \to g'}| \le \varepsilon \sigma_{r,m}^g \text{ a.e. in } \mathbb{R}. \end{array} \right. \quad (4) $$

It happens that the coefficient $\nu \sigma_f^g$ vanishes in some regions.

Hypothesis 4—(iii) is valid while modelling the core of a pressurized water reactor: the scattering cross-sections are weaker than the removal cross-sections of an order $0 < \varepsilon << 1$. Thus, the matrices ${}^t\mathbb{T}_{e,o}$ are strictly diagonally dominant matrices: they are invertible.

Let us set $\mathbb{D} = {}^t\mathbb{H}\mathbb{T}_o^{-1}\mathbb{H}$.

Problem 2 can be written as a set of coupled primal diffusion-like equations with single unknown $\phi \in \mathcal{V}$:

$$ \text{Solve in } \phi \in \mathcal{V} \mid -\mathrm{div} \left( \mathbb{D} \nabla \phi \right) + \mathbb{T}_e \phi = S_f. \quad (5) $$

The variational formulation of (5) writes:

$$ \text{Solve in } \phi \in \mathcal{V} \mid \forall \psi \in \mathcal{V}: c(\phi, \psi) = \ell(\psi), \quad (6) $$

where: $\left\{ \begin{array}{rcl} c: \mathcal{V} \times \mathcal{V} & \to & \mathbb{R} \\ c(\phi, \psi) & = & (\mathbb{D} \nabla \phi, \nabla \psi)_\mathcal{L} + (\mathbb{T}_e \phi, \psi)_\mathcal{L}, \end{array} \right.$ and $\left\{ \begin{array}{rcl} \ell: \mathcal{V} & \to & \mathbb{R} \\ \ell(\psi) & = & (S_f, \psi)_\mathcal{L}. \end{array} \right.$
• For $g' \neq g$:

($T_e$)$_{g,g'}$ := $-\mathbb{S}_e^{g' \to g} \in \mathbb{R}^{\hat{N} \times \hat{N}}$ denotes the even scattering matrix, such that:

$$ \mathbb{S}_e^{g' \to g} = \mathrm{diag} \left( t_0 \sigma_{s,0}^{g' \to g}, t_2 \sigma_{s,2}^{g' \to g}, \dots \right), $$

($T_o$)$_{g,g'}$ := $-\mathbb{S}_o^{g' \to g} \in \mathbb{R}^{\hat{N} \times \hat{N}}$ denotes the odd scattering matrix, such that:

$$ \mathbb{S}_o^{g' \to g} = \mathrm{diag} \left( t_1 \sigma_{s,1}^{g' \to g}, t_3 \sigma_{s,3}^{g' \to g}, \dots \right), $$

where $\sigma_{s,m}^{g' \to g}$ are the $P_N$ moments of the macroscopic scattering cross sections from energy group $g'$ to energy group $g$.

• ($M_f$)$_{g,g'}$ := $\chi^g M_f^{g'} \in \mathbb{R}^{\hat{N} \times \hat{N}}$ is such that $M_f^{g'} \phi^{g'} = {}^t(\nu \sigma_f^{g'} \phi_0^{g'}, 0, \dots)$ where the coefficient $\nu \sigma_f^{g'}$ is the product of the number of neutrons emitted per fission times the macroscopic fission cross section; and the coefficient $\chi_g$ is the fission spectrum of energy group $g$.

The coefficients of the matrices $T_{e,o}$, $M_f$ are supposed to be such that:

$$ \left\{ \begin{array}{ll} (0) & \forall g, g' \in \mathcal{I}_G, \forall m \in \mathcal{I}_{e,o}: \\ & (\sigma_{r,m}^g, \sigma_{s,m}^{g' \to g}, \nu \sigma_f^g) \in \mathcal{P}W^{1,\infty}(\mathcal{R}) \times L^\infty(\mathcal{R}) \times L^\infty(\mathcal{R}). \\ (i) & \exists (\sigma_{r,(e,o)})_*, (\sigma^{r,(e,o)})^* > 0 \mid \forall g \in \mathcal{I}_G, \forall m \in \mathcal{I}_{e,o}: \\ & (\sigma_{r,(e,o)})_* \le t_m \sigma_{r,m}^g \le (\sigma^{r,(e,o)})^* \text{ a.e. in } \mathcal{R}. \\ (ii) & \exists (\nu \sigma_f)^* > 0 \mid \forall g \in \mathcal{I}_G, 0 \le \nu \sigma_f^g \le (\nu \sigma_f)^* \text{ a.e. in } \mathcal{R} \text{ and } \exists g' \mid \nu \sigma_f^{g'} \ne 0. \\ (iii) & \exists 0 < \varepsilon < \frac{1}{G-1} \mid \forall m \in \mathcal{I}_{e,o}, \forall g, g' \in \mathcal{I}_G, g' \ne g, \\ & |\sigma_{s,m}^{g \to g'}| \le \varepsilon \sigma_{r,m}^g \text{ a.e. in } \mathcal{R}. \end{array} \right. \quad (4) $$

It happens that the coefficient $\nu \sigma_f^g$ vanishes in some regions.

Hypothesis 4—(iii) is valid while modelling the core of a pressurized water reactor: the scattering cross-sections are weaker than the removal cross-sections of an order $0 < \varepsilon < \ll 1$. Thus, the matrices $^tT_{e,o}$ are strictly diagonally dominant matrices: they are invertible.

Let us set $\mathbb{D} = {}^t\mathbb{H}\mathbb{T}_o^{-1}\mathbb{H}$.

Problem 2 can be written as a set of coupled primal diffusion-like equations with single unknown $\phi \in \mathcal{V}$:

$$ \text{Solve in } \phi \in \mathcal{V} \mid -\operatorname{div} (\mathbb{D} \nabla \phi) + T_e \phi = S_f. \quad (5) $$

The variational formulation of (5) writes:

$$ \text{Solve in } \phi \in \mathcal{V} \mid \forall \psi \in \mathcal{V}: c(\phi, \psi) = \ell(\psi), \quad (6) $$

where: $\begin{cases} c: \mathcal{V} \times \mathcal{V} \to \mathbb{R} \\ c(\phi, \psi) = (\mathbb{D} \nabla \phi, \nabla \psi)_L + (T_e \phi, \psi)_L \end{cases}$, and $\begin{cases} \ell: \mathcal{V} \to \mathbb{R} \\ \ell(\psi) = (S_f, \psi)_L \end{cases}$.
**Theorem 1** Suppose that $\mathbb{D}$ is positive definite. For a given source term $S_f \in \underline{\mathbb{L}}$, it exists a unique $\phi \in \underline{\mathbb{V}}$ that solves Problem 6. In addition, it holds: $\|\phi\|_{\underline{\mathbb{V}}} \lesssim \|S_f\|_{\underline{\mathbb{L}}}$. □

**Proof:** The bilinear form $c$ and the linear form $\ell$ are continuous and under the hypothesis on $\mathbb{D}$, the bilinear form $c$ is coercive: we can apply Lax-Milgram theorem to conclude. □

In the same way, Problem 3 can be written as:

$$ \text{Solve in } (\lambda, \phi) \in \mathbb{R}^* \times \underline{\mathbb{V}}\setminus\{0\} \mid -\text{div} \left( \mathbb{D} \vec{\nabla} \phi \right) + \mathbb{T}_e \phi = \lambda^{-1} \mathbb{M}_f \phi. \quad (7) $$

The variational formulation of (7) writes:

$$ \text{Solve in } (\lambda, \phi) \in \mathbb{R}^* \times \underline{\mathbb{V}}\setminus\{0\} \mid \forall \psi \in \underline{\mathbb{V}}: c(\phi, \psi) = \lambda^{-1} \ell_f(\phi, \psi), \quad (8) $$

$$ \text{where: } \left\{ \begin{array}{l} \ell_f: \underline{\mathbb{L}} \times \underline{\mathbb{L}} \rightarrow \mathbb{R} \\ \ell_f(\phi, \psi) = (\mathbb{M}_f \phi, \psi)_{\underline{\mathbb{L}}} \end{array} \right. . $$

**Theorem 2** Suppose that $\mathbb{D}$ is positive definite. There exists a unique compact operator $T_f: \underline{\mathbb{L}} \rightarrow \underline{\mathbb{L}}$ such that $\forall (\phi, \psi) \in \underline{\mathbb{L}} \times \underline{\mathbb{V}}: c(T_f\phi, \psi) = \ell_f(\phi, \psi)$. □

**Proof:** The bilinear form $c$ is a continuous and under the hypothesis on $\mathbb{D}$, it is coercive onto $\underline{\mathbb{V}} \times \underline{\mathbb{V}}$. The bilinear form $\ell_f$ is a continuous onto $\underline{\mathbb{L}} \times \underline{\mathbb{V}}$. Finally, $\underline{\mathbb{V}}$ is a subset of $\underline{\mathbb{L}}$ with a compact embedding. We can then apply the work of Babuška and Osborn in [2]. □

Thus, the couple $(\phi, \lambda^{-1})$ is a solution to Problem 8 iff the couple $(\phi, \lambda)$ is an eigenpair of operator $T_f$. Moreover, Problem 8 admits a countable number of eigenvalues.

We propose first to derive conditions on the macroscopic cross sections so that Problems 5 and 7 are well-posed. Then we obtain a priori error estimates for a discretization performed with some $H^1$-conforming FEM and a Discontinuous Galerkin method, namely the Symmetric Interior Penalty Galerkin method (SIPG) [8, Chapter 4]. The outline is as follows: in Section 3, we exhibit some conditions so that the matrix $\mathbb{T}_o^{-1}$ and $\mathbb{T}_e$ are positive definite. Then we study the discretization of the source problem (5) in Section 5, and the discretization of the eigenproblem in Section 6. Finally, we perform in Section 7 a numerical study of convergence on a benchmark representative of a nuclear core.

### 3 Properties of $\mathbb{T}_e$ and $\mathbb{T}_o^{-1}$

Consider the diagonal matrix containing the even (resp. odd) removal macroscopic cross sections: $\mathbb{T}_{r,(e,o)} = \text{diag}(\mathbb{T}_{e,o}^1, \dots, \mathbb{T}_{e,o}^G)$. We split $\mathbb{T}_{e,o}$ so that: $\mathbb{T}_{e,o} = \mathbb{T}_{r,(e,o)}(\mathbb{I} - \varepsilon \mathbb{U}_{e,o})$, where $\mathbb{I} \in (\mathbb{R}^{\hat{N} \times \hat{N}})^{G \times G}$ is the identity matrix, and:

$$ \begin{align*} \forall g, g' \in \mathcal{I}_G, g' \neq g, \quad (\mathbb{U}_{e,o})_{g,g'} &= \text{diag} \left( \left( \frac{\sigma_{s,m}^{g' \to g}}{\varepsilon \sigma_{r,m}^g} \right)_{m \in \mathcal{I}_{e,o}} \right) \in \mathbb{R}^{\hat{N} \times \hat{N}}; \\ \forall g \in \mathcal{I}_G, \quad (\mathbb{U}_{e,o})_{g,g} &= 0 \in \mathbb{R}^{\hat{N} \times \hat{N}}. \end{align*} $$
**Theorem 1** Suppose that $\mathbb{D}$ is positive definite. For a given source term $S_f \in \underline{\mathbb{L}}$, it exists a unique $\phi \in \underline{\mathbb{V}}$ that solves Problem 6. In addition, it holds: $\|\phi\|_{\underline{\mathbb{V}}} \lesssim \|S_f\|_{\underline{\mathbb{L}}}$. □

**Proof:** The bilinear form $c$ and the linear form $\ell$ are continuous and under the hypothesis on $\mathbb{D}$, the bilinear form $c$ is coercive: we can apply Lax-Milgram theorem to conclude. □

In the same way, Problem 3 can be written as:

$$ \text{Solve in } (\lambda, \phi) \in \mathbb{R}^* \times \underline{\mathbb{V}}\setminus\{0\} \mid -\text{div} \left( \mathbb{D}\vec{\nabla} \phi \right) + \mathbb{T}_e \phi = \lambda^{-1}\mathbb{M}_f \phi. \quad (7) $$

The variational formulation of (7) writes:

$$ \text{Solve in } (\lambda, \phi) \in \mathbb{R}^* \times \underline{\mathbb{V}}\setminus\{0\} \mid \forall \psi \in \underline{\mathbb{V}}: c(\phi, \psi) = \lambda^{-1}\ell_f(\phi, \psi), \quad (8) $$

$$ \text{where: } \left\{ \begin{array}{l} \ell_f: \underline{\mathbb{L}} \times \underline{\mathbb{L}} \rightarrow \mathbb{R} \\ \ell_f(\phi, \psi) = (\mathbb{M}_f \phi, \psi)_{\underline{\mathbb{L}}} \end{array} \right. . $$

**Theorem 2** Suppose that $\mathbb{D}$ is positive definite. There exists a unique compact operator $T_f: \underline{\mathbb{L}} \rightarrow \underline{\mathbb{L}}$ such that $\forall (\phi, \psi) \in \underline{\mathbb{L}} \times \underline{\mathbb{V}}: c(T_f\phi, \psi) = \ell_f(\phi, \psi)$. □

**Proof:** The bilinear form $c$ is a continuous and under the hypothesis on $\mathbb{D}$, it is coercive onto $\underline{\mathbb{V}} \times \underline{\mathbb{V}}$. The bilinear form $\ell_f$ is a continuous onto $\underline{\mathbb{L}} \times \underline{\mathbb{V}}$. Finally, $\underline{\mathbb{V}}$ is a subset of $\underline{\mathbb{L}}$ with a compact embedding. We can then apply the work of Babuška and Osborn in [2]. □

Thus, the couple $(\phi, \lambda^{-1})$ is a solution to Problem 8 iff the couple $(\phi, \lambda)$ is an eigenpair of operator $T_f$. Moreover, Problem 8 admits a countable number of eigenvalues.

We propose first to derive conditions on the macroscopic cross sections so that Problems 5 and 7 are well-posed. Then we obtain a priori error estimates for a discretization performed with some $H^1$-conforming FEM and a Discontinuous Galerkin method, namely the Symmetric Interior Penalty Galerkin method (SIPG) [8, Chapter 4]. The outline is as follows: in Section 3, we exhibit some conditions so that the matrix $\mathbb{T}_o^{-1}$ and $\mathbb{T}_e$ are positive definite. Then we study the discretization of the source problem (5) in Section 5, and the discretization of the eigenproblem in Section 6. Finally, we perform in Section 7 a numerical study of convergence on a benchmark representative of a nuclear core.

# 3 Properties of $\mathbb{T}_e$ and $\mathbb{T}_o^{-1}$

Consider the diagonal matrix containing the even (resp. odd) removal macroscopic cross sections: $\mathbb{T}_{r,(e,o)} = \text{diag}(\mathbb{T}_{e,o}^1, \dots, \mathbb{T}_{e,o}^G)$. We split $\mathbb{T}_{e,o}$ so that: $\mathbb{T}_{e,o} = \mathbb{T}_{r,(e,o)}(\mathbb{I} - \varepsilon\mathbb{U}_{e,o})$, where $\mathbb{I} \in (\mathbb{R}^{\hat{N} \times \hat{N}})^{G \times G}$ is the identity matrix, and:

$$ \forall g, g' \in \mathcal{I}_G, g' \neq g, \quad (\mathbb{U}_{e,o})_{g,g'} = \text{diag} \left( \left( \frac{\sigma_{s,m}^{g' \to g}}{\varepsilon \sigma_{r,m}^g} \right)_{m \in \mathcal{I}_{e,o}} \right) \in \mathbb{R}^{\hat{N} \times \hat{N}}; $$

$$ \forall g \in \mathcal{I}_G, \quad (\mathbb{U}_{e,o})_{g,g} = 0 \in \mathbb{R}^{\hat{N} \times \hat{N}}. $$
We have then: $||U_{e,o}||_2 \lesssim \frac{\alpha_{s,(e,o)}}{\varepsilon}$ where: $\alpha_{s,(e,o)} := (G-1) \max_{m \in \mathcal{I}_{e,o}} \max_{g \neq g' \in \mathcal{I}_G} \sup_{\vec{x} \in \mathbb{R}} \frac{|\sigma_{s,m}^{g' \to g}(\vec{x})|}{\sigma_{r,m}^g(\vec{x})}.$

Let us set $\alpha_{r,(e,o)} = \frac{(\sigma_{r,(e,o)})^*}{(\sigma_{r,(e,o)})_*} > 1$. We have the following properties.

**Property 3** Suppose that $\alpha_{s,e} < \frac{1}{\alpha_{r,e}}$. The matrix $T_e$ is such that:

$$ \forall X \in \mathbb{R}^{\hat{N} \times G} \quad (T_e X | X) \geq \tau_e ||X||_2^2 \quad \text{where } \tau_e = (\sigma_{r,e})_* (1 - \alpha_{r,e} \alpha_{s,e}). \quad (9) $$

**Proof:** We have: $\forall X \in \mathbb{R}^{\hat{N} \times G}$, $(T_e X | X) = (T_{r,e} X | X) - \varepsilon(U_e X | T_{r,e} X)$, so that: $((T_e X | X) \geq (\sigma_{r,e})_* - \varepsilon||U_e||_2 ||T_{r,e}||_2) ||X||_2$, where $||T_{r,e}||_2 \leq (\sigma_{r,e})^*$. $\square$

**Property 4** Suppose that $\alpha_{s,o} < \frac{1}{\alpha_{r,o} + 1}$, the matrix $T_o^{-1}$ is such that:

$$ \forall X \in \mathbb{R}^{\hat{N} \times G} \quad (T_o^{-1} X | X) \geq \tau_o ||X||_2^2 \quad \text{where } \tau_o = \frac{1}{(\sigma_{r,o})^*} \left( 1 - \frac{\alpha_{r,o} \alpha_{s,o}}{1 - \alpha_{s,o}} \right). \quad (10) $$

**Proof:** The Taylor expansion of $T_o^{-1}$ writes: $T_o^{-1} = (I + \sum_{l>0} \varepsilon^l U_o^l) T_{r,o}^{-1}$.

We get that $\forall X \in \mathbb{R}^{\hat{N} \times G}$:

$$ \begin{aligned} (T_o^{-1} X | X) &= (T_{r,o}^{-1} X | X) + \sum_{l>0} \varepsilon^l (U_o^l T_{r,o}^{-1} X | X) \\ &\geq \frac{1}{(\sigma_{r,o})^*} \left( 1 - \alpha_{r,o} \sum_{l>0} \varepsilon^l ||U_o||_2^l \right) ||X||_2^2, \\ &\geq \frac{1}{(\sigma_{r,o})^*} \left( 1 - \alpha_{r,o} \frac{\varepsilon ||U_o||_2}{1 - \varepsilon ||U_o||_2} \right) ||X||_2^2, \\ &\geq \frac{1}{(\sigma_{r,o})^*} \left( 1 - \frac{\alpha_{r,o} \alpha_{s,o}}{1 - \alpha_{s,o}} \right) ||X||_2^2. \end{aligned} $$

<figure><img src="image_1.png" /></figure>

Under assumptions of Properties 3 and 4 the matrices $T_e$ and $T_o^{-1}$ are positive definite. Moreover, one can show that $||H \nabla \phi||_L \gtrsim ||\nabla \phi||_L$ [12]. We infer that the matrix $D$ is positive definite and that there exists a constant $C_D > 0$ such that for all $\xi \in \mathbb{R}^{\hat{N} \times G}$,

$$ (D\xi|D\xi) \leq C_D ||\xi||_2^2. \quad (11) $$

From now on, we suppose that this property holds.
We have then: $\|U_{e,o}\|_2 \lesssim \frac{\alpha_{s,(e,o)}}{\varepsilon}$ where: $\alpha_{s,(e,o)} := (G-1) \max_{m \in \mathcal{I}_{e,o}} \max_{g \neq g' \in \mathcal{I}_G} \sup_{\vec{x} \in \mathbb{R}} \frac{|\sigma_{s,m}^{g' \to g}(\vec{x})|}{\sigma_{r,m}^g(\vec{x})}.$

Let us set $\alpha_{r,(e,o)} = \frac{(\sigma_{r,(e,o)})^*}{(\sigma_{r,(e,o)})_*} > 1$. We have the following properties.

**Property 3** Suppose that $\alpha_{s,e} < \frac{1}{\alpha_{r,e}}$. The matrix $T_e$ is such that:

$$ \forall X \in \mathbb{R}^{\hat{N} \times G} \quad (T_e X | X) \ge \tau_e \|X\|_2^2 \quad \text{where } \tau_e = (\sigma_{r,e})_* (1 - \alpha_{r,e} \alpha_{s,e}). \quad (9) $$

**Proof:** We have: $\forall X \in \mathbb{R}^{\hat{N} \times G}$, $(T_e X | X) = (T_{r,e} X | X) - \varepsilon(U_e X | T_{r,e} X)$, so that: $((T_e X | X) \ge (\sigma_{r,e})_* - \varepsilon\|U_e\|_2 \|T_{r,e}\|_2) \|X\|_2$, where $\|T_{r,e}\|_2 \le (\sigma_{r,e})^*$. $\square$

**Property 4** Suppose that $\alpha_{s,o} < \frac{1}{\alpha_{r,o} + 1}$, the matrix $T_o^{-1}$ is such that:

$$ \forall X \in \mathbb{R}^{\hat{N} \times G} \quad (T_o^{-1} X | X) \ge \tau_o \|X\|_2^2 \quad \text{where } \tau_o = \frac{1}{(\sigma_{r,o})^*} \left(1 - \frac{\alpha_{r,o} \alpha_{s,o}}{1 - \alpha_{s,o}}\right). \quad (10) $$

**Proof:** The Taylor expansion of $T_o^{-1}$ writes: $T_o^{-1} = \left(\mathbb{I} + \sum_{l>0} \varepsilon^l U_o^l\right) T_{r,o}^{-1}$.

We get that $\forall X \in \mathbb{R}^{\hat{N} \times G}$:

$$ \begin{aligned} (T_o^{-1} X | X) &= (T_{r,o}^{-1} X | X) + \sum_{l>0} \varepsilon^l (U_o^l T_{r,o}^{-1} X | X) \\ &\ge \frac{1}{(\sigma_{r,o})^*} \left(1 - \alpha_{r,o} \sum_{l>0} \varepsilon^l \|U_o\|_2^l\right) \|X\|_2^2, \\ &\ge \frac{1}{(\sigma_{r,o})^*} \left(1 - \alpha_{r,o} \frac{\varepsilon \|U_o\|_2}{1 - \varepsilon \|U_o\|_2}\right) \|X\|_2^2, \\ &\ge \frac{1}{(\sigma_{r,o})^*} \left(1 - \frac{\alpha_{r,o} \alpha_{s,o}}{1 - \alpha_{s,o}}\right) \|X\|_2^2. \end{aligned} $$

$\square$

Under assumptions of Properties 3 and 4 the matrices $T_e$ and $T_o^{-1}$ are positive definite. Moreover, one can show that $\|\mathbb{H} \vec{\nabla} \phi\|_{\mathbb{L}} \gtrsim \|\vec{\nabla} \phi\|_{\mathbb{L}}$ [12]. We infer that the matrix $\mathbb{D}$ is positive definite and that there exists a constant $C_{\mathbb{D}} > 0$ such that for all $\xi \in \mathbb{R}^{\hat{N} \times G}$,

$$ (\mathbb{D} \xi | \mathbb{D} \xi) \le C_{\mathbb{D}} \|\xi\|_2^2. \quad (11) $$

From now on, we suppose that this property holds.
# 4 Discretizations

Let $\mathcal{T}_h$ be a shape-regular mesh of $\mathcal{R}$, with mesh size $h$. We denote by $K$ its elements and $F$ its facets. To simplify the presentation, we assume that the meshes are such that in every element, the cross-sections are regular. We define by $\mathcal{F}_h^i$ the set of interior faces of $\mathcal{T}_h$, $\mathcal{F}_h^b$ the set of boundary facets and $\mathcal{F}_h = \mathcal{F}_h^i \cup \mathcal{F}_h^b$. We denote by $N_\partial$ the maximum number of mesh faces composing the boundary of mesh elements

$$N_\partial := \max_{K \in \mathcal{T}_h} \mathrm{Card}\{F \in \mathcal{F}_h, F \subset \partial K\}.$$

We will first consider an $H^1$-conforming finite element method (FEM). For $k \in \mathbb{N}^*$, $V_h^k \subset V$ and $\underline{V}_h^k \subset \underline{V}$ are the finite dimension spaces defined by:

$$V_h^k = \{v_h \in V, \forall K \in \mathcal{T}_h, v_h|_K \in \mathbb{P}_k\}, \quad \underline{V}_h^k := (V_h^k)^{\hat{N} \times G}.$$

The discrete variational formulation associated to Problem (6) writes:

$$\text{Solve in } \phi_h \in \underline{V}_h^k \mid \forall \psi_h \in \underline{V}_h^k : c(\phi_h, \psi_h) = \ell(\psi_h), \qquad (12)$$

Similarly, the discrete variational formulation associated to Problem (7) writes:

$$\text{Solve in } (\lambda_h, \phi_h) \in \mathbb{R}^* \times \underline{V}_h^k \setminus \{0\} \mid \forall \psi \in \underline{V}_h^k : c(\phi_h, \psi_h) = \lambda_h^{-1} \ell_f(\phi_h, \psi_h). \quad (13)$$

Then, we will consider a non-conforming FEM. We define the broken spaces:

$$V_{\mathrm{NC}} = \{v \in L^2(\mathcal{R}) \mid \forall K \in \mathcal{T}_h, v \in H^1(K)\}, \quad \underline{V}_{\mathrm{NC}} = (V_{\mathrm{NC}})^{\hat{N} \times G}.$$

For $(\phi, \psi) \in \underline{V}_{\mathrm{NC}} \times \underline{V}_{\mathrm{NC}}$, and $\mathbb{T} \in \mathbb{R}^{\hat{N} \times G}$, we set:

$$(\mathbb{D} \vec{\nabla}_h \phi, \vec{\nabla}_h \psi)_{\mathcal{T}_h} = \sum_{K \in \mathcal{T}_h} (\mathbb{D} \vec{\nabla}_h \phi, \vec{\nabla}_h \psi)_{\mathbf{L}(K)}, \quad \text{and} \quad \|\vec{\nabla}_h \psi\|_{\mathcal{T}_h} = (\vec{\nabla}_h \psi, \vec{\nabla}_h \psi)_{\mathcal{T}_h}^{1/2}.$$

For $F \in \mathcal{F}_h^i$ such that $F = \partial K_1 \cap \partial K_2$, we define the average $\{\mathbb{D} \vec{\nabla}_h \psi\}$ and the jump $\llbracket \psi \rrbracket$ as:

$$\begin{aligned} \{\mathbb{D} \vec{\nabla}_h \psi\}_F &= \frac{1}{2} \left( (\mathbb{D}_1 \vec{\nabla}_h \psi_1)|_F + (\mathbb{D}_2 \vec{\nabla}_h \psi_2)|_F \right) \in (\mathbb{R}^{\hat{N} \times G})^d, \\ \llbracket \psi \rrbracket|_F &= \psi_1|_F \mathbf{n}_1 + \psi_2|_F \mathbf{n}_2 \in (\mathbb{R}^{\hat{N} \times G})^d. \end{aligned}$$

where $\mathbf{n}_i$ is the unit outward normal to $K_i$ at face $F$ and $\mathbb{D}_i = \mathbb{D}|_{K_i}$, $\psi_i = \psi|_{K_i}$.

For $F \in \mathcal{F}_h^b$ such that $F \in K$, we set $\{\mathbb{D} \vec{\nabla}_h \psi\}_F = \mathbb{D}|_K \vec{\nabla}_h \psi|_K$ and $\llbracket \psi \rrbracket|_F = (\psi_K)|_F \mathbf{n}$, where $\psi_K = \psi|_K$ and $\mathbf{n}$ is the unit outward normal to $K$ at face $F$.

For $k \in \mathbb{N}^*$, $V_{h,\mathrm{NC}}^k \subset H^1(\mathcal{T}_h)$ and $\underline{V}_{h,\mathrm{NC}}^k$ are the finite dimension spaces defined by:

$$V_{h,\mathrm{NC}}^k = \{v_h \in L^1(\mathcal{R}) ; \forall K \in \mathcal{T}_h, v_h|_K \in \mathbb{P}_k\}, \quad \underline{V}_{h,\mathrm{NC}}^k := (V_{h,\mathrm{NC}}^k)^{\hat{N} \times G}.$$
## 4 Discretizations

Let $\mathcal{T}_h$ be a shape-regular mesh of $\mathcal{R}$, with mesh size $h$. We denote by $K$ its elements and $F$ its facets. To simplify the presentation, we assume that the meshes are such that in every element, the cross-sections are regular. We define by $\mathcal{F}_h^i$ the set of interior faces of $\mathcal{T}_h$, $\mathcal{F}_h^b$ the set of boundary facets and $\mathcal{F}_h = \mathcal{F}_h^i \cup \mathcal{F}_h^b$. We denote by $N_\partial$ the maximum number of mesh faces composing the boundary of mesh elements

$$N_\partial := \max_{K \in \mathcal{T}_h} \mathrm{Card}\{F \in \mathcal{F}_h, F \subset \partial K\}.$$

We will first consider an $H^1$-conforming finite element method (FEM). For $k \in \mathbb{N}^*$, $V_h^k \subset V$ and $\underline{V}_h^k \subset \underline{V}$ are the finite dimension spaces defined by:

$$V_h^k = \{v_h \in V, \forall K \in \mathcal{T}_h, v_h|_K \in \mathbb{P}_k\}, \quad \underline{V}_h^k := (V_h^k)^{\hat{N} \times G}.$$

The discrete variational formulation associated to Problem (6) writes:

$$\text{Solve in } \phi_h \in \underline{V}_h^k \mid \forall \psi_h \in \underline{V}_h^k : c(\phi_h, \psi_h) = \ell(\psi_h), \quad (12)$$

Similarly, the discrete variational formulation associated to Problem (7) writes:

$$\text{Solve in } (\lambda_h, \phi_h) \in \mathbb{R}^* \times \underline{V}_h^k \setminus \{0\} \mid \forall \psi \in \underline{V}_h^k : c(\phi_h, \psi_h) = \lambda_h^{-1} \ell_f(\phi_h, \psi_h). \quad (13)$$

Then, we will consider a non-conforming FEM. We define the broken spaces:

$$V_{\mathrm{NC}} = \{v \in L^2(\mathcal{R}) \mid \forall K \in \mathcal{T}_h, v \in H^1(K)\}, \quad \underline{V}_{\mathrm{NC}} = (V_{\mathrm{NC}})^{\hat{N} \times G}.$$

For $(\phi, \psi) \in \underline{V}_{\mathrm{NC}} \times \underline{V}_{\mathrm{NC}}$, and $\mathbb{T} \in \mathbb{R}^{\hat{N} \times G}$, we set:

$$(\mathbb{D} \vec{\nabla}_h \phi, \vec{\nabla}_h \psi)_{\mathcal{T}_h} = \sum_{K \in \mathcal{T}_h} (\mathbb{D} \vec{\nabla}_h \phi, \vec{\nabla}_h \psi)_{\mathbb{L}(K)}, \quad \text{and} \quad \|\vec{\nabla}_h \psi\|_{\mathcal{T}_h} = (\vec{\nabla}_h \psi, \vec{\nabla}_h \psi)_{\mathcal{T}_h}^{1/2}.$$

For $F \in \mathcal{F}_h^i$ such that $F = \partial K_1 \cap \partial K_2$, we define the average $\{\mathbb{D} \vec{\nabla}_h \psi\}$ and the jump $\llbracket \psi \rrbracket$ as:

$$\{\mathbb{D} \vec{\nabla}_h \psi\}_F = \frac{1}{2} \left( (\mathbb{D}_1 \vec{\nabla}_h \psi_1)|_F + (\mathbb{D}_2 \vec{\nabla}_h \psi_2)|_F \right) \in (\mathbb{R}^{\hat{N} \times G})^d, \\ \llbracket \psi \rrbracket|_F = \psi_1|_F \mathbf{n}_1 + \psi_2|_F \mathbf{n}_2 \in (\mathbb{R}^{\hat{N} \times G})^d.$$

where $\mathbf{n}_i$ is the unit outward normal to $K_i$ at face $F$ and $\mathbb{D}_i = \mathbb{D}|_{K_i}$, $\psi_i = \psi|_{K_i}$.

For $F \in \mathcal{F}_h^b$ such that $F \in K$, we set $\{\mathbb{D} \vec{\nabla}_h \psi\}_F = \mathbb{D}|_K \vec{\nabla}_h \psi|_K$ and $\llbracket \psi \rrbracket|_F = (\psi_K)|_F \mathbf{n}$, where $\psi_K = \psi|_K$ and $\mathbf{n}$ is the unit outward normal to $K$ at face $F$.

For $k \in \mathbb{N}^*$, $V_{h,\mathrm{NC}}^k \subset H^1(\mathcal{T}_h)$ and $\underline{V}_{h,\mathrm{NC}}^k$ are the finite dimension spaces defined by:

$$V_{h,\mathrm{NC}}^k = \{v_h \in L^1(\mathcal{R}) ; \forall K \in \mathcal{T}_h, v_h|_K \in \mathbb{P}_k\}, \quad \underline{V}_{h,\mathrm{NC}}^k := (V_{h,\mathrm{NC}}^k)^{\hat{N} \times G}.$$
For $\phi_h, \psi_h \in \underline{V}_{h,NC}^k$, we set: $$ (\{\mathbb{D}\vec{\nabla}_h\phi_h\}, [[\psi_h]])_{\mathcal{F}_h^i} = \sum_{F \in \mathcal{F}_h^i} (\{\mathbb{D}\vec{\nabla}_h\phi_h\}, [[\psi_h]])_{\underline{\mathcal{L}}(F)} $$

Let us set

$$ c_h(\phi_h, \psi_h) = c_{\mathcal{T}_h}(\phi_h, \psi_h) + c_{\mathcal{F}_h}(\phi_h, \psi_h), \quad (14) $$

with

$$ c_{\mathcal{T}_h}(\phi_h, \psi_h) = (\mathbb{D}\vec{\nabla}_h\phi_h, \vec{\nabla}_h\psi_h)_{\mathcal{T}_h} + (\mathbb{T}_e\phi_h, \psi_h)_{\underline{\mathcal{L}}}, $$

$$ c_{\mathcal{F}_h}(\phi_h, \psi_h) = \sum_{F \in \mathcal{F}_h} \frac{\alpha}{h_F} ([[\phi_h]], [[\psi_h]])_{\underline{\mathcal{L}}(F)} - (\{\mathbb{D}\vec{\nabla}_h\psi_h\}, [[\phi_h]])_{\mathcal{F}_h^i} - (\{\mathbb{D}\vec{\nabla}_h\phi_h\}, [[\psi_h]])_{\mathcal{F}_h^i}, $$

where $\alpha$ is a stabilization parameter.

The Symmetric Interior Penalty Galerkin method (SIPG) associated to Problem (6) writes:

$$ \text{Solve in } \phi_h \in \underline{V}_{h,NC}^k \mid \forall \psi_h \in \underline{V}_{h,NC}^k : c_h(\phi_h, \psi_h) = \ell(\psi_h). \quad (15) $$

Similarly, the SIPG method associated to Problem (8) writes:

$$ \text{Solve in } (\lambda_h, \phi_h) \in \mathbb{R}^* \times \underline{V}_{h,NC}^k \setminus \{0\} \mid \forall \psi_h \in \underline{V}_{h,NC}^k : \\ c_h(\phi_h, \psi_h) = \lambda_h^{-1} \ell_f(\phi_h, \psi_h). \quad (16) $$

# 5 The source problem

## 5.1 Conforming discretization

**Theorem 5** Suppose that there exists $r_{\max}$ in $[0, 1]$ such that $\forall r \in [0, r_{\max}[$, $\phi \in (H^{1+r}(\mathcal{R}))^{\hat{N}\times G}$ ([6], Proposition 1). Let us set $\mu = \min(r_{\max}, k)$. The solution of (12), $\phi_h$ is such that: $\|\phi - \phi_h\|_{\underline{V}} \lesssim h^\mu \|S_f\|_{\underline{L}}$ and $\|\phi - \phi_h\|_{\underline{L}} \lesssim h^{2\mu} \|S_f\|_{\underline{L}}$.

**Proof:** From Céa's lemma and Aubin-Nitsche lemma as detailed in ([10], §2.3).□

## 5.2 SIPG discretization

**Assumption 5.1 (Regularity of exact solution and space $V^*$)** Let us denote by $W^{2,p}(\mathcal{T}_h)$ the broken Sobolev space spanned by those functions $v$ such that for all $K \in \mathcal{T}_h$, $v|_K \in W^{2,p}(K)$. We set $\underline{W}^{2,p}(\mathcal{T}_h) = (W^{2,p}(\mathcal{T}_h))^{\hat{N}\times G}$. We assume that $d \ge 2$ and that there is $2d/(d+2) < p \le 2$ such that, for the exact solution $\phi \in \underline{V}^* := \underline{V} \cap \underline{W}^{2,p}(\mathcal{T}_h)$. This holds for our assumptions on the coefficients, which are piecewise constant with respect to the triangulation [15].

This assumption requires $p > 1$ for $d = 2$ and $p > 6/5$ for $d = 3$. In particular, we observe that, in two space dimensions, $\phi \in \underline{W}^{2,p}(\mathcal{T}_h)$ in polygonal domains. Moreover, using Sobolev embeddings [4, Sect. IX.3], this implies

$$ \phi \in (H^{1+\alpha_p}(\mathcal{R}))^{\hat{N}\times G}, \quad \alpha_p = \frac{d+2}{2} - \frac{d}{p} > 0. $$

We state the following lemma [8, Lemma 1.46, p.27].
For $\phi_h, \psi_h \in \underline{V}_{h,NC}^k$, we set: $\{\mathbb{D}\vec{\nabla}_h\phi_h\}, [[\psi_h]]\}_{\mathcal{F}_h^i} = \sum_{F \in \mathcal{F}_h^i} \{\mathbb{D}\vec{\nabla}_h\phi_h\}, [[\psi_h]]\}_{\underline{\mathcal{L}}(F)}$.

Let us set

$$c_h(\phi_h, \psi_h) = c_{\mathcal{T}_h}(\phi_h, \psi_h) + c_{\mathcal{F}_h}(\phi_h, \psi_h), \quad (14)$$

with

$$c_{\mathcal{T}_h}(\phi_h, \psi_h) = (\mathbb{D}\vec{\nabla}_h\phi_h, \vec{\nabla}_h\psi_h)_{\mathcal{T}_h} + (\mathbb{T}_e\phi_h, \psi_h)_{\underline{\mathcal{L}}},$$

$$c_{\mathcal{F}_h}(\phi_h, \psi_h) = \sum_{F \in \mathcal{F}_h} \frac{\alpha}{h_F} ([[\phi_h]], [[\psi_h]])_{\underline{\mathcal{L}}(F)} - (\{\mathbb{D}\vec{\nabla}_h\psi_h\}, [[\phi_h]])_{\mathcal{F}_h^i} - (\{\mathbb{D}\vec{\nabla}_h\phi_h\}, [[\psi_h]])_{\mathcal{F}_h^i},$$

where $\alpha$ is a stabilization parameter.

The Symmetric Interior Penalty Galerkin method (SIPG) associated to Problem (6) writes:

$$\text{Solve in } \phi_h \in \underline{V}_{h,NC}^k \mid \forall \psi_h \in \underline{V}_{h,NC}^k : c_h(\phi_h, \psi_h) = \ell(\psi_h). \quad (15)$$

Similarly, the SIPG method associated to Problem (8) writes:

$$\text{Solve in } (\lambda_h, \phi_h) \in \mathbb{R}^* \times \underline{V}_{h,NC}^k \setminus \{0\} \mid \forall \psi_h \in \underline{V}_{h,NC}^k : \\ c_h(\phi_h, \psi_h) = \lambda_h^{-1} \ell_f(\phi_h, \psi_h). \quad (16)$$

# 5 The source problem

## 5.1 Conforming discretization

**Theorem 5** Suppose that there exists $r_{\max}$ in $[0, 1]$ such that $\forall r \in [0, r_{\max}[$, $\phi \in (H^{1+r}(\mathcal{R}))^{\hat{N}\times G}$ ([6], Proposition 1). Let us set $\mu = \min(r_{\max}, k)$. The solution of (12), $\phi_h$ is such that: $\|\phi - \phi_h\|_{\underline{V}} \lesssim h^\mu \|S_f\|_{\underline{L}}$ and $\|\phi - \phi_h\|_{\underline{L}} \lesssim h^{2\mu} \|S_f\|_{\underline{L}}$.

**Proof:** From Céa's lemma and Aubin-Nitsche lemma as detailed in ([10], §2.3).□

## 5.2 SIPG discretization

**Assumption 5.1 (Regularity of exact solution and space $V^*$)** Let us denote by $W^{2,p}(\mathcal{T}_h)$ the broken Sobolev space spanned by those functions $v$ such that for all $K \in \mathcal{T}_h$, $v|_K \in W^{2,p}(K)$. We set $\underline{W}^{2,p}(\mathcal{T}_h) = (\underline{W}^{2,p}(\mathcal{T}_h))^{\hat{N}\times G}$. We assume that $d \ge 2$ and that there is $2d/(d+2) < p \le 2$ such that, for the exact solution $\phi \in \underline{V}^* := \underline{V} \cap \underline{W}^{2,p}(\mathcal{T}_h)$. This holds for our assumptions on the coefficients, which are piecewise constant with respect to the triangulation [15].

This assumption requires $p > 1$ for $d = 2$ and $p > 6/5$ for $d = 3$. In particular, we observe that, in two space dimensions, $\phi \in \underline{W}^{2,p}(\mathcal{T}_h)$ in polygonal domains. Moreover, using Sobolev embeddings [4, Sect. IX.3], this implies

$$\phi \in (H^{1+\alpha_p}(\mathcal{R}))^{\hat{N}\times G}, \quad \alpha_p = \frac{d+2}{2} - \frac{d}{p} > 0.$$

We state the following lemma [8, Lemma 1.46, p.27].
**Lemma 6** Suppose that $(\mathcal{T}_h)_h$ is a shape- and contact-regular mesh sequence. Then, we have for all $h > 0$:

$$ \forall \psi_h \in \underline{V}_{h,\text{NC}}^k, \forall K \in \mathcal{T}_h, \forall F \in \partial K, \quad h_K^{1/2} \|\psi_h\|_{\underline{L}^2(F)} \le C_{tr} \|\psi_h\|_{\underline{L}^2(K)}, \quad (17) $$

where $h_K$ is the diameter of element K.

We aim at asserting the discrete coercivity using the following norm:

$$ \forall \psi_h \in \underline{V}_{h,\text{NC}}^k, \quad \|\psi_h\|_{sip}^2 := c_{\mathcal{T}_h}(\psi_h, \psi_h) + \|\psi_h\|_J^2 $$

with the jump semi-norm

$$ \|\psi_h\|_J^2 := \sum_{F \in \mathcal{F}_h} \frac{1}{h_F} \|[[\psi_h]]\|_{\underline{L}(F)}^2. $$

Under assumption (4), there exists $\beta > 0$ we have for all $\psi_h \in \underline{V}_{h,\text{NC}}^k$

$$ c_{\mathcal{T}_h}(\psi_h, \psi_h) \ge \beta \left( \|\vec{\nabla}_h \psi_h\|_{\mathcal{T}_h}^2 + \|\psi_h\|_{\underline{L}}^2 \right), \quad (18) $$

so that

$$ \|\psi_h\|_{sip}^2 \ge \beta \left( \|\vec{\nabla}_h \psi_h\|_{\mathcal{T}_h}^2 + \|\psi_h\|_{\underline{L}}^2 + \|\psi_h\|_J^2 \right). $$

**Lemma 7 (Discrete coercivity)** Let $\underline{\alpha} := C_{tr}^2 N_{\partial} \frac{C_{\mathbb{D}}}{\beta}$ where

* $C_{tr}$ results from the discrete trace inequality (17),

* $N_{\partial}$ is defined in Section 4,

* $C_{\mathbb{D}}$ is defined in (11).

For all $\alpha \ge \underline{\alpha}$, the SIP bilinear form defined by (14) is coercive on $\underline{V}_{h,\text{NC}}^k$ with respect to the $\|\cdot\|_{sip}$-norm, i.e.,

$$ c_h(\psi_h, \psi_h) \ge C_{\alpha} \|\psi_h\|_{sip}^2, $$

$$ \text{with } C_{\alpha} := \left( \alpha - C_{tr}^2 N_{\partial} \frac{C_{\mathbb{D}}}{\beta} \right) \min \left\{ \frac{1}{2}, \beta \left( \alpha + C_{tr}^2 N_{\partial} \frac{C_{\mathbb{D}}}{\beta} \right)^{-1} \right\}. $$

**Proof:** We follow the proof of [8, Lemma 4.12]. For all $\psi_h \in \underline{V}_{h,\text{NC}}^k$,

$$ \begin{align*} c_h(\psi_h, \psi_h) & = c_{\mathcal{T}_h}(\psi_h, \psi_h) + c_{\mathcal{F}_h}(\psi_h, \psi_h) \\ & = c_{\mathcal{T}_h}(\psi_h, \psi_h) + \sum_{F \in \mathcal{F}_h} \frac{\alpha}{h_F} \|[[\psi_h]]\|_{\underline{L}(F)}^2 - 2 \left( \{\mathbb{D} \vec{\nabla}_h \psi_h\}, [[\psi_h]] \right)_{\mathcal{F}_h^i} \\ & \ge c_{\mathcal{T}_h}(\psi_h, \psi_h) + \alpha \|\psi_h\|_J^2 - 2C_{tr}(N_{\partial})^{1/2} \|\mathbb{D} \vec{\nabla}_h \psi_h\|_{\mathcal{T}_h} \|\psi_h\|_J \end{align*} $$
**Lemma 6** Suppose that $(\mathcal{T}_h)_h$ is a shape- and contact-regular mesh sequence. Then, we have for all $h > 0$:

$$ \forall \psi_h \in \underline{V}_{h,\text{NC}}^k, \forall K \in \mathcal{T}_h, \forall F \in \partial K, \quad h_K^{1/2} \|\psi_h\|_{\underline{L}^2(F)} \le C_{tr} \|\psi_h\|_{\underline{L}^2(K)}, \quad (17) $$

where $h_K$ is the diameter of element K.

We aim at asserting the discrete coercivity using the following norm:

$$ \forall \psi_h \in \underline{V}_{h,\text{NC}}^k, \quad \|\psi_h\|_{sip}^2 := c_{\mathcal{T}_h}(\psi_h, \psi_h) + \|\psi_h\|_J^2 $$

with the jump semi-norm

$$ \|\psi_h\|_J^2 := \sum_{F \in \mathcal{F}_h} \frac{1}{h_F} \|[[\psi_h]]\|_{\underline{L}(F)}^2. $$

Under assumption (4), there exists $\beta > 0$ we have for all $\psi_h \in \underline{V}_{h,\text{NC}}^k$

$$ c_{\mathcal{T}_h}(\psi_h, \psi_h) \ge \beta \left( \|\vec{\nabla}_h \psi_h\|_{\mathcal{T}_h}^2 + \|\psi_h\|_{\underline{L}}^2 \right), \quad (18) $$

so that

$$ \|\psi_h\|_{sip}^2 \ge \beta \left( \|\vec{\nabla}_h \psi_h\|_{\mathcal{T}_h}^2 + \|\psi_h\|_{\underline{L}}^2 + \|\psi_h\|_J^2 \right). $$

**Lemma 7 (Discrete coercivity)** Let $\underline{\alpha} := C_{tr}^2 N_{\partial} \frac{C_{\mathbb{D}}}{\beta}$ where

* $C_{tr}$ results from the discrete trace inequality (17),

* $N_{\partial}$ is defined in Section 4,

* $C_{\mathbb{D}}$ is defined in (11).

For all $\alpha \ge \underline{\alpha}$, the SIP bilinear form defined by (14) is coercive on $\underline{V}_{h,\text{NC}}^k$ with respect to the $\|\cdot\|_{sip}$-norm, i.e.,

$$ c_h(\psi_h, \psi_h) \ge C_{\alpha} \|\psi_h\|_{sip}^2, $$

$$ \text{with } C_{\alpha} := \left( \alpha - C_{tr}^2 N_{\partial} \frac{C_{\mathbb{D}}}{\beta} \right) \min \left\{ \frac{1}{2}, \beta \left( \alpha + C_{tr}^2 N_{\partial} \frac{C_{\mathbb{D}}}{\beta} \right)^{-1} \right\}. $$

**Proof:** We follow the proof of [8, Lemma 4.12]. For all $\psi_h \in \underline{V}_{h,\text{NC}}^k$,

$$ \begin{align*} c_h(\psi_h, \psi_h) & = c_{\mathcal{T}_h}(\psi_h, \psi_h) + c_{\mathcal{F}_h}(\psi_h, \psi_h) \\ & = c_{\mathcal{T}_h}(\psi_h, \psi_h) + \sum_{F \in \mathcal{F}_h} \frac{\alpha}{h_F} \|[[\psi_h]]\|_{\underline{L}(F)}^2 - 2 \left( \{\mathbb{D} \vec{\nabla}_h \psi_h\}, [[\psi_h]] \right)_{\mathcal{F}_h^i} \\ & \ge c_{\mathcal{T}_h}(\psi_h, \psi_h) + \alpha \|\psi_h\|_J^2 - 2C_{tr}(N_{\partial})^{1/2} \|\mathbb{D} \vec{\nabla}_h \psi_h\|_{\mathcal{T}_h} \|\psi_h\|_J \end{align*} $$
where we used Cauchy-Schwarz and Lemma 6 in the last line. Using the inequality $2ab \le \varepsilon a + \varepsilon^{-1}b$ for any $\varepsilon > 0$, we obtain

$$2C_{tr}(N_{\partial})^{1/2} \left\| \mathbb{D} \vec{\nabla}_h \psi_h \right\|_{\mathcal{T}_h} \|\psi_h\|_J \le \varepsilon C_{tr}^2 N_{\partial} \left\| \mathbb{D} \vec{\nabla}_h \psi_h \right\|_{\mathcal{T}_h}^2 + \varepsilon^{-1} \|\psi_h\|_J^2 \\ \le \varepsilon C_{tr}^2 N_{\partial} C_{\mathbb{D}} \left\| \vec{\nabla}_h \psi_h \right\|_{\mathcal{T}_h}^2 + \varepsilon^{-1} \|\psi_h\|_J^2.$$

Using (18), we obtain that there exists a constant $\beta > 0$ such that

$$c_h(\psi_h, \psi_h) \\ \ge \beta(1 - \varepsilon\alpha) \left\| \vec{\nabla}_h \psi_h \right\|_{\mathcal{T}_h}^2 + \beta \|\psi_h\|_{\underline{L}}^2 + (\alpha - \varepsilon^{-1}) \|\psi_h\|_J^2.$$

Choosing $\varepsilon = 2(\alpha + \alpha)^{-1}$ yields the assertion. $\square$

Thus, it only remains to prove boundedness. To this purpose, we need to define $\underline{V}^{*,h} = \underline{V}^* + \underline{V}_{h,NC}^k$ and the following norm

$$\|\psi\|_{sip,\star} := \left( \|\psi\|_{sip}^p + \sum_{K \in \mathcal{T}_h} h_K^{1+\gamma_p} \|\vec{\nabla} \psi|_K \cdot \mathbf{n}_K \|_{\underline{L}^p(\partial K)} \right)^{1/p},$$

where $\gamma_p = \frac{d(p-2)}{2}$ and $\mathbf{n}_K$ is the unit outward normal to $K$. Following [8, Section 4.2], we obtain the following results.

**Lemma 8 (Boundedness)** There is $C_{bnd}$, independent of $h$, such that for all $(\phi, \psi_h) \in \underline{V}^{*,h} \times \underline{V}_h$

$$c_h(\phi, \psi_h) \le C_{bnd} \|\phi\|_{sip,\star} \|\psi_h\|_{sip}$$

**Theorem 9 (Convergence)** Suppose that there exists $r_{max}$ in $(0, 1]$ such that $\forall r \in [0, r_{max}]$, $\phi \in (H^{1+r}(\mathbb{R}))^{\overline{N} \times G}$ ([6], Proposition 1). Then the solution of (15), $\phi_h$ is such that:

$$\|\phi - \phi_h\|_{sip} \lesssim C \inf_{\psi_h \in \underline{V}_{h,NC}} \|\phi - \psi_h\|_{sip,\star},$$

where $C$ is a constant independent of $h$. Moreover, under Assumption 5.1, there holds

$$\|\phi - \phi_h\|_{sip} \le C |\phi|_{W^{2,p}(\mathcal{T}_h)} h^{\mu},$$

where $\mu = r_{max}$, $C$ is a constant independent of $h$ and $p$ is such that $\mu = \frac{d+2}{2} - \frac{d}{p}$.

**Theorem 10 (L²-norm estimate)** Suppose that there exists $r_{max}$ in $(0, 1]$ such that $\forall r \in [0, r_{max}]$, $\phi_m^g \in H^{1+r}(\mathbb{R})$ ([6], Proposition 1). Under Assumption 5.1, the solution of (15), $\phi_h$ is such that: $\|\phi - \phi_h\|_{\underline{L}} \lesssim h^{2\mu} \|S_f\|_{\underline{L}}$, where $\mu = r_{max}$.

**Proof:** We apply the Aubin-Nitsche similarly as in [8, Theorem 4.25]. $\square$
where we used Cauchy-Schwarz and Lemma 6 in the last line. Using the inequality $2ab \le \varepsilon a + \varepsilon^{-1}b$ for any $\varepsilon > 0$, we obtain

$$2C_{tr}(N_{\partial})^{1/2} \left\| \mathbb{D} \vec{\nabla}_h \psi_h \right\|_{\mathcal{T}_h} \|\psi_h\|_J \le \varepsilon C_{tr}^2 N_{\partial} \left\| \mathbb{D} \vec{\nabla}_h \psi_h \right\|_{\mathcal{T}_h}^2 + \varepsilon^{-1} \|\psi_h\|_J^2 \\ \le \varepsilon C_{tr}^2 N_{\partial} C_{\mathbb{D}} \left\| \vec{\nabla}_h \psi_h \right\|_{\mathcal{T}_h}^2 + \varepsilon^{-1} \|\psi_h\|_J^2.$$

Using (18), we obtain that there exists a constant $\beta > 0$ such that

$$c_h(\psi_h, \psi_h) \\ \ge \beta(1 - \varepsilon\alpha) \left\| \vec{\nabla}_h \psi_h \right\|_{\mathcal{T}_h}^2 + \beta \|\psi_h\|_{\underline{L}}^2 + (\alpha - \varepsilon^{-1}) \|\psi_h\|_J^2.$$

Choosing $\varepsilon = 2(\alpha + \alpha)^{-1}$ yields the assertion. $\square$

Thus, it only remains to prove boundedness. To this purpose, we need to define $\underline{V}^{*,h} = \underline{V}^* + \underline{V}_{h,NC}^k$ and the following norm

$$\|\psi\|_{sip,\star} := \left( \|\psi\|_{sip}^p + \sum_{K \in \mathcal{T}_h} h_K^{1+\gamma_p} \|\vec{\nabla} \psi|_K \cdot \mathbf{n}_K \|_{\underline{L}^p(\partial K)} \right)^{1/p},$$

where $\gamma_p = \frac{d(p-2)}{2}$ and $\mathbf{n}_K$ is the unit outward normal to $K$. Following [8, Section 4.2], we obtain the following results.

**Lemma 8 (Boundedness)** There is $C_{bnd}$, independent of $h$, such that for all $(\phi, \psi_h) \in \underline{V}^{*,h} \times \underline{V}_h$

$$c_h(\phi, \psi_h) \le C_{bnd} \|\phi\|_{sip,\star} \|\psi_h\|_{sip}$$

**Theorem 9 (Convergence)** Suppose that there exists $r_{max}$ in $(0, 1]$ such that $\forall r \in [0, r_{max}]$, $\phi \in (H^{1+r}(\mathbb{R}))^{\bar{N} \times G}$ ([6], Proposition 1). Then the solution of (15), $\phi_h$ is such that:

$$\|\phi - \phi_h\|_{sip} \lesssim C \inf_{\psi_h \in \underline{V}_{h,NC}} \|\phi - \psi_h\|_{sip,\star},$$

where $C$ is a constant independent of $h$. Moreover, under Assumption 5.1, there holds

$$\|\phi - \phi_h\|_{sip} \le C |\phi|_{W^{2,p}(\mathcal{T}_h)} h^{\mu},$$

where $\mu = r_{max}$, $C$ is a constant independent of $h$ and $p$ is such that $\mu = \frac{d+2}{2} - \frac{d}{p}$.

**Theorem 10 (L²-norm estimate)** Suppose that there exists $r_{max}$ in $(0, 1]$ such that $\forall r \in [0, r_{max}]$, $\phi_m^g \in H^{1+r}(\mathbb{R})$ ([6], Proposition 1). Under Assumption 5.1, the solution of (15), $\phi_h$ is such that: $\|\phi - \phi_h\|_{\underline{L}} \lesssim h^{2\mu} \|S_f\|_{\underline{L}}$, where $\mu = r_{max}$.

**Proof:** We apply the Aubin-Nitsche similarly as in [8, Theorem 4.25]. $\square$
# 6 The eigenproblem

## 6.1 Conforming discretization

**Theorem 11** Let $\mu$ be the regularity of the eigenfunction $\varphi$ associated to $\lambda$, and $\omega = \min(\mu, k)$. Let $\lambda_h$ be the discrete eigenvalue associated to Problem (13). The following a priori error estimate holds: $|\lambda - \lambda_h| \lesssim h^{2\omega}$.

**Proof:** As in the continuous case (Theorem 2), since the discretization is conforming, there exists a unique compact operator $T_h: \underline{V}_h^k \to \underline{V}_h^k$ such that $\forall (\phi_h, \psi_h) \in \underline{V}_h^k \times \underline{V}_h^k: c(T_h\phi_h, \psi_h) = \ell_f(\phi_h, \psi_h)$. According to Thm. 5, the sequence of the operators $(T_h)_h$ is pointwise converging towards $T$. As $T_h$ and $T$ are a compact operators, the sequence of the operators $(T_h)_h$ is then converging in $\mathcal{L}(\underline{V})$ towards $T: \|T_h - T\|_{\mathcal{L}(\underline{V})} \to 0$. The norm convergence guarantees that there is no spectral pollution (see [17]). Moreover, we can apply Theorem 8.3 in [2] to state the error estimate on the eigenvalue. We remark that $(M_f\phi, \phi)_L$ is a norm over $\underline{V}_\lambda := \{\phi \in \underline{V} | \forall \psi \in \underline{V}, c(\phi, \psi) = \lambda \ell_f(\phi, \psi)\}$ [12, Section 5.2.2 p. 78]. $\square$

## 6.2 SIPG discretization

We recall that, in this section, we work under the assumption 5.1.

**Theorem 12** Let $\mu$ be the regularity of the eigenfunction $\varphi$ associated to $\lambda$, and $\omega = \min(\mu, k)$. Let $\lambda_h$ be the discrete eigenvalue associated to Problem (16). The following a priori error estimate holds: $|\lambda - \lambda_h| \lesssim h^{2\omega}$.

**Proof:** We apply the theory developed in [1]. The proof is decomposed as follows. We first show that there is no spectral pollution. Then, we derive the error estimate.

Let $E: \underline{V} + \underline{V}_{h,NC}^k \to \underline{V} + \underline{V}_{h,NC}^k$ be the continuous spectral projector relative to $\lambda$ defined by

$$E = \frac{1}{2\pi i} \int_{\Gamma} \left( z - T|_{\underline{V} + \underline{V}_{h,NC}^k} \right)^{-1} dz,$$

where $\Gamma$ is a circle in the complex plane centred at $\lambda$ which lies in $\rho(T|_{\underline{V} + \underline{V}_{h,NC}^k})$ and encloses no other points of $\sigma(T|_{\underline{V} + \underline{V}_{h,NC}^k})$. The absence of spectral pollution relies on two properties. First, using interpolation results [8, Assumption 4.31] we have for all $\phi \in E(\underline{V} + \underline{V}_{h,NC}^k)$,

$$\inf_{\psi_h \in \underline{V}_{h,NC}^k} \|\phi - \psi_h\|_{sip} \le Ch^\mu,$$
# 6 The eigenproblem

## 6.1 Conforming discretization

**Theorem 11** Let $\mu$ be the regularity of the eigenfunction $\varphi$ associated to $\lambda$, and $\omega = \min(\mu, k)$. Let $\lambda_h$ be the discrete eigenvalue associated to Problem (13). The following a priori error estimate holds: $|\lambda - \lambda_h| \lesssim h^{2\omega}$.

**Proof:** As in the continuous case (Theorem 2), since the discretization is conforming, there exists a unique compact operator $T_h: \underline{V}_h^k \to \underline{V}_h^k$ such that $\forall (\phi_h, \psi_h) \in \underline{V}_h^k \times \underline{V}_h^k: c(T_h\phi_h, \psi_h) = \ell_f(\phi_h, \psi_h)$. According to Thm. 5, the sequence of the operators $(T_h)_h$ is pointwise converging towards $T$. As $T_h$ and $T$ are a compact operators, the sequence of the operators $(T_h)_h$ is then converging in $\mathcal{L}(\underline{V})$ towards $T: \|T_h - T\|_{\mathcal{L}(\underline{V})} \to 0$. The norm convergence guarantees that there is no spectral pollution (see [17]). Moreover, we can apply Theorem 8.3 in [2] to state the error estimate on the eigenvalue. We remark that $(M_f\phi, \phi)_L$ is a norm over $\underline{V}_\lambda := \{\phi \in \underline{V} | \forall \psi \in \underline{V}, c(\phi, \psi) = \lambda \ell_f(\phi, \psi)\}$ [12, Section 5.2.2 p. 78]. $\square$

## 6.2 SIPG discretization

We recall that, in this section, we work under the assumption 5.1.

**Theorem 12** Let $\mu$ be the regularity of the eigenfunction $\varphi$ associated to $\lambda$, and $\omega = \min(\mu, k)$. Let $\lambda_h$ be the discrete eigenvalue associated to Problem (16). The following a priori error estimate holds: $|\lambda - \lambda_h| \lesssim h^{2\omega}$.

**Proof:** We apply the theory developed in [1]. The proof is decomposed as follows. We first show that there is no spectral pollution. Then, we derive the error estimate.

Let $E: \underline{V} + \underline{V}_{h,NC}^k \to \underline{V} + \underline{V}_{h,NC}^k$ be the continuous spectral projector relative to $\lambda$ defined by

$$E = \frac{1}{2\pi i} \int_{\Gamma} \left( z - T|_{\underline{V} + \underline{V}_{h,NC}^k} \right)^{-1} dz,$$

where $\Gamma$ is a circle in the complex plane centred at $\lambda$ which lies in $\rho(T|_{\underline{V} + \underline{V}_{h,NC}^k})$ and encloses no other points of $\sigma(T|_{\underline{V} + \underline{V}_{h,NC}^k})$. The absence of spectral pollution relies on two properties. First, using interpolation results [8, Assumption 4.31] we have for all $\phi \in E(\underline{V} + \underline{V}_{h,NC}^k)$,

$$\inf_{\psi_h \in \underline{V}_{h,NC}^k} \|\phi - \psi_h\|_{sip} \le Ch^\mu,$$
where $C$ is a constant independent of $h$. Second, we have for all $\phi_h \in \underline{V}_{h,NC}^k$,

$$ \begin{aligned} \|(T - T_h)\phi_h\|_{sip} &\le Ch^\mu |T\phi_h|_{\underline{W}^{2,p}(T_h)} \\ &\le Ch^\mu \|T\phi_h\|_{(H^{1+\alpha_p}(\mathbb{R}))^{\tilde{N}\times G}} \\ &\le Ch^\mu \|\phi_h\|_{\underline{L}} \\ &\le Ch^\mu \|\phi_h\|_{sip}, \end{aligned} $$

where we used Theorem 9 in the second line and regularity results [15] in the third line. Applying [1, Theorem 3.7], we obtain that there is no spectral pol- lution.

Moreover, we apply [1, Theorem 3.14] to state the error estimate on the eigenvalue,

$$ |\lambda - \lambda_h| \le C \delta_h \delta_{*,h}, $$

where

$$ \begin{aligned} \delta_h &= \gamma_h + \left\| (T - T_h) |_{E(\underline{V} + \underline{V}_{h,NC}^k)} \right\|_{sip} \\ \delta_{*,h} &= \gamma_{*,h} + \left\| (T_* - T_{*,h}) |_{E(\underline{V} + \underline{V}_{h,NC}^k)} \right\|_{sip}, \end{aligned} $$

with

$$ \begin{aligned} \gamma_h &= \delta(E(V + \underline{V}_{h,NC}^k), \underline{V}_{h,NC}^k), \\ \gamma_{*,h} &= \delta(E_*(V + \underline{V}_{h,NC}^k), \underline{V}_{h,NC}^k), \end{aligned} $$

where

$$ \delta(Y, Z) = \sup_{y \in Y, \|y\|_{sip} = 1} \left( \inf_{z \in Z} \|y - z\|_{sip} \right), $$

and $E_* : \underline{V} + \underline{V}_{h,NC}^k \to \underline{V} + \underline{V}_{h,NC}^k$ is the continuous spectral projector of the adjoint operator $T_*|_{\underline{V}+\underline{V}_{h,NC}^k}$ relative to $\bar{\lambda}$.

Using again elliptic regularity results [15] and Theorem 9, we obtain

$$ \begin{aligned} \left\| (T - T_h) |_{E(\underline{V} + \underline{V}_{h,NC}^k)} \right\|_{sip} &\le Ch^\mu \\ \left\| (T_* - T_{*,h}) |_{E(\underline{V} + \underline{V}_{h,NC}^k)} \right\|_{sip} &\le Ch^\mu. \end{aligned} $$

Using elliptic regularity results, we get

$$ \|\varphi\|_{(H^{1+\alpha_p}(\mathbb{R}))^{\tilde{N}\times G}} \le C \|\varphi\|_{\underline{L}} \le C \|\varphi\|_{\underline{V}}. $$
where C is a constant independent of h. Second, we have for all $\phi_h \in \underline{V}_{h,NC}^k$,

$$ \begin{aligned} \|(T - T_h)\phi_h\|_{sip} &\le Ch^\mu |T\phi_h|_{\underline{W}^{2,p}(T_h)} \\ &\le Ch^\mu \|T\phi_h\|_{(H^{1+\alpha_p}(\mathbb{R}))^{\tilde{N}\times G}} \\ &\le Ch^\mu \|\phi_h\|_{\underline{L}} \\ &\le Ch^\mu \|\phi_h\|_{sip}, \end{aligned} $$

where we used Theorem 9 in the second line and regularity results [15] in the third line. Applying [1, Theorem 3.7], we obtain that there is no spectral pol- lution.

Moreover, we apply [1, Theorem 3.14] to state the error estimate on the eigenvalue,

$$ |\lambda - \lambda_h| \le C \delta_h \delta_{*,h}, $$

where

$$ \begin{aligned} \delta_h &= \gamma_h + \left\| (T - T_h) |_{E(\underline{V} + \underline{V}_{h,NC}^k)} \right\|_{sip} \\ \delta_{*,h} &= \gamma_{*,h} + \left\| (T_* - T_{*,h}) |_{E(\underline{V} + \underline{V}_{h,NC}^k)} \right\|_{sip}, \end{aligned} $$

with

$$ \begin{aligned} \gamma_h &= \delta(E(V + \underline{V}_{h,NC}^k), \underline{V}_{h,NC}^k), \\ \gamma_{*,h} &= \delta(E_*(V + \underline{V}_{h,NC}^k), \underline{V}_{h,NC}^k), \end{aligned} $$

where

$$ \delta(Y, Z) = \sup_{y \in Y, \|y\|_{sip} = 1} \left( \inf_{z \in Z} \|y - z\|_{sip} \right), $$

and $E_* : \underline{V} + \underline{V}_{h,NC}^k \to \underline{V} + \underline{V}_{h,NC}^k$ is the continuous spectral projector of the adjoint operator $T_*|_{\underline{V}+\underline{V}_{h,NC}^k}$ relative to $\bar{\lambda}$.

Using again elliptic regularity results [15] and Theorem 9, we obtain

$$ \begin{aligned} \left\| (T - T_h) |_{E(\underline{V} + \underline{V}_{h,NC}^k)} \right\|_{sip} &\le Ch^\mu \\ \left\| (T_* - T_{*,h}) |_{E(\underline{V} + \underline{V}_{h,NC}^k)} \right\|_{sip} &\le Ch^\mu. \end{aligned} $$

Using elliptic regularity results, we get

$$ \|\varphi\|_{(H^{1+\alpha_p}(\mathbb{R}))^{\tilde{N}\times G}} \le C \|\varphi\|_{\underline{L}} \le C \|\varphi\|_{\underline{V}}. $$
Applying Theorem 9, we infer that

$$ \gamma_h \le Ch^\mu \\ \gamma_{*,h} \le Ch^\mu. $$

This concludes the proof.

<figure><img src="image_2.png" /></figure>

# 7 Numerical Results

We consider the test case Model 2, case 1 from the benchmark of Takeda and Ikeda [19]. The geometry of the core is three-dimensional and the domain is $\{(x, y, z) \in \mathbb{R}^3, 0 \le x \le 140 \text{ cm}; 0 \le y \le 140 \text{ cm}; 0 \le z \le 150 \text{ cm}\}$. This test is defined with 4 energy groups, isotropic scattering and vacuum boundary conditions. Figure 1 represents the cross-sectional geometry on the plane $z = 75 \text{ cm}$.

Since the scattering is isotropic, the $SP_3$ formulation can easily be reformulated as a multigroup diffusion problem with 8 energy groups and an isotropic albedo boundary condition [3]. We then made the computations with the PRIAM solver from the code CRONOS2 [14] for the conforming case and with the MINARET solver [13] from the APOLLO3® code [18] for the SIPG discretization.

<figure><img src="image_3.png" /><figcaption>Figure 1: Cross-sectional view of the core ($z = 75 \text{ cm}$).</figcaption></figure>

In Figure 2, we consider the convergence of the fundamental mode where we used the $SP_3$ formulation with $Q^1$ finite elements and a regular cartesian mesh of size $h$. The approximated order of convergence is 2.22.

In Figure 3, we consider the convergence of the fundamental mode for different the $SP_N$ formulations with discontinuous $P^1$ finite elements and a prismatic mesh of size $h$. The approximated orders of convergence are given in Table 1.
Applying Theorem 9, we infer that

$$ \gamma_h \le Ch^\mu \\ \gamma_{*,h} \le Ch^\mu. $$

This concludes the proof.

<figure><img src="image_3.png" /></figure>

# 7 Numerical Results

We consider the test case Model 2, case 1 from the benchmark of Takeda and Ikeda [19]. The geometry of the core is three-dimensional and the domain is $\{(x, y, z) \in \mathbb{R}^3, 0 \le x \le 140 \text{ cm}; 0 \le y \le 140 \text{ cm}; 0 \le z \le 150 \text{ cm}\}$. This test is defined with 4 energy groups, isotropic scattering and vacuum boundary conditions. Figure 1 represents the cross-sectional geometry on the plane $z = 75 \text{ cm}$.

Since the scattering is isotropic, the $SP_3$ formulation can easily be reformulated as a multigroup diffusion problem with 8 energy groups and an isotropic albedo boundary condition [3]. We then made the computations with the PRIAM solver from the code CRONOS2 [14] for the conforming case and with the MINARET solver [13] from the APOLLO3® code [18] for the SIPG discretization.

<figure><img src="image_4.png" /><figcaption>Figure 1: Cross-sectional view of the core ($z = 75 \text{ cm}$).</figcaption></figure>

In Figure 2, we consider the convergence of the fundamental mode where we used the $SP_3$ formulation with $Q^1$ finite elements and a regular cartesian mesh of size $h$. The approximated order of convergence is 2.22.

In Figure 3, we consider the convergence of the fundamental mode for different the $SP_N$ formulations with discontinuous $P^1$ finite elements and a prismatic mesh of size $h$. The approximated orders of convergence are given in Table 1.
<figure><img src="image_4.png" /><figcaption>Figure 2: Error on the discrete eigenvalue for the $SP_3$ formulation with $Q^1$ finite elements</figcaption></figure>

<table><caption>Table 1: Approximated order of convergence associated to Figure 3</caption><tr><td>$SP_3$</td><td>$SP_5$</td><td>$SP_7$</td></tr><tr><td>1.88</td><td>1.96</td><td>1.92</td></tr></table>

# 8 Conclusion

We did the numerical analysis of the approximation with an $H^1$-conforming finite element method of the neutron multigroup $SP_N$ equations. We also studied the numerical analysis of the approximation with the Symmetric Interior Penalty Galerkin method of the neutron multigroup $SP_N$ equations. We then illustrated numerically the convergence results on a benchmark representative of a nuclear core. Those results can be extended to a mixed finite element method, see [5] for the diffusion case with an $H^1$-conforming finite element method.

# Acknowledgements

The authors gratefully acknowledge P. Ciarlet for fruitful discussions.

# References

[1] A. Alonso and A. D. Russo. Spectral approximation of variationally-posed eigenvalue problems by nonconforming methods. *Journal of Computational*
<figure><img src="image_5.png" /><figcaption>Figure 2: Error on the discrete eigenvalue for the $SP_3$ formulation with $Q^1$ finite elements</figcaption></figure>

<table><caption>Table 1: Approximated order of convergence associated to Figure 3</caption><tr><td>$SP_3$</td><td>$SP_5$</td><td>$SP_7$</td></tr><tr><td>1.88</td><td>1.96</td><td>1.92</td></tr></table>

# 8 Conclusion

We did the numerical analysis of the approximation with an $H^1$-conforming finite element method of the neutron multigroup $SP_N$ equations. We also studied the numerical analysis of the approximation with the Symmetric Interior Penalty Galerkin method of the neutron multigroup $SP_N$ equations. We then illustrated numerically the convergence results on a benchmark representative of a nuclear core. Those results can be extended to a mixed finite element method, see [5] for the diffusion case with an $H^1$-conforming finite element method.

# Acknowledgements

The authors gratefully acknowledge P. Ciarlet for fruitful discussions.

# References

[1] A. Alonso and A. D. Russo. Spectral approximation of variationally-posed eigenvalue problems by nonconforming methods. *Journal of Computational*
<figure><img src="image_5.png" /><figcaption>Figure 3: Error on the discrete eigenvalue for the $SP_3$ formulation with discontinuous linear finite elements</figcaption></figure>

and *Applied Mathematics*, 223(1):177-197, 2009.

[2] I. Babuska and J. Osborn. Eigenvalue problems, in handbook of numerical analysis, vol. II. *PG Ciarlet and JL Lions Editions, Elsevier science Publishers BV (North-Holland)*, pages 645-785, 1991.

[3] A.-M. Baudron and J.-J. Lautard. Simplified $P_N$ transport core calculations in the Apollo3 system. *International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering (M&C 2011)*, Rio de Janeiro, RJ, Brazil, May 8-12, 2011.

[4] H. Brézis. Analyse fonctionnelle, théorie et applications, 1983.

[5] P. Ciarlet, L. Giret, E. Jamelot, and F.D. Kpadonou. Numerical analysis of the mixed finite element method for the neutron diffusion eigenproblem with heterogeneous coefficients. *ESAIM: Mathematical Modelling and Numerical Analysis*, 52(5):2003-2035, 2018.

[6] P. Ciarlet Jr, E. Jamelot, and F.D. Kpadonou. Domain decomposition methods for the diffusion equation with low-regularity solution. *Computers & Mathematics with Applications*, 74(10):2369-2384, 2017.

[7] R. Dautray and J.-L. Lions. *Analyse mathématique et calcul numérique pour les sciences et les techniques*. Masson, 1985.

[8] D. A. Di Pietro and A. Ern. *Mathematical aspects of discontinuous Galerkin methods*, volume 69. Springer Science & Business Media, 2011.
<figure><img src="image_6.png" /><figcaption>Figure 3: Error on the discrete eigenvalue for the $SP_3$ formulation with discontinuous linear finite elements</figcaption></figure>

and *Applied Mathematics*, 223(1):177-197, 2009.

[2] I. Babuska and J. Osborn. Eigenvalue problems, in handbook of numerical analysis, vol. II. *PG Ciarlet and JL Lions Editions, Elsevier science Publishers BV (North-Holland)*, pages 645-785, 1991.

[3] A.-M. Baudron and J.-J. Lautard. Simplified $P_N$ transport core calculations in the Apollo3 system. *International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering (M&C 2011)*, Rio de Janeiro, RJ, Brazil, May 8-12, 2011.

[4] H. Brézis. Analyse fonctionnelle, théorie et applications, 1983.

[5] P. Ciarlet, L. Giret, E. Jamelot, and F.D. Kpadonou. Numerical analysis of the mixed finite element method for the neutron diffusion eigenproblem with heterogeneous coefficients. *ESAIM: Mathematical Modelling and Numerical Analysis*, 52(5):2003-2035, 2018.

[6] P. Ciarlet Jr, E. Jamelot, and F.D. Kpadonou. Domain decomposition methods for the diffusion equation with low-regularity solution. *Computers & Mathematics with Applications*, 74(10):2369-2384, 2017.

[7] R. Dautray and J.-L. Lions. *Analyse mathématique et calcul numérique pour les sciences et les techniques*. Masson, 1985.

[8] D. A. Di Pietro and A. Ern. *Mathematical aspects of discontinuous Galerkin methods*, volume 69. Springer Science & Business Media, 2011.
[9] J. J. Duderstadt and L. J. Hamilton. *Nuclear reactor analysis*, volume 650. John Wiley & Sons, Inc., 1976.

[10] A. Ern and J.-L. Guermond. *Theory and practice of finite elements*, volume 159. Springer Science & Business Media, 2013.

[11] E. M. Gelbard. Application of spherical harmonics method to reactor problems. *Bettis Atomic Power Laboratory, West Mifflin, PA, Technical Report No. WAPD-BT-20*, 1960.

[12] L. Giret. *Non-Conforming Domain Decomposition for the Multigroup Neutron SPN Equation*. PhD thesis, Paris Saclay, 2018.

[13] J.-J. Lautard and J.-Y. Moller. Minaret, a deterministic neutron transport solver for nuclear core calculations. In M&C, Rio de Janeiro, 2011.

[14] J.J. Lautard, S. Loubière, and C. Fedon-Magnaud. Cronos: a modular computational system for neutronic core calculations. 1992.

[15] S. Nicaise and A.-M. Sändig. General interface problems. I, II. *Mathematical Methods in the Applied Sciences*, 17(6):395–429, 431–450, 1994.

[16] ie: for $1 \le g \le G$, $m \in \mathcal{I}_e$, $(\phi_m^g)|_{\partial\mathcal{R}} = 0$.

[17] J. E. Osborn. Spectral approximation for compact operators. *Mathematics of computation*, 29(131):712–725, 1975.

[18] D Schneider et al. APOLLO3®: CEA/DEN deterministic multi-purpose code for reactor physics analysis. PHYSOR-2016, May 1-5 2016, Sun Valley, Idaho, USA.

[19] T. Takeda and H. Ikeda. 3-d neutron transport benchmarks. *Journal of Nuclear Science and Technology*, 28(7):656–669, 1991.
[9] J. J. Duderstadt and L. J. Hamilton. *Nuclear reactor analysis*, volume 650. John Wiley & Sons, Inc., 1976.

[10] A. Ern and J.-L. Guermond. *Theory and practice of finite elements*, volume 159. Springer Science & Business Media, 2013.

[11] E. M. Gelbard. Application of spherical harmonics method to reactor problems. *Bettis Atomic Power Laboratory, West Mifflin, PA, Technical Report No. WAPD-BT-20*, 1960.

[12] L. Giret. *Non-Conforming Domain Decomposition for the Multigroup Neutron SPN Equation*. PhD thesis, Paris Saclay, 2018.

[13] J.-J. Lautard and J.-Y. Moller. Minaret, a deterministic neutron transport solver for nuclear core calculations. In M&C, Rio de Janeiro, 2011.

[14] J.J. Lautard, S. Loubière, and C. Fedon-Magnaud. Cronos: a modular computational system for neutronic core calculations. 1992.

[15] S. Nicaise and A.-M. Sändig. General interface problems. I, II. *Mathematical Methods in the Applied Sciences*, 17(6):395–429, 431–450, 1994.

[16] ie: for $1 \le g \le G$, $m \in \mathcal{I}_e$, $(\phi_m^g)|_{\partial\mathcal{R}} = 0$.

[17] J. E. Osborn. Spectral approximation for compact operators. *Mathematics of computation*, 29(131):712–725, 1975.

[18] D Schneider et al. APOLLO3®: CEA/DEN deterministic multi-purpose code for reactor physics analysis. PHYSOR-2016, May 1-5 2016, Sun Valley, Idaho, USA.

[19] T. Takeda and H. Ikeda. 3-d neutron transport benchmarks. *Journal of Nuclear Science and Technology*, 28(7):656–669, 1991.