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

HAL
open science

# A new AI based computation method of the Eddington factor in the M1-multigroup model

Gonzague Radureau, Claire Michaut

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

## To cite this version:

Gonzague Radureau, Claire Michaut. A new AI based computation method of the Eddington factor in the M1-multigroup model. 2024. hal-04661994

HAL Id: hal-04661994

https://hal.science/hal-04661994v1

Preprint submitted on 25 Jul 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.

Public Domain
# A new AI based computation method of the Eddington factor in the M1-multigroup model

Gonzague Radureauª, Claire Michautª

ªUniversité Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Nice, 06300, France

## Abstract

The M1-multigroup model intricately captures the interaction between light and matter, incorporating the spectral behaviour of photons. However, a critical aspect of this model lies in determining the Eddington factor, which is used in the closure relation linking radiative pressure to radiative energy. Despite lacking an analytical expression, our investigation unveils that the Eddington factor depends solely on three parameters: the *radiative temperature*, the *reduced flux*, and the *frequency bounds' ratio* of a considered group. To address this challenge, we have devised a novel approach leveraging neural networks and polynomials, enabling rapid and accurate estimation of this quantity.

Our method showcases significant advantages over existing techniques. It demonstrates speeds up to 3,000 times faster than the most precise method utilising a line search algorithm and achieves precision levels up to 1,000 times higher than those relying on the M1-grey model's Eddington factor expression. Moreover, unlike interpolation-based methods, our approach eliminates the need for prior knowledge of radiative quantities. Consequently, our method emerges as one of the most efficient means to precisely compute the Eddington factor in the M1-multigroup model, offering a potent tool for advancing radiative hydrodynamics simulations.

**Keywords:** Radiation Hydrodynamics, Radiative Transfer, Artificial Intelligence, Closure Relation, Eddington Factor

## 1. Introduction

Radiative hydrodynamics models the coupling between a hypersonic hot plasma's dynamics and self-produced or external radiation. Almost all numerical codes use simplified models, such as the cooling functions or the diffusion approximation, that are, in most cases, limited in the general study of the interaction of photons with matter. To accurately model the photon transport, the HADES 2D code was specifically developed [1, 2, 3]. Such a code is indispensable for studying astrophysical objects, in which optically intermediate regions are still poorly modelled yet commonly encountered within such phenomena.

This code couples the hydrodynamics with the M1-grey and M1-multigroup models [4, 5, 6] for radiative transfer. The M1-multigroup model considers the spectral behaviour of light by partitioning the electromagnetic spectrum into groups. These models are based on the first two moments of the radiative transfer equations and require a closure relation between the second moment (radiative pressure) and the first (radiative energy) via the so-called Eddington tensor. This tensor depends on a scalar quantity called the Eddington factor [7, 8, 9], which has an analytical expression in the M1-grey model [10] but not in the M1-multigroup model. To compute the Eddington factor in the multigroup model, there exist three methods in the literature: the search algorithms [1], which are computationally expensive but precise; the interpolation of precalculated values [1, 11], which requires prior knowledge of the radiative quantities; and the expression of the M1-grey model [12],

which the fastest method, but is not valid, as we will see in this paper.

Historically, the HADES code implemented two methods to compute the Eddington factor. The first method is based on search algorithms such as line search and Dichotomy-Newton's method [1], offers high precision but come with significant computational costs, thus limiting access to simulations using the M1-multigroup model. Additionally, this method occasionally fails to converge, resulting in instabilities within the simulations. The second method utilises interpolations [1, 11], which require prior knowledge of radiative quantities, a requirement that is not always met in radiative hydrodynamics simulations.

Neural networks have previously been employed in supernova numerical simulations to compute the Eddington tensor of the M1 model, which is used to model the interactions between neutrinos and matter [13]. Notably, the closure relation needs to have the required accuracy for this purpose. This study has shown the efficiency of neural networks in achieving reasonably precise results while keeping computational costs manageable. This finding underscores the potential of neural networks to significantly improve the accuracy and efficiency of supernova simulations compared to traditional methods.

In this paper, we introduce a new AI-based method developed within HADES, which computes precise values of the Eddington factor while effectively reducing computational expenses. First, we will analyse the dependencies of the Eddington factor in the M1-multigroup model (Section 2). Next, we will explore the intricacies of our devised approach (Section 3). Subsequently, we will conduct a comparative analysis, evaluating the performance of our method against other existing meth-

Email address: gonzague.radureau@oca.eu (Gonzague Radureau)
ods (Section 4). Finally, we will present our conclusions (Section 5). As the mathematical developments are complex, we detail the calculations in the seven appendices.

## 2. Analysis of the Eddington factor in the M1-multigroup model

We define all the variables used in Table 1.

### 2.1. Definition of the radiative quantities

The radiation is defined by its specific intensity $I_v$, which depends on the frequency of the photons and their propagation direction. It reflects the quantity of energy transported by radiation per surface unit. One can derive respectively the radiative energy, the radiative flux and the radiative pressure as the two first moments of the specific intensity:

$$E_v = \frac{1}{c} \int_{4\pi} I_v(\Omega) \, d\Omega \; ; \qquad (1)$$

$$\mathbf{F}_v = \int_{4\pi} \Omega I_v(\Omega) \, d\Omega \; ; \qquad (2)$$

$$\mathbb{P}_v = \frac{1}{c} \int_{4\pi} \Omega \otimes \Omega I_v(\Omega) \, d\Omega \; . \qquad (3)$$

Where c is the light velocity, v is the photon frequency, Ω is a given direction, and ⊗ is the dyadic product. The radiative pressure is related to the radiative energy via the Eddington tensor D_g according to the following expression:

$$\mathbb{P}_v = \mathbb{D}_v E_v \; . \qquad (4)$$

According to the M1 model, which assumes that the radiation propagates symmetrically according to the direction of the radiative flux, the Eddington tensor is further related to a scalar quantity called the Eddington factor $\chi_v$. They are related according to the following relation:

$$\mathbb{D}_v = \frac{1 - \chi_v}{2} \mathbb{I} + \frac{3\chi_v - 1}{2} \frac{\mathbf{F}_v \otimes \mathbf{F}_v}{\|\mathbf{F}_v\|^2} \; . \qquad (5)$$

Of course, it is impossible to consider the entirety of the electromagnetic spectrum, which is continuous. Thus there exist two models: one that averages the behaviour of photons in the electromagnetic spectrum (M1-grey model) and another that discretises the electromagnetic spectrum to achieve more precision (M1-multigroup model).

#### 2.1.1. M1-grey model

The M1-grey model considers all the photons as a whole. Consequently, radiative quantities are integrated over all frequencies according to the following expressions:

$$\left\{ \begin{aligned} \mathbb{E}_R &= \int_0^\infty \mathbb{E}_v \, d\nu \\ \mathbf{F}_R &= \int_0^\infty \mathbf{F}_v \, d\nu \; . \\ \mathbb{P}_R &= \int_0^\infty \mathbb{P}_v \, d\nu \end{aligned} \right. \qquad (6)$$

In this case, the M1 closure relation (4) holds, and the Eddington tensor can still be expressed as a function of the Eddington factor following Eq. (5). In this case, the Eddington

factor admits an analytical expression [10], depending solely on the reduced flux, which is defined by the relation:

$$f_R = \frac{\|\mathbf{F}_R\|}{c \mathbb{E}_R} \qquad (7)$$

Its values lie in [0, 1] and indicate the degree of anisotropy of the radiation ($f_R = 0$: isotropic radiation, $f_R = 1$: anisotropic radiation propagating in one direction). This expression is:

$$\chi_{\text{grey}} = \frac{3 + 4f_g^2}{5 + 2\sqrt{4 - 3f_g^2}} \; . \qquad (8)$$

We will refer to this as the *grey case Eddington factor*.

#### 2.1.2. M1-multigroup model

By cutting it into groups, the M1-multigroup takes into account more precisely the spectral behaviour in the electromagnetic spectrum. We can thus define G groups, each having bound frequencies [$v_g, v_{g+1}$]. To each group are associated radiative quantities of each group defined as:

$$\forall g \in [1, G], \left\{ \begin{aligned} \mathbb{E}_g &= \int_{v_g}^{v_{g+1}} \mathbb{E}_v \, d\nu \\ \mathbf{F}_g &= \int_{v_g}^{v_{g+1}} \mathbf{F}_v \, d\nu \\ \mathbb{P}_g &= \int_{v_g}^{v_{g+1}} \mathbb{P}_v \, d\nu \end{aligned} \right. \qquad (9)$$

Once again, Eq. (4) gives the closure relation of each group, and the Eddington tensor is related to the Eddington factor via Eq. (5). However, in this case, the Eddington factor lacks an analytical expression. To begin, we will determine the parameters on which it depends.

### 2.2. Dependencies of the Eddington factor in the M1-multigroup model

As seen in Eqs. (1), (2) and (3), we must first determine the expression of the specific intensity to compute the radiative quantities. By minimising the radiative entropy in the case of the M1-multigroup model [5], it can be expressed as:

$$\mathcal{I}_v(\Omega) = \sum_{g=1}^{G} \mathbb{I}_g(v) \mathcal{I}_{v,g}(\Omega) \; . \qquad (10)$$

Where $\mathbb{I}_g(v)$ is the indicator function, worth 1 for $v \in [v_g, v_{g+1}]$ and 0 everywhere else. $\mathcal{I}_{v,g}$ is the specific intensity of the group g and can be expressed as:

$$\mathcal{I}_{v,g}(\Omega) = \frac{2h\nu^3}{c^2} \left[ \exp\left(\frac{h\nu}{k_B}(\alpha_{0,g} + \alpha_{1,g}\mathbf{n}_g \cdot \mathbf{\Omega})\right) - 1 \right]^{-1} \qquad (11)$$

Where h and $k_B$ are, respectively, the Planck and Boltzmann constants, and $\mathbf{n}_g$ is the direction of the radiative flux of group g. If we introduce $\mathbf{n}_g \cdot \mathbf{\Omega} = \cos(\theta) = \mu$ and if we consider a frame in which $\theta = 0$ is the direction of the radiative flux, the radiative
<table><caption>Table 1: Variables used in this article</caption><thead><tr><th>Notation</th><th>Signification</th></tr></thead><tbody><tr><td><i>I</i><sub>v</sub></td><td>Spectral specific intensity (M1-multigroup model: Eq. 10)</td></tr><tr><td><i>I</i><sub>v,g</sub></td><td>Spectral specific intensity of group g (Eq. 11)</td></tr><tr><td><i>v</i><sub>g</sub>, <i>v</i><sub>g+1</sub></td><td>Frequency bounds of group g</td></tr><tr><td><i>v</i><sub>m,g</sub></td><td>Geometric mean of group g's frequency bounds, <i>v</i><sub>m,g</sub> = &radic;<span style="text-decoration: overline;"><i>v</i><sub>g</sub> <i>v</i><sub>g+1</sub></span></td></tr><tr><td><i>&delta;</i><sub>g</sub></td><td>Group g's frequency bounds' ratio, <i>&delta;</i><sub>g</sub> = <i>v</i><sub>g</sub>/<i>v</i><sub>g+1</sub> &isin; [0, 1]</td></tr><tr><td><i>E</i><sub>g</sub></td><td>Radiative energy of group g (Eq. 9)</td></tr><tr><td><b>F</b><sub>g</sub></td><td>Radiative flux vector of group g (Eq. 9). Components: <b>F</b><sub>g,x</sub>, <b>F</b><sub>g,y</sub></td></tr><tr><td><b>P</b><sub>g</sub></td><td>Radiative pressure tensor of group g (Eq. 9). Components: <b>P</b><sub>g,xx</sub>, <b>P</b><sub>g,yy</sub>, <b>P</b><sub>g,xy</sub></td></tr><tr><td><i>P</i><sub>g</sub></td><td>Component of the radiative pressure of group g in the direction of the radiative flux, <i>P</i><sub>g</sub> = <i>&chi;</i><sub>g</sub><i>E</i><sub>g</sub></td></tr><tr><td><i>&chi;</i><sub>grey</sub>, <i>&chi;</i><sub>g</sub></td><td>Eddington factor of the M1-grey model (Eq. 8) and of group g in the M1-multigroup model</td></tr><tr><td><i>f</i><sub>g</sub></td><td>Reduced flux of group g, <i>f</i><sub>g</sub> = ||<b>F</b><sub>g</sub>||/(c <i>E</i><sub>g</sub>) &isin; [0, 1]</td></tr><tr><td><i>T</i><sub>g</sub></td><td>Radiative temperature of group g, <i>T</i><sub>g</sub> = (<i>E</i><sub>g</sub>/<i>a</i><sub>R</sub>)<sup>1/4</sup></td></tr><tr><td><i>T</i><sub>g</sub><sup><i>v</i><sub>g+1</sub></sup></td><td>Adimensionalised radiative temperatures using <i>v</i><sub>m,g</sub> (Eq. 21)</td></tr><tr><td><i>T</i><sub>g</sub><sup><i>v</i><sub>g</sub></sup></td><td>Adimensionalised radiative temperatures using <i>v</i><sub>g+1</sub> (Eq. 23)</td></tr><tr><td><i>T</i><sub>g</sub><sup><i>v</i><sub>g</sub></sup></td><td>Adimensionalised radiative temperatures using <i>v</i><sub>g</sub> (Eq. 22)</td></tr><tr><td><i>T</i><sub>g</sub><sup><i>s</i></sup></td><td>Shifted radiative temperature (see Appendix C)</td></tr><tr><td><i>&alpha;</i><sub>0,g</sub>, <i>&alpha;</i><sub>0,g</sub></td><td>1<sup>st</sup> Lagrange multiplier and its adimensionalised form</td></tr><tr><td><i>&alpha;</i><sub>1,g</sub>, <i>&alpha;</i><sub>1,g</sub></td><td>2<sup>nd</sup> Lagrange multiplier and its adimensionalised form</td></tr><tr><td><i>&beta;</i><sub>g</sub></td><td>Anisotropy parameter, <i>&beta;</i><sub>g</sub> = <i>&alpha;</i><sub>1,g</sub>/<i>&alpha;</i><sub>0,g</sub> &isin; ] - 1, 1[</td></tr></tbody></table>

quantities of group g are expressed as:

$$E_g = \frac{2\pi}{c} \int_{v_g}^{v_{g+1}} \int_{-1}^{1} \frac{2h\nu^3/c^2}{e^{h\nu/k_B(\alpha_{0,g}+\alpha_{1,g}\mu)} - 1} \mathrm{d}\nu \mathrm{d}\mu ; \quad (12)$$

$$\|\mathbf{F}_g\| = 2\pi \int_{v_g}^{v_{g+1}} \int_{-1}^{1} \frac{2h\nu^3\mu/c^2}{e^{h\nu/k_B(\alpha_{0,g}+\alpha_{1,g}\mu)} - 1} \mathrm{d}\nu \mathrm{d}\mu ; \quad (13)$$

$$\mathcal{P}_g = \frac{2\pi}{c} \int_{v_g}^{v_{g+1}} \int_{-1}^{1} \frac{2h\nu^3\mu^2/c^2}{e^{h\nu/k_B(\alpha_{0,g}+\alpha_{1,g}\mu)} - 1} \mathrm{d}\nu \mathrm{d}\mu . \quad (14)$$

Where $\mathcal{P}_g$ is an eigenvector of the radiative pressure tensor and is related to the Eddington factor and the radiative energy as $\mathcal{P}_g = \chi_g E_g$. $\|\mathbf{F}_g\|$ is the norm of the radiative flux and $\alpha_{0,g}$ and $\alpha_{1,g}$ are called the *Lagrange multipliers* of the group g. It has been demonstrated that if $(E_g, \|\mathbf{F}_g\|)$ is physically realistic, then there exist a unique couple $(\alpha_{0,g}, \alpha_{1,g})$ [5], which implies that the Eddington factor depends only on the radiative energy $E_g$, the norm of the radiative flux $\|\mathbf{F}_g\|$ and the frequency bounds $v_g$ and $v_{g+1}$. However, in the M1-grey model, we observed that the Eddington factor only depends on the reduced flux $f_g$, indicating that further analysis is required to identify the parameters significantly impacting the Eddington factor.

Supposing that $v_g > 0$ and $v_{g+1} < \infty$ if we do the change of variable $\tilde{\nu} = \nu/v_{m,g}$, where $v_{m,g}$ is the geometric mean of the group's frequency bounds. We obtain the following expres-

sions:

$$\frac{k_B^4 E_g}{h^4 \nu_{m,g}^4} = a_R \frac{15}{\pi^4} \int_{\delta_g^{1/2}}^{\delta_g^{-1/2}} \int_{-1}^{1} \frac{\tilde{\nu}^3}{e^{\tilde{\nu}(\alpha_{0,g}+\alpha_{1,g}\mu)} - 1} \mathrm{d}\nu \mathrm{d}\mu ; \quad (15)$$

$$\frac{k_B^4 \|\mathbf{F}_g\|}{h^4 \nu_{m,g}^4} = a_R c \frac{15}{\pi^4} \int_{\delta_g^{1/2}}^{\delta_g^{-1/2}} \int_{-1}^{1} \frac{\tilde{\nu}^3 \mu}{e^{\tilde{\nu}(\alpha_{0,g}+\alpha_{1,g}\mu)} - 1} \mathrm{d}\tilde{\nu} \mathrm{d}\mu ; \quad (16)$$

$$\frac{k_B^4 \mathcal{P}_g}{h^4 \nu_{m,g}^4} = a_R \frac{15}{\pi^4} \int_{\delta_g^{1/2}}^{\delta_g^{-1/2}} \int_{-1}^{1} \frac{\tilde{\nu}^3 \mu^2}{e^{\tilde{\nu}(\alpha_{0,g}+\alpha_{1,g}\mu)} - 1} \mathrm{d}\mu . \quad (17)$$

Where $\tilde{\alpha}_{i,g} = \alpha_{i,g} h\nu_{m,g}/k_B$ are adimensionalised Lagrange multipliers, $\delta_g = v_g/v_{g+1}$ is the the group's frequency bounds' ratio, and $a_R$ is the radiative constant. We then can express the radiative temperature, the reduced flux and the Eddington factor as:

$$\frac{k_B T_g}{h \nu_{m,g}} = \left( \frac{15}{\pi^4} \int_{\delta_g^{1/2}}^{\delta_g^{-1/2}} \int_{-1}^{1} \frac{\tilde{\nu}^3}{e^{\tilde{\nu}(\alpha_{0,g}+\alpha_{1,g}\mu)} - 1} \mathrm{d}\nu \mathrm{d}\mu \right)^{1/4} ; \quad (18)$$

$$f_g = \frac{\int_{\delta_g^{1/2}}^{\delta_g^{-1/2}} \int_{-1}^{1} \frac{\tilde{\nu}^3 \mu}{e^{\tilde{\nu}(\alpha_{0,g}+\alpha_{1,g}\mu)} - 1} \mathrm{d}\tilde{\nu} \mathrm{d}\mu}{\int_{\delta_g^{1/2}}^{\delta_g^{-1/2}} \int_{-1}^{1} \frac{\tilde{\nu}^3}{e^{\tilde{\nu}(\alpha_{0,g}+\alpha_{1,g}\mu)} - 1} \mathrm{d}\tilde{\nu} \mathrm{d}\mu} ; \quad (19)$$

$$\chi_g = \frac{\int_{\delta_g^{1/2}}^{\delta_g^{-1/2}} \int_{-1}^{1} \frac{\tilde{\nu}^3 \mu^2}{e^{\tilde{\nu}(\alpha_{0,g}+\alpha_{1,g}\mu)} - 1} \mathrm{d}\tilde{\nu} \mathrm{d}\mu}{\int_{\delta_g^{1/2}}^{\delta_g^{-1/2}} \int_{-1}^{1} \frac{\tilde{\nu}^3}{e^{\tilde{\nu}(\alpha_{0,g}+\alpha_{1,g}\mu)} - 1} \mathrm{d}\tilde{\nu} \mathrm{d}\mu} . \quad (20)$$
From these expressions, it becomes evident that the Eddington factor in the M1-multigroup model is solely dependent on three quantities:

1. the adimensionalised radiative temperature $\mathcal{T}_g = k_B T_g/(h \nu_{m,g})$: comparison of the group's radiation energy to the energy level of the photons,

2. the reduced flux $f_g$: degree of anisotropy of the radiation in group g. Its values are in [0, 1] ($f_g = 0$: isotropic radiation, $f_g = 1$: directional radiation),

3. the frequency bounds' ratio $\delta_g = \nu_g/\nu_{g+1}$: indication of the separation of the group's frequency bounds. Its values are in [0, 1] ($\delta_g \to 0$: frequency bounds far from each other, $\delta_g \to 1$: frequency bounds very close to each other).

One can show in the same way that we can choose the adimensionalised temperature between the three following possibilities<fnref n="1" />:

$$ \mathcal{T}_g = k_B T_g/(h \nu_{m,g}) ; \qquad (21) $$

$$ \mathcal{T}_g^{\nu_g} = k_B T_g/(h \nu_g) = \mathcal{T}_g / \sqrt{\delta_g} ; \qquad (22) $$

$$ \mathcal{T}_g^{\nu_{g+1}} = k_B T_g/(h \nu_{g+1}) = \mathcal{T}_g \sqrt{\delta_g} . \qquad (23) $$

Depending on the context, it may be more appropriate to use one definition of the adimensionalised temperature over another. Therefore, we will always specify which adimensionalised radiative temperature we are using. For simplicity, we will refer to the adimensionalised radiative temperature simply as *radiative temperature* $\mathcal{T}$.

In Appendix H, we detail the cases $\nu_g = 0$ and $\nu_{g+1} \to \infty$. To efficiently compute the Eddington factor, we will analyse the significance of its dependence on all three parameters.

### 2.2.1. Influence of the reduced flux

First, we display curves of the Eddington factor $\chi_g$ as a function of the reduced flux $f_g$ in Fig. 1, from which we draw two conclusions:

1. Whatever $\mathcal{T}$, $\chi_g$ always evolves from $1/3$ at $f_g = 0$ to 1 at $f_g = 1$ and has a slope equal to 0 at $f_g = 0$,

2. $\mathcal{T}$ only influences the evolution of $\chi_g$ at intermediate values of $f_g$.

This shows that $\chi_g$ always depends strongly on $f_g$. Let us further analyse the dependence on $\mathcal{T}$ and $\delta_g$.

### 2.2.2. Influence of the radiative temperature and the frequency bounds' ratio

The ratio of the Eddington factor $\chi_g$ in the M1-multigroup model relative to its grey case expression is displayed as a function of the radiative temperature $\mathcal{T}_g$ using different frequency bounds' ratios $\delta_g$ (Fig. 2). Qualitatively when the group's frequencies are significantly different ($\delta_g \ll 1$), we observe the

<figure><img src="image_3.png" /><figcaption>Figure 1: Evolution of the Eddington factor $\chi_g$ at different radiative temperatures $\mathcal{T}_g$ (Eq. 21), using a frequency bounds' ratio worth $\delta_g = 10^{-4}$.</figcaption></figure>

existence of a central “plateau” lying in the domain (Fig. 2a). As the frequencies get closer together, the size of this “plateau” diminishes until it disappears (Fig. 2b, 2c). When the frequencies are nearly identical, two distinct domains emerge: one at low and one at high temperatures. The transition between these domains occurs at lower radiative temperatures as the frequencies become closer ($\delta_g \to 1$, Fig. 2d). Overall, three main domains can be seen in Fig. 2:

1. **Domain L:** Low radiative temperature. In this domain, $\mathcal{T}_g \ll \sqrt{\delta_g}$. $\chi_g$ depends solely on $f_g$. A method to compute it is detailed in Appendix A,

2. **Domain I:** Intermediate radiative temperatures. $\chi_g$ depends strongly on the $\mathcal{T}_g$, $f_g$ and $\delta_g$,

3. **Domain H:** High radiative temperature. In this domain, $\mathcal{T}_g \gg 1/\sqrt{\delta_g}$. $\chi_g$ depends solely on $f_g$. A method to compute it is detailed in Appendix B.

The intermediate domain I, based on the separation of the group's frequencies, can be further categorised as follows:

* **Case $\delta_g \to 0$:** Bound frequencies far from each other. There exist three intermediate domains (Fig. 2a):

- $I_a$: $\chi_g$ strongly depends on $\mathcal{T}$. Empirically, we have noticed $\chi_g$ only depends on $f_g$ and $\mathcal{T}_g^{\nu_g}$,

- $I_b$: corresponds to the domain where the grey case Eddington factor is valid. The width of this domain is related to $\delta_g$,

- $I_c$: $\chi_g$ strongly depends on $\mathcal{T}$. $\chi_g$ only depends on $f_g$ and $\mathcal{T}_g^{\nu_{g+1}}$,

* **Intermediate $\delta_g$:** In the domains $I_d$ and $I_e$, $\chi_g$ depends on $\mathcal{T}_g$, $f_g$ and $\delta_g$ (Fig. 2b, 2c),

* **Case $\delta_g \to 1$:** Bound frequencies close to each other. In the domain $I_f$, we observed that $\chi_g$ depends only on $f_g$ and a *shifted radiative temperature* $\mathcal{T}_g^s$ (Fig. 2d). This shift accounts for the transition domains L to H moving towards lower radiative temperatures (see Appendix C).

We have chosen to use neural networks to efficiently compute the Eddington factor $\chi_g$ in domains where it depends on two or more parameters. In the next section, we detail our computation methodology.

<footnote n="1">¹The choice of $\mathcal{T}_g$ comes from the calculations we present in this paper. The choice of $\mathcal{T}_g^{\nu_g}$ and $\mathcal{T}_g^{\nu_{g+1}}$ comes from analogous calculations, taking $\bar{v} = v/\nu_g$ and $\bar{v} = v/\nu_{g+1}$, respectively.</footnote>
<figure><img src="image_4.png" /><figcaption>(a) $\delta_g = 1.10^{-4}$. L, H: Asymptotic domains, $\chi_g = \mathcal{F}(f_g)$; $I_a$: Low-intermediate radiative temperatures, $\chi_g = \mathcal{F}(f_g, T_g^{\gamma_g})$; $I_b$: Domain in which $\chi_{grey}$ is valid; $I_c$: High-intermediate radiative temperatures, $\chi_g = \mathcal{F}(f_g, T_g^{\gamma_g+1})$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) $\delta_g = 7.10^{-2}$. L, H: Asymptotic domains, $\chi_g = \mathcal{F}(f_g)$; $I_d$: Low-intermediate radiative temperatures, $\chi_g = \mathcal{F}(f_g, T_g^{\gamma_g+1})$; $I_e$: High-intermediate radiative temperatures, $\chi_g = \mathcal{F}(f_g, T_g^{\gamma_g+1}, \delta_g)$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) $\delta_g = 7.10^{-2}$. L, H: Asymptotic domains, $\chi_g = \mathcal{F}(f_g)$; $I_d$: Low-intermediate radiative temperatures, $\chi_g = \mathcal{F}(f_g, T_g^{\gamma_g+1})$; $I_e$: High-intermediate radiative temperatures, $\chi_g = \mathcal{F}(f_g, T_g^{\gamma_g+1}, \delta_g)$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(d) $\delta_g = 0.99$. L, H: Asymptotic domains, $\chi_g = \mathcal{F}(f_g)$; $I_f$: Intermediate radiative temperatures, $\chi_g = \mathcal{F}(f_g, T_g^{\gamma_g})$. $T_g^s$ is the shifted radiative temperature (see Appendix C).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 3: Multi-Layer Perceptron comprising 5 inputs and 3 outputs. It contains one input layer, one hidden layer and one output layer. The circles represent perceptrons.</figcaption></figure>

# 3. Details on the AI method

## 3.1. A short introduction to supervised neural networks

In this work, we use the neural network architecture Multi-Layer Perceptron (MLP). It is composed of units called perceptrons. Each of these perceptrons takes a number $n$ of inputs $x_1, x_2, ..., x_n$ to which are associated weights $w_1, w_2, ..., w_n$ and a bias $b$ and returns an output $y = f(\sum_{i=1}^n w_i x_i + b)$, where $f$ is called activation function and is a non-linear function.

Several perceptrons can then be arranged in layers. There are different types of layers: the input layer receives the input data, the output layer generates the output data, and the hidden

layers are the layers between the input and the output layers (see Fig. 3).

The first thing one must do is to define the architecture of the neural network: the number of layers, the number of perceptrons, the activation function functions used at each layer, ... Once this is done, one must train the neural network for it to solve a given problem. This is done by finding the optimal values of the weights and biases of each perceptron in the neural network so that it outputs the correct answers to several known examples. This requires choosing a function that will evaluate the error of the neural network, which is called the cost function. Finally, one must select the optimiser that updates the weights and biases that minimise the cost function. The number of iteration steps during the optimisation phase is called epoch.

One splits the data into two sets: the training set, which is the set of known examples on which one trains the neural network and the test set, that enables to check if the neural network has correctly learnt.

## 3.2. Limits of the domains

As we have seen in Sec. 2, three main domains exist: the low radiative temperature domain L, several intermediate domains I and the high radiative temperature domain H. There are different sub-domains in the domain I, in which the Eddington factor depends more or less on the radiative temperature and the frequency bounds (domains $I_a, I_b, I_c, I_d, I_e$ and $I_f$).
In machine learning, there are two possible approaches: using a large neural network to predict the Eddington factor for all cases or using smaller, specialised neural networks to compute the Eddington factor in specific domains. We have chosen the latter approach to limit the size of the neural networks required, thereby reducing computation time. Table 2 summarises the methods used in the different domains, discussed in Sec. 1. Here are the names of the neural networks used for specific domains:

* MLP<sub>a</sub>: For domain I<sub>a</sub>. It has two inputs,

* MLP<sub>c</sub>: For domain I<sub>c</sub>. It has two inputs,

* MLP<sub>e,1</sub> and MLP<sub>e,2</sub>: For domain I<sub>e</sub>. We use two neural networks to reduce the size required and thus the prediction time. MLP<sub>e,1</sub> is used for δ<sub>g</sub> ∈ ]0.03, 0.1] and MLP<sub>e,2</sub> for δ<sub>g</sub> ∈ ]0.1, 0.9]. These two neural networks have three inputs,

* MLP<sub>f</sub>: For domain I<sub>f</sub>. It has two inputs.

For domain I<sub>d</sub>, we use MLP<sub>a</sub> that only takes two parameters and predicts accurate Eddington factors. For all these neural networks, we use normalised inputs and impose known constraints on the output (Appendix D).

<table><thead><tr><th rowspan="2"></th><th colspan="2">Range δ<sub>g</sub></th><th colspan="2">Range T<sub>g</sub></th><th rowspan="2">Method<sup>**</sup></th></tr><tr><th>Min</th><th>Max</th><th>Min</th><th>Max</th></tr></thead><tbody><tr><td>L</td><td>0<br>10<sup>-2.6</sup><br>0.9</td><td>10<sup>-2.6</sup><br>0.9<br>1</td><td>-∞<br>-∞<br>-∞<sup>*</sup></td><td>10<sup>-5</sup>√δ<sub>g</sub><br>1<br>10<sup>-3</sup>/√δ<sub>g</sub><sup>*</sup></td><td>P<br>In: f<sub>g</sub></td></tr><tr><td>H</td><td>0<br>10<sup>-2.6</sup><br>0.9</td><td>10<sup>-2.6</sup><br>0.9<br>1</td><td>10<sup>0.6</sup>/√δ<sub>g</sub><br>1<br>10<sup>0.6</sup>/√δ<sub>g</sub><sup>*</sup></td><td>+∞<br>+∞<br>+∞<sup>*</sup></td><td>P<br>In: f<sub>g</sub></td></tr><tr><td>I<sub>a</sub></td><td>0</td><td>0.03</td><td>10<sup>-5</sup>√δ<sub>g</sub></td><td>10<sup>1.1</sup>√δ<sub>g</sub></td><td>MLP<br>In: f<sub>g</sub>, T<sub>g</sub><sup>y<sub>g</sub></sup></td></tr><tr><td>I<sub>b</sub></td><td>0</td><td>10<sup>-2.6</sup></td><td>10<sup>1.1</sup>√δ<sub>g</sub></td><td>10<sup>-1.5</sup>/√δ<sub>g</sub></td><td>χ<sub>grey</sub></td></tr><tr><td>I<sub>c</sub></td><td>0</td><td>0.03</td><td>10<sup>-3</sup>/√δ<sub>g</sub></td><td>10<sup>0.6</sup>/√δ<sub>g</sub></td><td>MLP<br>In: f<sub>g</sub>, T<sub>g</sub><sup>y<sub>g+1</sub></sup></td></tr><tr><td>I<sub>d</sub></td><td>0.03</td><td>0.9</td><td>10<sup>-5</sup>/√δ<sub>g</sub></td><td>10<sup>-3</sup>/√δ<sub>g</sub></td><td>MLP<br>In: f<sub>g</sub>, T<sub>g</sub><sup>y<sub>g+1</sub></sup></td></tr><tr><td>I<sub>e</sub></td><td>0.03</td><td>0.9</td><td>10<sup>-3</sup>/√δ<sub>g</sub></td><td>10<sup>0.6</sup>/√δ<sub>g</sub></td><td>MLP<br>In: f<sub>g</sub>, δ<sub>g</sub>, T<sub>g</sub><sup>y<sub>g+1</sub></sup></td></tr><tr><td>I<sub>f</sub></td><td>0.9</td><td>1</td><td>10<sup>-3</sup>/√δ<sub>g</sub><sup>*</sup></td><td>10<sup>0.6</sup>/√δ<sub>g</sub><sup>*</sup></td><td>MLP<br>In: f<sub>g</sub>, T<sub>g</sub><sup>s</sup></td></tr></tbody></table>

<caption>Table 2: Limits for each domain.</caption>

### 3.3. Training/validation data

We have generated training and test datasets using the line search algorithm we have developed (Appendix G). Table 3 summarises the parameters used to generate these datasets. Linear spacing has been used for all the quantities. We draw the reader's attention to the fact that the amount of data expected does not correspond to the amount of data obtained because the line search algorithm does not always converge.

<table><thead><tr><th></th><th></th><th>Training set</th><th>Test set</th><th>Range</th></tr></thead><tbody><tr><td rowspan="3">MLP<sub>a</sub></td><td>log<sub>10</sub>(T<sub>g</sub><sup>y<sub>g</sub></sup>)</td><td>500</td><td>1,000</td><td>[-5, 1.1]</td></tr><tr><td>f<sub>g</sub></td><td>100</td><td>1,000</td><td>[0, 1]</td></tr><tr><td># data</td><td>49,695</td><td>993,372</td><td>-</td></tr><tr><td rowspan="3">MLP<sub>c</sub></td><td>log<sub>10</sub>(T<sub>g</sub><sup>y<sub>g+1</sub></sup>)</td><td>300</td><td>1,000</td><td>[-1.5, 0.6]</td></tr><tr><td>f<sub>g</sub></td><td>100</td><td>1,000</td><td>[0, 1]</td></tr><tr><td># data</td><td>28,621</td><td>951,821</td><td>-</td></tr><tr><td rowspan="3">MLP<sub>e,1</sub></td><td>log<sub>10</sub>(T<sub>g</sub><sup>y<sub>g+1</sub></sup>)</td><td>100</td><td>150</td><td>[-5, 0.6]</td></tr><tr><td>f<sub>g</sub></td><td>100</td><td>100<sup>*</sup></td><td>[0, 1]</td></tr><tr><td>δ<sub>g</sub></td><td>60</td><td>100</td><td>[0.03, 0.1]</td></tr><tr><td rowspan="3">MLP<sub>e,2</sub></td><td># data</td><td>590,160</td><td>1,467,561</td><td>-</td></tr><tr><td>log<sub>10</sub>(T<sub>g</sub><sup>y<sub>g+1</sub></sup>)</td><td>100</td><td>150</td><td>[-5, 1.1]</td></tr><tr><td>f<sub>g</sub></td><td>100</td><td>100<sup>*</sup></td><td>[0, 1]</td></tr><tr><td rowspan="3">MLP<sub>f</sub></td><td>δ<sub>g</sub></td><td>60</td><td>100</td><td>[0.1, 0.9]</td></tr><tr><td># data</td><td>589,910</td><td>1,465,462</td><td>-</td></tr><tr><td>log<sub>10</sub>(T<sub>g</sub><sup>s</sup>)</td><td>300</td><td>1,000</td><td>[-3, 0.6]</td></tr><tr><td rowspan="3">MLP<sub>f</sub></td><td>f<sub>g</sub></td><td>100</td><td>1,000</td><td>[0, 1]</td></tr><tr><td># data</td><td>28,712</td><td>955,143</td><td>-</td></tr></tbody></table>

<caption>Table 3: Description of the datasets.</caption>

### 3.4. Architectures of the neural networks

All neural networks were trained using Julia's Flux library<fnref n="2" />, known for its versatility in neural network modelling. We selected architectures large enough to provide precise predictions yet small enough to ensure quick computations. The architecture of these MLPs is detailed in Table 4. Our selection methodology is explained in Appendix E.

### 3.5. Validation

Finally, we check the relative error of the MLPs' prediction on the test set (see Sec. 3.3). Table 5 summarises the average relative error on the Eddington factor and the components of the radiative pressure. These errors are low, which demonstrates the generalisability of these MLPs. The relative error can reach

<footnote>*Min/Max values for the shifted temperature T<sub>g</sub><sup>s</sup> (see Appendix C).</footnote>

<footnote>**P: polynomial, MLP: neural network, χ<sub>grey</sub>: grey case Eddington factor. The variables after “In:” are the inputs used for each method.</footnote>

<footnote>*The range f<sub>g</sub> ∈ [10<sup>-5</sup>, 0.999999] has instead been used in the test sets of MLP<sub>e,1</sub> and MLP<sub>e,2</sub>, to change the values of f<sub>g</sub>.</footnote>

<footnote n="2">²https://fluxml.ai/Flux.jl/stable/</footnote>
<table><caption>Table 4: Architecture of the MLPs used.</caption><thead><tr><th>Neural network</th><th>Total # of layers</th><th>Total # of perceptrons</th><th># of parameters</th></tr></thead><tbody><tr><td>MLP<sub>a</sub></td><td>3</td><td>15</td><td>77</td></tr><tr><td>MLP<sub>c</sub></td><td>3</td><td>15</td><td>77</td></tr><tr><td>MLP<sub>e,1</sub></td><td>4</td><td>24</td><td>146</td></tr><tr><td>MLP<sub>e,2</sub></td><td>3</td><td>19</td><td>116</td></tr><tr><td>MLP<sub>f</sub></td><td>3</td><td>15</td><td>77</td></tr></tbody></table>

100 % on the radiative pressure component $P_{g,xy}$ for very low values of $f_g$, which is not an issue since this component does not have a significant impact on the physics in such domain, due to its low value compared to the other components of the radiative pressure.

<table><caption>Table 5: Average relative errors of each MLP on the test set. The average errors on the components of the radiative pressure are calculated using Eqs. (F.7), (F.8) and (F.9) in Appendix F.</caption><thead><tr><th></th><th>χ<sub>g</sub></th><th>P<sub>g,xx,yy</sub></th><th>P<sub>g,xy</sub></th></tr></thead><tbody><tr><td>MLP<sub>a</sub></td><td>3.3 × 10<sup>-3</sup> %</td><td>6.9 × 10<sup>-3</sup> %</td><td>2.2 × 10<sup>-2</sup> %</td></tr><tr><td>MLP<sub>c</sub></td><td>8.8 × 10<sup>-3</sup> %</td><td>2.4 × 10<sup>-2</sup> %</td><td>3.7 × 10<sup>-2</sup> %</td></tr><tr><td>MLP<sub>e,1</sub></td><td>8.2 × 10<sup>-3</sup> %</td><td>1.8 × 10<sup>-2</sup> %</td><td>4.1 × 10<sup>-2</sup> %</td></tr><tr><td>MLP<sub>e,2</sub></td><td>2.0 × 10<sup>-2</sup> %</td><td>3.7 × 10<sup>-2</sup> %</td><td>9.8 × 10<sup>-2</sup> %</td></tr><tr><td>MLP<sub>f</sub></td><td>1.7 × 10<sup>-2</sup> %</td><td>3.9 × 10<sup>-2</sup> %</td><td>8.5 × 10<sup>-2</sup> %</td></tr></tbody></table>

### 3.6. Utilisation in the code HADES

To exploit the trained neural networks, we have converted our Julia models to Keras models<fnref n="3" /> and then translated these into C functions with the library keras2c<fnref n="4" />. The C functions are read by the code HADES, written in Fortran 90 (upgraded to version 2003).

## 4. Validation: 1D Marshak wave

To evaluate this new method's efficiency, we simulated a 1D Marshak wave and compared it to different methods for computing the Eddington factor.

### 4.1. Simulation description: Marshak wave

We consider a box of size $L = 0.1$ m for 1,000 cells, in which there is a fluid at rest, having a temperature of $T = 300$ K and a density of $\rho = 8.2 \times 10^{-9}$ kg/m<fnref n="3" />. On its left boundary, the fluid experiences heating from an imposed temperature of $T_l = 5780$ K, accompanied by an influx of radiative flux governed by the following expression:

$$ \forall g \in [1, G], \forall x \in \mathbb{R}^*, F_{g,x}(x) = c(E_g(x) - E_g(0)) ; \quad (24) $$

where G is the total number of groups. Initially, the radiation is in equilibrium with the gas. The gas under consideration

possesses an atomic mass of $m_a = 1$ u and an adiabatic index of $\gamma = 5/3$. The simulation spans a duration of 13.3 ns of physical time, during which various configurations of group numbers, as detailed in Table 6, are employed. To compute the Eddington factor, we considered four methods:

1. the line search algorithm described in Appendix G,

2. using the grey case Eddington factor (Eq. 8),

3. the AI method we have described in this paper,

4. the interpolation of the Eddington factor from precalculated values [1].

<table><caption>Table 6: Definition of the groups</caption><thead><tr><th>Number of groups</th><th>ν<sub>g</sub> [eV]</th><th>ν<sub>g+1</sub> [eV]</th><th>δ<sub>g</sub></th></tr></thead><tbody><tr><td rowspan="2">2 groups</td><td>- group 1</td><td>0</td><td>1.0 × 10<sup>-1</sup></td></tr><tr><td>- group 2</td><td>1.0 × 10<sup>-1</sup></td><td>∞</td></tr><tr><td rowspan="3">3 groups</td><td>- group 1</td><td>0</td><td>1.0 × 10<sup>-3</sup></td></tr><tr><td>- group 2</td><td>1.0 × 10<sup>-3</sup></td><td>1.0 × 10<sup>-1</sup></td></tr><tr><td>- group 3</td><td>1.0 × 10<sup>-1</sup></td><td>∞</td></tr><tr><td rowspan="4">4 groups</td><td>- group 1</td><td>0</td><td>1.0 × 10<sup>-3</sup></td></tr><tr><td>- group 2</td><td>1.0 × 10<sup>-3</sup></td><td>2.0 × 10<sup>-1</sup></td></tr><tr><td>- group 3</td><td>2.0 × 10<sup>-1</sup></td><td>4.0 × 10<sup>-1</sup></td></tr><tr><td>- group 4</td><td>4.0 × 10<sup>-1</sup></td><td>∞</td></tr><tr><td rowspan="5">5 groups</td><td>- group 1</td><td>0</td><td>1.0 × 10<sup>-3</sup></td></tr><tr><td>- group 2</td><td>1.0 × 10<sup>-3</sup></td><td>1.0 × 10<sup>0</sup></td></tr><tr><td>- group 3</td><td>1.0 × 10<sup>0</sup></td><td>1.1 × 10<sup>0</sup></td></tr><tr><td>- group 4</td><td>1.1 × 10<sup>0</sup></td><td>1.0 × 10<sup>1</sup></td></tr><tr><td>- group 5</td><td>1.0 × 10<sup>1</sup></td><td>∞</td></tr></tbody></table>

### 4.2. Computation time of different methods

<table><caption>Table 7: Comparison of the simulation times (CPU time).</caption><thead><tr><th>Number of groups</th><th>Grey case expression</th><th>AI method</th><th>Interpolation</th><th>Line search</th></tr></thead><tbody><tr><td>2 groups</td><td>3.11 min</td><td>4.18 min</td><td>3.70 min</td><td>3.58 days</td></tr><tr><td>3 groups</td><td>3.92 min</td><td>4.92 min</td><td>4.69 min</td><td>7.05 days</td></tr><tr><td>4 groups</td><td>4.55 min</td><td>10.37 min</td><td>5.65 min</td><td>22.6 days</td></tr><tr><td>5 groups</td><td>5.53 min</td><td>7.72 min</td><td>6.88 min</td><td>16.6 days</td></tr></tbody></table>

Table 7 presents a comprehensive comparison of simulation times across different multigroup configurations. Our developed method demonstrates remarkable efficiency, outperforming the line search algorithm by 1,000 to 3,000 times in computation time. However, it is slightly slower than the method using the grey case expression, i.e. the fastest method, by a factor of 1 to 3, its computation speed remains comparable to the fastest methods available.

It is worth noting that simulation times for the 4-group configuration are prolonged compared to other setups. This delay is attributed to utilising heavier neural networks, as seen in Sec. 3.4,

<footnote n="3">³https://keras.io</footnote>

<footnote n="4">⁴https://f0uriest.github.io/keras2c/index.html</footnote>
which necessitates more time for predictions. Consequently, our developed method exhibits accelerated performance, particularly evident with fewer groups avoiding the use of MLP<sub>e,1</sub> and MLP<sub>e,2</sub>, which are heavier.

### 4.3. Accuracy

At time t=0.133 ns, we present the radiative pressure $P_{R,xx}$ obtained using the four different methods for computing the Eddington factor. The simulation using the line search method is designated as the reference due to its superior accuracy. Analysis of the resulting curves reveals that the AI method exhibits errors of $10^{-3}$ %, whereas the grey case expression can yield errors up to 2.5 %. Additionally, the error in radiative pressure can reach 0.1 %, depending on the set of precalculated Eddington factors used. This disparity highlights the enhanced precision of our new method. Furthermore, we observe that the error associated with the grey case expression and the interpolation method increases as the number of groups increases.

## 5. Conclusion

To conclude, we have obtained a neural network approach that provides an average error on the prediction of the Eddington factor in the M1-multigroup case of $10^{-3}$ % to $10^{-2}$ % and is faster by a factor 3,000 compared to a line search algorithm and only slower by a factor 1 to 3 compared using the grey case Eddington factor. Moreover, this approach does not need prior knowledge of the radiative energy and reduced flux, which is an advantage compared to the method using interpolations of precalculated Eddington factors, making it one of the most efficient existing.

Yet, the simulations we have done so far have not revealed significant differences in the hydrodynamics quantities between the different methods, even using many groups or 2D simulations. The impact of the grey case Eddington factor's lack of precision could yet prove decisive in the case where the radiation dominates, and the reduced flux is close to 1.

However, we propose this method to serve as a reference to future work using the M1-multigroup, as it renders precise estimation of the Eddington factor computationally accessible.

## Acknowledgements

This work was supported by a French government grant managed by the Agence Nationale de la Recherche under the Investissement d'avenir program, reference ANR-19-P3IA-0002.

## Declaration of generative AI and AI-assisted technologies in the writing process.

During the preparation of this work we used ChatGPT and Grammarly in order to improve the clarity of the article and correct the grammar. After using this tool/service, the author(s) reviewed and edited the content as needed and take(s) full responsibility for the content of the published article.

## Appendix A. Low-temperature approximation for the Eddington factor (Domain L)

We have fitted a polynomial of order 11 at $T_g^{v_g} = 10^{-8}$, using Julia's least square regression library LsqFit<sup>5</sup>. The resulting expression is:

$$ P_L(f_g) = 1/3 + 2f_g^2/3 + \begin{aligned} & 0.3058945350580 (f_g^3 - f_g^2) - \\ & 3.5960491961545 (f_g^4 - f_g^2) + \\ & 24.1130825450229 (f_g^5 - f_g^2) - \\ & 92.2832025267787 (f_g^6 - f_g^2) + \\ & 220.1853600902954 (f_g^7 - f_g^2) - \\ & 329.3727903123292 (f_g^8 - f_g^2) + \\ & 299.9957746394199 (f_g^9 - f_g^2) - \\ & 151.3384873479037 (f_g^{10} - f_g^2) + \\ & 32.2668624654447 (f_g^{11} - f_g^2) . \end{aligned} \quad (A.1) $$

The form of this polynomial ensures that the predicted Eddington factor is $1/3$ at $f_g = 0$, $1$ at $f_g = 1$ and that the slope is zero at $f_g = 0$. The maximal relative error is $10^{-5}$ % on the Eddington factor and $10^{-3}$ % on the radiative pressure components.

## Appendix B. High-temperature approximation for the Eddington factor (Domain H)

If we consider the Lagrange multiplier $\alpha_0$ to be $\alpha_0 \ll k_B/(h\nu_{g+1})$ the specific intensity of group g can express $I_{v,g}$ as:

$$ I_{v,g}(\mu) = \frac{2k_B}{c^2} \frac{\nu^2}{\alpha_{0,g}(1 + \beta_g\mu)} . \quad (B.1) $$

Where $\beta_g = \alpha_{1,g}/\alpha_{0,g}$. The Eqs. (12), (13) and (14) can be rewritten:

$$ E_g = \frac{8\pi k_B (\nu_{g+1}^3 - \nu_g^3) \operatorname{arctanh}(\beta_g)}{3c^3 \alpha_{0,g} \beta_g \beta_g} ; \quad (B.2) $$

$$ \|\mathbf{F}_g\| = - \frac{8\pi k_B (\nu_{g+1}^3 - \nu_g^3) \operatorname{arctanh}(\beta_g) - \beta_g}{3c^3 \alpha_{0,g} \beta_g^2} ; \quad (B.3) $$

$$ P_g = \frac{8\pi k_B (\nu_{g+1}^3 - \nu_g^3) \operatorname{arctanh}(\beta_g) - \beta_g}{3c^3 \alpha_{0,g} \beta_g^3} . \quad (B.4) $$

One can verify using these expressions that these equations imply that the radiative temperature is such that $T_g \gg h\nu_{g+1}/k_B$, i.e. $T_g^{v_{g+1}} \gg 1$. The reduced flux and Eddington factor can then be expressed:

$$ f_g = - \frac{\operatorname{arctanh}(\beta_g) - \beta_g}{\beta_g \operatorname{arctanh}(\beta_g)} ; \quad (B.5) $$

$$ \chi_g = \frac{\operatorname{arctanh}(\beta_g) - \beta_g}{\beta_g^2 \operatorname{arctanh}(\beta_g)} . \quad (B.6) $$

<footnote n="5">⁵https://julianlsolvers.github.io/LsqFit.jl/latest/</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 4: Component xx of the radiative pressure at time t=0.133 ns. The relative error is given by the formula: $(P_{xx,method} - P_{xx,line search})/P_{xx,line search} \times 100$.</figcaption></figure>

Thus, the Eddington factor only depends on the reduced flux. From these expressions, one can show that the Eddington factor is 1/3 at $f_g = 0$, 1 at $f_g = 1$ and that the slope of the Eddington factor is zero at $f_g = 0$, which is consistent with what we saw in Sec. 2.2.1. Additionally, at $f_g = 1$ the slope of the Eddington factor is worth $\partial_{f_g}\chi_g(1) = 1$.

Since these expressions do not provide an expression of $\chi_g$ as a function of $f_g$, We fitted a polynomial that incorporates the properties mentioned above. Using a dichotomy algorithm, we calculated 100 values of $\chi_g$ for a reduced flux varying between 0 and 1. Then, we fitted a polynomial $P_H$ of degree 12 using Julia's library LsqFit⁵. The resulting polynomial expression is:

$$ P_H(f_g) = 1/3 + f_g^2 - f_g^3/3 + \begin{aligned} & 6.5813886320063 (f_g^4 - 2f_g^3 + f_g^2) - \\ & 44.5593808930324 (f_g^5 - 3f_g^3 + 2f_g^2) + \\ & 179.9251066153469 (f_g^6 - 4f_g^3 + 3f_g^2) - \\ & 463.2145547920471 (f_g^7 - 5f_g^3 + 4f_g^2) + \\ & 776.0741150675088 (f_g^8 - 6f_g^3 + 5f_g^2) - \\ & 841.0048410765069 (f_g^9 - 7f_g^3 + 6f_g^2) + \\ & 566.6515584347650 (f_g^{10} - 8f_g^3 + 7f_g^2) - \\ & 215.1853394583736 (f_g^{11} - 9f_g^3 + 8f_g^2) + \\ & 35.1130440991836 (f_g^{12} - 10f_g^3 + 9f_g^2) . \end{aligned} \quad (B.7) $$

The maximum relative error of this expression is worth $10^{-5}\%$ on the Eddington factor and $10^{-3}\%$ on the radiative pressure components.

## Appendix C. Expression of the shifted radiative temperature (Case $\delta_g \to 1$)

For values of the bound frequencies that are very close ($\delta_g \ge 0.9$), we observed that the topology of the Eddington factor remains unchanged but shifts towards lower radiative temperatures (Fig. C.5). To determine the extent of this temperature shift, we examined how the evolution of the position of the most significant gradient, $\mathcal{T}_g^{\nabla,max}$ according to $\delta_g$. We see that $\log_{10}(\mathcal{T}_g^{\nabla,max})$ varies linearly with $X = \log_{10}(1 - \delta_g)$ as $\delta_g \to 1$ (Fig. C.6). Therefore, we fitted a linear function for $\delta_g > 0.99$ (i.e. $X < -2$) and obtained the following expression:

$$ \log_{10}(\mathcal{T}_g^{\nabla,max}) = 0.250114 X - 0.193613 . \quad (C.1) $$

For $\delta_g \in [0.9, 0.99]$ (i.e. $X \in [-2, -1]$), we fitted a second-degree polynomial as the relationship starts to deviate from linearity. We have obtained the following expression:

$$ \log_{10}(\mathcal{T}_g^{\nabla,max}) = 0.007385 X^2 + 0.277981 X - 0.166850 . \quad (C.2) $$

We introduce a *shifted radiative temperature*, using the radiative temperature at $\delta_g = 0.99$ as a reference (i.e., $\mathcal{T}_g^s$ equals $\mathcal{T}_g$ when $\delta_g = 0.99$). For $\delta_g > 0.99$, it is defined as:

$$ \log_{10}(\mathcal{T}_g^s) = \log_{10}(\mathcal{T}_g) - 0.250114 (X + 2) . \quad (C.3) $$

For $\delta_g \in [0.9, 0.99]$, it is defined as:

$$ \log_{10}(\mathcal{T}_g^s) = \log_{10}(\mathcal{T}_g) - 0.007385 (X^2 - 4) - 0.277981 (X + 2) . \quad (C.4) $$
<figure><img src="image_10.png" /><figcaption>Figure C.5: Shift of the curve as $\delta_g$ gets closer to 1.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure C.6: Relation of the radiative temperature at the position of the most significant gradient as a function of $X = \log_{10}(1 - \delta_g)$.</figcaption></figure>

## Appendix D. Input/Outputs of neural networks

As inputs for the neural networks, normalised values are used to enhance learning efficiency:

* Inputs for MLP<sub>a</sub>:

$$ 
\begin{aligned}
x_1 &= \frac{\log_{10}(\mathcal{T}_g^{\nu_g}) - \log_{10}(\mathcal{T}_{g,min}^{\nu_g})}{\log_{10}(\mathcal{T}_{g,max}^{\nu_g}) - \log_{10}(\mathcal{T}_{g,min}^{\nu_g})} ; \\
x_2 &= f_g .
\end{aligned}
$$

* Inputs for MLP<sub>c</sub>:

$$ 
\begin{aligned}
x_1 &= \frac{\log_{10}(\mathcal{T}_g^{\nu_{g+1}}) - \log_{10}(\mathcal{T}_{g,min}^{\nu_{g+1}})}{\log_{10}(\mathcal{T}_{g,max}^{\nu_{g+1}}) - \log_{10}(\mathcal{T}_{g,min}^{\nu_{g+1}})} ; \\
x_2 &= f_g .
\end{aligned}
$$

* Inputs for MLP<sub>e,1</sub> and MLP<sub>e,2</sub>:

$$ 
\begin{aligned}
x_1 &= \frac{\log_{10}(\mathcal{T}_g^{\nu_{g+1}}) - \log_{10}(\mathcal{T}_{g,min}^{\nu_{g+1}})}{\log_{10}(\mathcal{T}_{g,max}^{\nu_{g+1}}) - \log_{10}(\mathcal{T}_{g,min}^{\nu_{g+1}})} ; \\
x_2 &= f_g ; \\
x_3 &= \frac{\delta_g - \delta_{g,min}}{\delta_{g,max} - \delta_{g,min}} .
\end{aligned}
$$

* Inputs for MLP<sub>f</sub>:

$$ 
\begin{aligned}
x_1 &= \frac{\log_{10}(\mathcal{T}_g^s) - \log_{10}(\mathcal{T}_{g,min}^s)}{\log_{10}(\mathcal{T}_{g,max}^s) - \log_{10}(\mathcal{T}_{g,min}^s)} ; \\
x_2 &= f_g .
\end{aligned}
$$

Where the $X_{min}$ and $X_{max}$ quantities represent the minimum and maximum values within the considered domain, and $x_i$ is an input of the neural networks.

Additionally, the MLPs do not directly predict the Eddington factor $\chi_g$. Instead, they predict a related quantity, described by the expression:

$$ \chi_g = 1/3 + 2/3 f_g^2 (1 + (1 - f_g)(3/2 y - 1)) . \quad (\text{D.1}) $$

Where $y$ is the output of the MLP, which has been empirically verified to fall within the range [0, 1]. This technique ensures that the predicted Eddington factor is $1/3$ at $f_g = 0$, $1$ at $f_g = 1$, and that the slope is zero at $f_g = 0$.

## Appendix E. Selection of the architecture of the neural network

### Appendix E.1. Training strategy

To train the neural networks, we utilise Julia's library Flux<sup>6</sup>. We employ the Light Broyden-Fletcher-Goldfarb-Shanno (LBFGS) optimiser from the Julia library Optim<sup>7</sup>, chosen for its rapid convergence and efficient minimisation of the cost function. Although LBFGS lacks stochastic properties, potentially affecting generalisation, our comparison with the vanilla Adam optimiser [14] shows that neural networks trained with LBFGS remain generalisable in our specific scenarios. For the training, we have defined two distinct cost functions:

$$ \mathcal{L}_1(x_1, x_2, \chi_{g,\text{data}}, \chi_{g,\text{pred}}) = \left( \frac{\chi_{g,\text{pred}} - \chi_{g,\text{data}}}{\alpha_f \alpha_T} \right)^2 ; \quad (\text{E.1}) $$

$$ \mathcal{L}_2(x_1, x_2, x_3, \chi_{g,\text{data}}, \chi_{g,\text{pred}}) = \left( \frac{\chi_{g,\text{pred}} - \chi_{g,\text{data}}}{\alpha_f \alpha_T \alpha_v} \right)^2 . \quad (\text{E.2}) $$

Where $\chi_{g,\text{data}}$ and $\chi_{g,\text{pred}}$ are, respectively, the label and predicted values of the Eddington factor, $x_1$, $x_2$, $x_3$ are the inputs of the neural network (see Appendix D), and $\alpha_f$, $\alpha_T$, $\alpha_v$ are factors defined by the expressions:

$$ \alpha_T = \min(x_1, 1 - x_1) + \epsilon_T ; \quad (\text{E.3}) $$

$$ \alpha_f = \min(1 - x_2, \chi_{g,\text{data}} - 1/3) + \epsilon_f ; \quad (\text{E.4}) $$

$$ \alpha_v = \min(x_3, 1 - x_3) + \epsilon_v . \quad (\text{E.5}) $$

Where $\epsilon_T$, $\epsilon_f$ and $\epsilon_v$ are constants that prevent these factors to be worth 0. We have respectively set these constants to 0.1, $10^{-3}$ and 0.1. These factors serve two purposes:

<footnote n="6"><sup>6</sup>https://fluxml.ai/Flux.jl/stable/</footnote>

<footnote n="7"><sup>7</sup>https://julianlsolvers.github.io/Optim.jl/v0.9.3/algo/lbfgs/</footnote>
<figure><img src="image_12.png" /><figcaption>Figure E.7: Learning curves of the different architectures tested. L represents the number of hidden layers, and N is the total number of perceptrons. For example, architecture L1N6 has 1 hidden layer and 6 perceptrons.</figcaption></figure>

* $\alpha_T$ and $\alpha_v$ reduce the discontinuities at the interfaces of the domains defined in Sec. 3.2,

* $\alpha_f$ emphasises the importance of domains in which an error on the Eddington factor has a significant impact on the radiative pressure (see Appendix F).

We use $\mathcal{L}_1$ to train the neural networks MLP$_a$, MLP$_c$ and MLP$_f$ and $\mathcal{L}_2$ for MLP$_{e,1}$ and MLP$_{e,2}$.

The models are trained until the weights and biases of the MLP exhibit variations below $10^{-32}$; the predicted values show variations below $10^{-32}$, or the gradient falls below $10^{-8}$.

## Appendix E.2. Selection of the architecture

For all the neural networks, we use different activation functions:

* input layer: linear function,

* hidden layers: hyperbolic tangent function tanh,

* output layer: sigmoid function $\sigma$.

Let us consider the neural network MLP$_a$ as an example. We trained various neural networks with different numbers of hidden layers and perceptrons that are simple enough to ensure fast prediction times. We have used all the training set data specified in Sec. 3.3 to train the neural networks. The architecture we have selected is the one whose cost function is the lowest at the end of the training, which is the architecture having 3 hidden layers and 15 neurons in the case of MLP$_a$ (Fig. E.7). We have employed the same methodology for all the other neural networks.

## Appendix F. Relative error of the radiative pressure

Using Eq. (5), one can express each component of the radiative pressure tensor as:

$$ \mathbb{P}_{g,xx} = \left\{ \frac{1 - \chi_g}{2} + \frac{3\chi_g - 1}{2} \frac{F_{g,x}^2}{\|\mathbf{F}_g\|^2} \right\} E_g ; \quad (F.1) $$

$$ \mathbb{P}_{g,yy} = \left\{ \frac{1 - \chi_g}{2} + \frac{3\chi_g - 1}{2} \frac{F_{g,y}^2}{\|\mathbf{F}_g\|^2} \right\} E_g ; \quad (F.2) $$

$$ \mathbb{P}_{g,xy} = \frac{3\chi_g - 1}{2} \frac{F_{g,x}F_{g,y}}{\|\mathbf{F}_g\|^2} E_g . \quad (F.3) $$

Considering only an error $\delta\chi_g$ in the determination of the Eddington factor, The relative error on the radiative pressure components can be expressed:

$$ \left| \frac{\delta \mathbb{P}_{g,xx}}{\mathbb{P}_{g,xx}} \right| = \left| \frac{3x^2 - 1}{(3x^2 - 1)\chi_g + 1 - x^2} \right| \delta\chi_g ; \quad (F.4) $$

$$ \left| \frac{\delta \mathbb{P}_{g,yy}}{\mathbb{P}_{g,yy}} \right| = \left| \frac{3y^2 - 1}{(3y^2 - 1)\chi_g + 1 - y^2} \right| \delta\chi_g ; \quad (F.5) $$

$$ \left| \frac{\delta \mathbb{P}_{g,xy}}{\mathbb{P}_{g,xy}} \right| = \frac{\delta\chi_g}{\chi_g - 1/3} . \quad (F.6) $$

Where $x = F_{g,x}/\|\mathbf{F}_g\|$ and $y = F_{g,y}/\|\mathbf{F}_g\|$. $\|\delta\mathbb{P}_{g,xx}/\mathbb{P}_{g,xx}\|$ and $\|\delta\mathbb{P}_{g,yy}/\mathbb{P}_{g,yy}\|$ admit a maximum relative error when respectively $x = 0$ or $1$ and $y = 0$ or $1$. In this case, we can express the maximum relative errors as:

$$ \left| \frac{\delta \mathbb{P}_{g,xx}}{\mathbb{P}_{g,xx}} \right|_{\max} = \frac{\delta\chi_g}{\min\{\chi_g, 1 - \chi_g\}} ; \quad (F.7) $$

$$ \left| \frac{\delta \mathbb{P}_{g,yy}}{\mathbb{P}_{g,yy}} \right|_{\max} = \frac{\delta\chi_g}{\min\{\chi_g, 1 - \chi_g\}} ; \quad (F.8) $$

$$ \left| \frac{\delta \mathbb{P}_{g,xy}}{\mathbb{P}_{g,xy}} \right|_{\max} = \frac{\delta\chi_g}{\chi_g - 1/3} . \quad (F.9) $$

These expressions reveal the importance of accurate predictions of the Eddington factor at values of the reduced flux $f_g$ close to 0 or 1 to reduce the relative errors on the components of the radiative pressure tensor.

## Appendix G. A short description of the 2D line search algorithm developed

According to Turpault 2002 [5], it is established that a physically realistic pair ($E_g$, $f_g$) corresponds uniquely to a set of Lagrange multipliers ($\alpha_{0,g}$, $\alpha_{1,g}$). Therefore, using a search algorithm to derive these Lagrange multipliers from the radiative quantities is feasible. Instead of directly employing the pairs ($\alpha_{0,g}$, $\alpha_{1,g}$) and ($E_g$, $f_g$), we’ve opted for inputs ($x_1 = \log(\alpha_{0,g})$, $x_2 = 100\sigma^{-1}(\beta_g)$), where $\sigma^{-1}$ denotes the inverse of the sigmoid function, and outputs ($y_1 = \log(E_g)$, $y_2 = \log(f_g)$). This choice facilitates the utilisation of a line search algorithm and ensures quicker convergence. The computation of $\log(E_g)$ and $\log(f_g)$ follows the methodology outlined in Chinh Hyunh Nguyen’s
<figure><img src="image_13.png" /><figcaption>Figure H.8: Ratio $\chi_g/\chi_{grey}$ with $f_g = 0.65$ $f_g = 0.65$ and $v_g = 0$. $I_b$: Domain in which $\chi_{grey}$ is valid; $I_c$: High-intermediate radiative, $\chi_g = \mathcal{F}(f_g, \mathcal{T}_g^{v_g+1})$; H: Asymptotically high radiative temperatures, $\chi_g = \mathcal{F}(f_g)$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure H.9: Ratio $\chi_g/\chi_{grey}$ with $f_g = 0.65$ and $v_{g+1} \to \infty$. L: Asymptotically low radiative temperatures, $\chi_g = \mathcal{F}(f_g)$; $I_a$: Low-intermediate radiative temperatures, $\chi_g = \mathcal{F}(f_g, \mathcal{T}_g^{v_g})$; $I_b$: Domain in which $\chi_{grey}$ is valid.</figcaption></figure>

thesis [1]. As for the line search algorithm, we’ve implemented a Fortran version based on the code from Numerical Recipes for Fortran 90 [15]. The Lagrange multipliers are searched until the error on the radiative energy and the reduced flux is calculated with the Eqs. (12) and (19) are less than $10^{-6}$ %.

# Appendix H. Case where $v_g = 0$ or $v_{g+1} \to \infty$

## Appendix H.1. Case $v_g = 0$

In this scenario, the only well-defined adimensionalised radiative temperature is $\mathcal{T}_g^{v_g+1}$. Therefore, the Eddington factor is such that $\chi_g = \mathcal{F}(f_g, \mathcal{T}_g^{v_g+1})$. Fig. H.8 illustrates the Eddington factor’s behaviour as a function of $\mathcal{T}_g^{v_g+1}$. In this case, we retrieve domains $I_b$, $I_c$, and $H$ as discussed in the limit where $\delta_g \to 0$ (Sec. 2.2.2).

## Appendix H.2. Case $v_{g+1} \to \infty$

Here, the only well-defined adimensionalised radiative temperature is $\mathcal{T}_g^{v_g}$. Consequently, the Eddington factor is such that $\chi_g = \mathcal{F}(f_g, \mathcal{T}_g^{v_g})$. Figure H.9 illustrates the Eddington factor’s behaviour as a function of $\mathcal{T}_g^{v_g}$. In this case, we retrieve domains $L$, $I_a$, and $I_b$ as discussed in the limit where $\delta_g \to 0$ (Sec. 2.2.2).

# References

[1] H. C. Nguyen, Simulation de modèles hydrodynamiques et de transfert radiatif intervenant dans la description d'écoulements astrophysiques, Ph.D. thesis, Paris 11 (2011).

[2] C. Michaut, H. Nguyen, L. Di Menza, Computational radiation hydrodynamics, Astrophysics and Space Science 336 (1) (2011) 175-181.

[3] C. Michaut, L. Di Menza, H. Nguyen, S. Bouquet, M. Mancini, Hades code for numerical simulations of high-mach number astrophysical radiative flows, High Energy Density Physics 22 (2017) 77-89.

[4] B. Dubroca, J.-L. Feugeas, Etude théorique et numérique d'une hiérarchie de modèles aux moments pour le transfert radiatif, Comptes Rendus de l'Académie des Sciences - Series I - Mathematics 329 (10) (1999) 915-920.

[5] R. Turpault, Construction d'un modèle m1-multigroupe pour les équations du transfert radiatif, Comptes Rendus Mathematique 334 (4) (2002) 331-336.

[6] R. Turpault, A consistent multigroup model for radiative transfer and its underlying mean opacities, Journal of Quantitative Spectroscopy and Radiative Transfer 94 (3) (2005) 357-371.

[7] B. Freeman, L. Hauser, J. Palmer, S. Pickard, G. Simmons, D. Williston, J. Zerkle, Dasa report no. 2135, Vol. I (La Jolla: Systems, Science, and Software, Inc.) (1968).

[8] G. Spillman, Formulation of the eddington factor for use in one-dimensional nonequilibrium diffusion calculations, Los Alamos National Laboratory Office Memorandum, Los Alamos National Laboratory, Los Alamos, NM 87545 (1968) 25.

[9] D. Mihalas, B. Weibel-Mihalas, Foundations of radiation hydrodynamics, Oxford University Press, 1984.

[10] C. Levermore, Relating eddington factors to flux limiters, Journal of Quantitative Spectroscopy and Radiative Transfer 31 (2) (1984) 149-160.

[11] R. Turpault, Modelisation, approximation numerique et applications du transfert radiatif en desequilibre spectral couple avec l'hydrodynamique, Theses, Université Sciences et Technologies - Bordeaux I (Dec. 2003).

[12] N. Vaytet, E. Audit, B. Dubroca, F. Delahaye, A numerical model for multigroup radiation hydrodynamics, Journal of Quantitative Spectroscopy and Radiative Transfer 112 (8) (2011) 1323-1335.

[13] A. Harada, S. Nishikawa, S. Yamada, Deep learning of the eddington tensor in core-collapse supernova simulation, The Astrophysical Journal 925 (2) (2022) 117.

[14] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, CoRR abs/1412.6980 (2014).

[15] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical recipes in Fortran 90 (2nd ed.): the art of parallel scientific computing, Cambridge University Press, USA, 1996.