UNIVERSITÉ
CÔTE D'AZUR

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

ÉCOLE DOCTORALE
SCIENCES
FONDAMENTALES
ET APPLIQUÉES

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

$$p \left( \frac{\partial v}{\partial t} + v \cdot \nabla v \right) = - \nabla p + \nabla \cdot T + f$$

$$e^{i\pi} + 1 = 0$$

# THÈSE DE DOCTORAT

## Machine learning for radiative hydrodynamics in astrophysics

Gonzague Paul-Constantin Hugues-Antoine RADUREAU

Laboratoire J.-L. Lagrange

Translated by the author, based on the official version archived on theses.fr (ID 2025COAZ5041).

### Presented with a view to obtaining the degree of doctor in

Sciences de la Planète et de l'Univers d'Université Côte d'Azur

Supervised by: Claire MICHAUT

defended the: 30 Septembre 2025

### Before the jury composed of:

Jérôme Breil,

Research Director CEA, CEA-CESTA, Le Barp, France

Emmanuel Franck,

Research Scientist, INRIA, Strasbourg, France

Stéphane Mazevet,

Director of OCA, Nice, France

Florence Hubert,

Professor, I2M, Marseille, France

Serge Bouquet,

Emeritus Research Director CEA, Paris, France

Claire Michaut,

Research Director, OCA, Nice, France

Andrew Comport,

Research Scientist, I3S, Sophia-Antipolis, France

UNIVERSITÉ
CÔTE D'AZUR

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

3iA Côte d'Azur
Interdisciplinary Institute
for Artificial Intelligence
<figure><img src="image_4.png" /></figure>
MACHINE LEARNING FOR RADIATIVE
HYDRODYNAMICS IN ASTROPHYSICS

JURY:

CHAIR OF THE JURY

Stéphane Mazevet,
*Director of OCA, Nice, France*

REVIEWERS

Jérôme Breil,
*Research Director, CEA, Le Barp, France*

Emmanuel Franck,
*Research Scientist, INRIA, Strasbourg, France*

EXAMINERS

Florence Hubert,
*Professor, I2M, Marseille, France*

Serge Bouquet,
*Emeritus Research Director, CEA, Paris, France*

SUPERVISOR

Claire Michaut,
*Research Director, OCA, Nice, France*

CO-SUPERVISOR

Andrew Comport,
*Research Scientist, I3S, Sophia-Antipolis, France*
# Machine learning for radiative hydrodynamics in astrophysics

## Abstract

Radiation hydrodynamics studies the interaction between the motion of a high-temperature hyper-sonic plasma and the radiation it emits or absorbs. This coupling is central to many astrophysical phenomena, particularly those linked to the processes of accretion and ejection. The HADES code has been developed to model such systems. It provides a realistic model of radiative transport, which is particularly crucial for optically intermediate media. The code HADES couples the hydrodynamic equations with M1-gray or M1-multigroup radiative transfer models.

However, radiative hydrodynamics simulations remain extremely costly in terms of computing time, due to two main limitations. Firstly, the M1-multigroup model is based on a closure relation that has no analytic form, requiring costly numerical estimates. Secondly, the Courant-Friedrichs-Lewy condition strongly constrains the explicit schemes used in HADES, imposing very small time steps. In order to overcome these limitations, two complementary strategies based on Artificial Intelligence were developed during the course of this thesis.

The first involves training a neural network of the Multi-Layer Perceptron type to approximate the closure relation of the M1-multigroup model. This innovative approach achieves excellent accuracy while reducing the computational cost by a factor 3 000, making it to date the most efficient method known for this type of calculation. This substantial increase in performance has enabled us to carry out simulations of radiative shocks - shock waves in which the radiation directly influences the structure and dynamics of the front - with an unprecedented level of accuracy. In particular, we have been able to quantify the influence of a detailed spectral description of the radiation: the more detailed the description, the more the shock is slowed down, and the larger the size of the radiative precursor.

The second approach is based on the use of Physics-Informed Neural Networks to solve the radiative hydrodynamics equations directly and extrapolate the simulations beyond their initial time domain. These networks, which are also based on Multi-Layer Perceptrons, explicitly integrate the physical equations into their cost function, thereby guiding learning towards solutions that are consistent with the fundamental laws of physics. Initial tests, carried out on purely hydrodynamic shock configurations, show that this method can deal effectively with discontinuities and faithfully reconstruct the evolution of the shock at later times. On the other hand, the application to radiative shocks is proving more difficult: the networks are struggling to extrapolate the profiles correctly, due to the high residuals associated with the radiative hydrodynamics equations, which indicates that the physical constraints are not being properly respected. This last phase still requires in-depth investigation in order to understand the origin of these difficulties and to propose appropriate improvement strategies.

**Keywords:** Radiative hydrodynamics, Artificial Intelligence, Radiative shocks, Numerical simulation
# Apprentissage machine pour l'hydrodynamique radiative en astrophysique

## Résumé

L'hydrodynamique radiative étudie l'interaction entre le mouvement d'un plasma hypersonique souvent à haute température et le rayonnement qu'il émet ou absorbe. Ce couplage est central dans de nombreux phénomènes astrophysiques, notamment ceux liés aux processus d'accrétion et d'éjection. Pour modéliser de tels systèmes, le code HADES a été développé. Il offre une modélisation réaliste du transport radiatif, particulièrement cruciale pour les milieux optiquement intermédiaires. Le code HADES couple les équations de l'hydrodynamique avec des modèles de transfert radiatif de type M1-gris ou M1-multigroupe.

Cependant, les simulations en hydrodynamique radiative demeurent extrêmement coûteuses en temps de calcul, en raison de deux limitations principales. D'une part, le modèle M1-multigroupe repose sur une relation de fermeture qui n'admet pas de forme analytique, ce qui impose des estimations numériques coûteuses. D'autre part, la condition de Courant-Friedrichs-Lewy contraint fortement les schémas explicites utilisés dans HADES, imposant des pas de temps très réduits. Afin de surmonter ces limitations, deux stratégies complémentaires basées sur l'Intelligence Artificielle ont été développées au cours de cette thèse.

La première consiste à entraîner un réseau de neurones de type Multi-Layer Perceptron pour approcher la relation de fermeture du modèle M1-multigroupe. Cette approche novatrice atteint une excellente précision tout en réduisant le coût de calcul d'un facteur 3 000, faisant d'elle à ce jour la méthode la plus efficace connue pour ce type de calcul. Ce gain substantiel de performance nous a permis de réaliser des simulations de chocs radiatifs, des ondes de choc où le rayonnement influence directement la structure et la dynamique du front, avec un niveau de précision inédit. En particulier, nous avons pu quantifier l'influence d'une description spectrale fine du rayonnement : plus cette description est détaillée, plus le choc est ralenti, et plus la taille du précurseur radiatif augmente.

La seconde approche repose sur l'utilisation des Physics-Informed Neural Networks pour résoudre directement les équations de l'hydrodynamique radiative et extrapoler les simulations au-delà de leur domaine temporel initial. Ces réseaux, également basés sur des Multi-Layer Perceptrons, intègrent explicitement les équations physiques dans leur fonction de coût, orientant ainsi l'apprentissage vers des solutions cohérentes avec les lois fondamentales de la physique. Les premiers tests, menés sur des configurations de chocs purement hydrodynamiques, montrent que cette méthode permet de traiter efficacement les discontinuités et de reconstituer fidèlement l'évolution du choc à des temps ultérieurs. En revanche, l'application aux chocs radiatifs s'avère plus délicate : les réseaux peinent à extrapoler correctement les profils, en raison de résidus élevés associés aux équations de l'hydrodynamique radiative, ce qui indique un mauvais respect des contraintes physiques. Cette dernière phase nécessite encore des investigations approfondies afin de comprendre l'origine de ces difficultés et de proposer des stratégies d'amélioration adaptées.

**Mots clés** : Hydrodynamique radiative, Intelligence Artificielle, Choc radiatifs, Simulation numérique
ACKNOWLEDGEMENTS

Je souhaite tout d'abord exprimer ma profonde gratitude à Claire Michaut, pour m'avoir proposé ce sujet de thèse et pour m'avoir initié à l'univers de l'hydrodynamique radiative. Merci pour ton encadrement, ta disponibilité et tous tes précieux conseils tout au long de ces années. J'adresse également mes remerciements à Andrew Comport, qui, malgré des difficultés personnelles, m'a grandement aidé à développer les méthodes présentées dans ce manuscrit, en partageant avec moi de nombreuses idées et en m'orientant vers des pistes existantes à explorer.

Je remercie chaleureusement mes rapporteurs, Jérôme Breil et Emmanuel Franck, pour avoir accepté d'évaluer ce travail et pris le temps d'en rédiger un rapport. J'adresse également mes remerciements à l'ensemble du jury, Stéphane Mazevet, Florence Hubert et Serge Bouquet, pour leur présence et leur participation à ma soutenance.

Je suis reconnaissant aux membres de mon comité de suivi de thèse, André Ferrari et Frédéric Paletou, pour le temps qu'ils ont consacré chaque année à suivre l'avancement de mes travaux et pour leur écoute bienveillante. Je souhaite remercier tout particulièrement Frédéric Paletou : sans lui, je n'aurais probablement jamais entrepris une thèse en hydrodynamique radiative. Ses cours de M2 sur le transfert radiatif et ses nombreuses digressions sur leurs applications m'ont profondément inspiré et donné l'envie d'approfondir ce domaine. Mes remerciements vont également à Malik Chami, avec qui j'ai pu avoir des échanges enrichissants sur l'avenir de l'intelligence artificielle dans nos disciplines, et qui m'a plusieurs fois débloqué sur des questions techniques, malgré la distance entre nos champs de recherche.

Je tiens à remercier le professeur Vladimir Mikhailovitch Lipunov, qui m'a donné mes tout premiers cours de transfert radiatif à l'Université de Moscou. Malgré l'exigence de son enseignement, ces cours ont suscité en moi une véritable passion pour la discipline. Je remercie également le professeur Alexeï Nikolaevitch Rubtsov, qui m'a fait découvrir les méthodes d'intelligence artificielle appliquées à l'astrophysique et éveillé un vif intérêt pour ce domaine en plein essor. Без вашего преподавания эта диссертация была бы невозможна. Я хотел бы выразить вам самую искреннюю благодарность.

Je souhaite exprimer ma reconnaissance à Andrea Chiavassa, pour son aide précieuse dans la recherche d'opportunités après la thèse, et pour nos nombreuses discussions sur l'hydrodynamique radiative et les codes de simulation. J'espère vivement que nos chemins se recroiseront à l'avenir et que nous aurons l'occasion de collaborer.

Je remercie l'Institut 3IA Côte d'Azur et l'Université Côte d'Azur pour avoir financé cette thèse et rendu possible tout le travail présenté dans ce manuscrit.
Mes remerciements s'adressent également au frère Manuel-Marie, pour ses conseils éclairés sur le sens du travail de recherche et la gestion du temps. Je tiens à exprimer ma gratitude à Sophie et Loïc, qui animent avec bienveillance l'aumônerie de Nice et ont été des hôtes exceptionnels au foyer Saint Jean-Paul II et bienheureux Christian Chessel.

Je souhaite remercier de tout cœur ma famille, et plus particulièrement ma mère, pour son soutien indéfectible tout au long de ces années. J'adresse également un chaleureux merci à Alessandra Marelli, pour sa présence et son soutien constant.

Je suis profondément reconnaissant envers Khaled, qui régale chaque jour les membres de l'Observatoire avec ses plats délicieux à la cantine, et qui crée des moments de détente inoubliables grâce à ses blagues et à la bonne humeur qu'il nous transmet quotidiennement. Sans ton enthousiasme et ton talent, la vie à l'Observatoire de la Côte d'Azur ne serait pas la même !

Je remercie également Coretin Pecontal, qui m'a grandement aidé lors des tests de PINNs durant les derniers mois de ma thèse, alors que je rédigeais ce manuscrit. J'ai réellement apprécié ton indépendance, ta curiosité et ton enthousiasme : j'espère que nous aurons l'occasion de collaborer de nouveau !

Enfin, j'adresse mes remerciements les plus chaleureux à tous les collègues et amis ren-
contrés à l'Observatoire : Paul, Louise, Manon, Elisa, Chloé, Gabriel, Marie, Katherine, Kate,
Adrien, Fabiola et Carlo. Votre présence et votre camaraderie ont rythmé et enrichi ces trois
années de thèse. Un grand merci également à tous mes amis niçois pour les soirées dansantes,
les balades et tous les souvenirs partagés. Je pense en particulier à Carlotta, Campe, Mariam,
Ségolène, Benoît, Guillaume, Alban, Claire, Timothée et Pierre : merci pour tous ces mo-
ments de joie. Je souhaite enfin remercier chaleureusement Anna Roussel, dont le soutien et
les nombreuses idées m'ont accompagné une grande partie de cette thèse.
"Ce que vous faites dans la vie résonne dans l'éternité."

Film

"Oui, sans cesse un monde se noie
Dans les feux d'un nouveau soleil,
Les cieux sont toujours dans la joie;
Toujours un astre a son réveil,
Partout où s'abaisse ta vue,
Un soleil levant te salue,
Les cieux sont un hymne sans fin! "

Hymne du matin, Harmonies poétiques et religieuses, Lamartine
TABLE OF CONTENT

<table><tr><td>Abstract</td><td>iv</td></tr><tr><td>Acknowledgements</td><td>vi</td></tr><tr><td>Table of Content</td><td>ix</td></tr><tr><td>List of Figures</td><td>xii</td></tr><tr><td>List of Tables</td><td>xviii</td></tr><tr><td>List of Acronyms</td><td>xx</td></tr><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 Radiative hydrodynamics</td><td>11</td></tr><tr><td>2.1 Hydrodynamics</td><td>12</td></tr><tr><td>2.1.1 Mass conservation</td><td>13</td></tr><tr><td>2.1.2 Momentum conservation</td><td>14</td></tr><tr><td>2.1.3 Total energy conservation</td><td>15</td></tr><tr><td>2.1.4 Euler's equations</td><td>16</td></tr><tr><td>2.2 Radiative transfer</td><td>21</td></tr><tr><td>2.2.1 Radiative transfer equation</td><td>22</td></tr><tr><td>2.2.2 The moment equations</td><td>26</td></tr><tr><td>2.2.3 The Eddington factor</td><td>32</td></tr><tr><td>2.3 Coupling</td><td>42</td></tr><tr><td>2.3.1 Lorentz transformation</td><td>42</td></tr><tr><td>2.3.2 Radiative hydrodynamics equations</td><td>45</td></tr><tr><td>2.3.3 Asymptotic cases</td><td>46</td></tr><tr><td>2.4 Synthesis</td><td>50</td></tr></table>
<table><tbody><tr><td><b>3</b></td><td><b>Calculation of the Eddington factor with Artificial Intelligence</b></td><td><b>51</b></td></tr><tr><td>3.1</td><td>Steps to generate a functional neural network</td><td>53</td></tr><tr><td>3.1.1</td><td>Multi-Layer Perceptrons (MLP)</td><td>55</td></tr><tr><td>3.1.2</td><td>Training</td><td>58</td></tr><tr><td>3.1.3</td><td>Underfitting and overfitting</td><td>62</td></tr><tr><td>3.2</td><td>Eddington factor in the M1-multigroup model</td><td>64</td></tr><tr><td>3.2.1</td><td>“Historical” methods</td><td>66</td></tr><tr><td>3.2.2</td><td>Description of the new method based on Artificial Intelligence</td><td>68</td></tr><tr><td>3.2.3</td><td>Performance of the proposed method</td><td>85</td></tr><tr><td>3.3</td><td>Synthesis</td><td>89</td></tr><tr><td><b>4</b></td><td><b>Study of radiative shocks</b></td><td><b>91</b></td></tr><tr><td>4.1</td><td>Introduction to radiative shocks</td><td>92</td></tr><tr><td>4.2</td><td>Influence of the spectral nature of light</td><td>97</td></tr><tr><td>4.2.1</td><td>Simulation description</td><td>98</td></tr><tr><td>4.2.2</td><td>Results</td><td>100</td></tr><tr><td>4.3</td><td>Synthesis</td><td>108</td></tr><tr><td><b>5</b></td><td><b>Extrapolation of radiative shock simulations with Physics-Informed Neural Networks</b></td><td><b>109</b></td></tr><tr><td>5.1</td><td>Introduction to PINNs</td><td>110</td></tr><tr><td>5.1.1</td><td>Cost function</td><td>112</td></tr><tr><td>5.1.2</td><td>Selecting data</td><td>116</td></tr><tr><td>5.2</td><td>Extrapolation of radiative shock simulations with Physics-Informed Neural Networks</td><td>118</td></tr><tr><td>5.2.1</td><td>Hydrodynamic tests</td><td>120</td></tr><tr><td>5.2.2</td><td>Radiative hydrodynamic tests</td><td>128</td></tr><tr><td>5.3</td><td>Synthesis</td><td>139</td></tr><tr><td><b>6</b></td><td><b>Conclusions and perspectives</b></td><td><b>141</b></td></tr><tr><td colspan="2"><b>Conclusions and perspectives</b></td><td><b>141</b></td></tr><tr><td><b>A</b></td><td><b>Basic concepts on the Planck function</b></td><td><b>145</b></td></tr><tr><td>A.1</td><td>Derivatives of the Planck function</td><td>145</td></tr><tr><td>A.2</td><td>Antiderivatives of the Planck function</td><td>146</td></tr><tr><td>A.3</td><td>Functions derived from the Planck function</td><td>148</td></tr></tbody></table>
<table><tr><td>B</td><td>Numerical calculation of the radiative quantities</td><td>157</td></tr><tr><td></td><td>B.1 “General” case</td><td>161</td></tr><tr><td></td><td>B.2 Case close to isotropy</td><td>162</td></tr><tr><td>C</td><td>Computing the Eddington factor using search algorithms</td><td>165</td></tr><tr><td></td><td>C.1 Search algorithms</td><td>166</td></tr><tr><td></td><td>C.2 Comparison of the search algorithms</td><td>169</td></tr><tr><td>D</td><td>Hydrodynamic shock on a wall</td><td>171</td></tr><tr><td></td><td>Bibliography</td><td>175</td></tr></table>
LIST OF FIGURES

<table><thead><tr><th>FIGURE</th><th></th><th>Page</th></tr></thead><tbody><tr><td>1.1</td><td>Examples of applications in the field of radiative hydrodynamics.</td><td>3</td></tr><tr><td>2.1</td><td>Schematic representation of a radiation beam crossing the surface dS and propagating within the solid angle dΩ.</td><td>22</td></tr><tr><td>2.2</td><td>Types of interactions between light and matter.</td><td>23</td></tr><tr><td>2.3</td><td>Electromagnetic spectrum. Example of discretization into G groups: the 1st group includes gamma rays, the 2nd includes X-rays, and the last includes long-wavelength radiation.</td><td>29</td></tr><tr><td>2.4</td><td>Evolution of the Eddington factor χ<sub>g</sub> at different radiative temperatures T<sub>g</sub>, using a group narrowness δ<sub>g</sub> = 10<sup>-4</sup>.</td><td>36</td></tr><tr><td>2.5</td><td>Ratio between the Eddington factor of the M1-multigroup model, χ<sub>g</sub>, and that of the M1-gray model, χ<sub>R</sub>, for a reduced flux f<sub>g</sub> = 0.65, at different values of the group finesse δ<sub>g</sub>.</td><td>38</td></tr><tr><td>2.6</td><td>Evolution of the ratio between the Eddington factor of the M1-multigroup model, χ<sub>g</sub>, and that of the M1-gray model, χ<sub>R</sub>, as a function of δ<sub>g</sub>, illustrating the shift of the curve as δ<sub>g</sub> approaches 1.</td><td>38</td></tr><tr><td>2.7</td><td>Ratio between the Eddington factor of the M1-multigroup model, χ<sub>g</sub>, and that of the M1-gray model, χ<sub>R</sub>, for f<sub>g</sub> = 0.65 and v<sub>g-1/2</sub> = 0. I<sub>b</sub>: domain in which χ<sub>R</sub> is valid; I<sub>c</sub>: intermediate-to-high radiative temperatures; H: asymptotically high radiative temperatures.</td><td>40</td></tr><tr><td>2.8</td><td>Ratio between the Eddington factor of the M1-multigroup model, χ<sub>g</sub>, and that of the M1-gray model, χ<sub>R</sub>, for f<sub>g</sub> = 0.65 and v<sub>g+1/2</sub> → ∞. L: asymptotically low radiative temperatures; I<sub>a</sub>: low-to-intermediate radiative temperatures; I<sub>b</sub>: domain in which χ<sub>R</sub> is valid.</td><td>41</td></tr><tr><td>3.1</td><td>Diagram illustrating the reasoning process involved in developing an artificial intelligence strategy.</td><td>55</td></tr></tbody></table>
<table>
  <tbody>
    <tr>
      <td>3.2</td>
      <td>Diagram of a perceptron</td>
      <td>56</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Examples of activation functions that can be used in neural networks.</td>
      <td>57</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>Example of an MLP with five inputs and three outputs. The network consists of an input layer, two hidden layers, and an output layer. The red circles represent the neurons (or perceptrons), and σ denotes the activation function applied to each of them. The same activation function may be used for all neurons, or different functions may be employed depending on the specific needs of the model.</td>
      <td>58</td>
    </tr>
    <tr>
      <td>3.5</td>
      <td>Different possible outcomes after the training of a neural network.</td>
      <td>63</td>
    </tr>
    <tr>
      <td>3.6</td>
      <td>Schematic overview of the developed method, for the case where v<sub>g-1/2</sub> &gt; 0 and v<sub>g+1/2</sub> &lt; +∞. χ<sub>g,L</sub> and χ<sub>g,H</sub> correspond respectively to the polynomials defined by equations (3.6) and (3.5), used to compute χ<sub>g</sub> in domains L and H, while χ<sub>R</sub> denotes the Eddington factor of the M1-gray model (see equation (2.71)), used in domain I<sub>b</sub>. The MLPs are used to approximate the closure in the intermediate domains I<sub>a</sub>, I<sub>c</sub>, I<sub>d</sub>, I<sub>e</sub>, and I<sub>f</sub>. The boundaries of domains L, H, and I<sub>a</sub> to I<sub>f</sub> are defined in table 3.1. The brown lines delineate the domains of applicability of the different approaches.</td>
      <td>69</td>
    </tr>
    <tr>
      <td>3.7</td>
      <td>Ratio between the Eddington factor of the M1-multigroup model, χ<sub>g</sub>, and that of the M1-gray model, χ<sub>R</sub>, for a reduced flux f<sub>g</sub> = 0.65, in the case where v<sub>g-1/2</sub> &gt; 0 and v<sub>g+1/2</sub> &lt; +∞, for different values of the group narrowness δ<sub>g</sub>.</td>
      <td>70</td>
    </tr>
    <tr>
      <td>3.8</td>
      <td>Ratio between the Eddington factor of the M1-multigroup model, χ<sub>g</sub>, and that of the M1-gray model, χ<sub>R</sub>, for a reduced flux f<sub>g</sub> = 0.65, in the cases where v<sub>g-1/2</sub> = 0 and where v<sub>g+1/2</sub> = +∞.</td>
      <td>71</td>
    </tr>
    <tr>
      <td>3.9</td>
      <td>Comparison of the Eddington factor computed using the polynomial expression χ<sub>g,H</sub>, defined by equation (3.5), with the reference data obtained by solving equations (3.3) and (3.4) using a bisection algorithm.</td>
      <td>75</td>
    </tr>
    <tr>
      <td>3.10</td>
      <td>Comparison of the computation of the Eddington factor between the polynomial expression χ<sub>g,L</sub>, defined by equation (3.6), and reference data obtained by solving equations (2.72), (2.73), and (2.74) using a line-search algorithm, for T<sub>g</sub> = 10<sup>-8</sup> and δ<sub>g</sub> = 10<sup>-4</sup>. This comparison also includes results obtained by solving the equations proposed by Minerbo (1978) using a Newton algorithm.</td>
      <td>77</td>
    </tr>
  </tbody>
</table>
3.11 MLP architectures used. lin, tanh, and $\sigma$ denote respectively the identity, hyperbolic tangent, and sigmoid activation functions. The “Normalization” and “Transformation” layers respectively apply equations (3.12) to (3.15), and (3.16). $\mathcal{T}$ is the dimensionless radiative temperature, defined differently depending on the MLP considered: $\mathcal{T}_g^{\nu_{g-1/2}}$ for MLP$_a$, $\mathcal{T}_g^{\nu_{g+1/2}}$ for MLP$_c$, MLP$_{e,1}$, and MLP$_{e,2}$, and $\mathcal{T}_g^s$ for MLP$_f$. 81

3.12 Learning curves of the different architectures tested for MLP$_a$. $L$ denotes the number of hidden layers and $N$ the total number of perceptrons. For example, the architecture L1N6 consists of 1 hidden layer with 6 perceptrons. 82

3.13 Comparison of the distribution of the relative error of the MLP$_a$ prediction for the training set and the validation set. 84

3.14 Initial and boundary conditions for the Marshak wave simulation. 86

3.15 Component $\mathbb{P}_{R,xx}$ of the total radiative pressure at time $t = 0.133$ ns. The relative error is given by the following formula: $(\mathbb{P}_{xx,method} - \mathbb{P}_{xx,line search})/\mathbb{P}_{xx,line search} \times 100$. 88

4.1 Different types of radiative shock structures. The red line represents the temperature, and the blue line represents the density. The position $x = 0$ corresponds to the location of the shock. 95

4.2 Types of radiative precursor structures in the case of optically thick shocks. The position $x = 0$ corresponds to the location of the shock. 97

4.3 Schematic illustration of the simulation of shock formation at the interface with a rigid wall. The simulation starts with homogeneous hydrodynamic conditions: density $\rho_0$, temperature $T_0$, and velocity $v_0$. At later times $t_1 < t_2 < t_3$, a shock front forms and gradually propagates away from the wall. 98

4.4 Evolution of the shock position as a function of simulated time. 101

4.5 Profiles of the physical quantities for the “Mach 4” simulation. The position $x = 0$ cm corresponds to the location of the shock. All these data are expressed in the shock frame. 102

4.6 Profiles of the physical quantities for the “Mach 8” simulation. The position $x = 0$ cm corresponds to the location of the shock. All these data are expressed in the shock frame. 103

4.7 Evolution of the shock velocity as a function of the group narrowness $\delta_g$. The case with $\delta_g = 0$ corresponds to the simulation using the M1-gray model. 104
LIST OF TABLES

<table><thead><tr><th>TABLE</th><th></th><th>Page</th></tr></thead><tbody><tr><td>3.1</td><td>Limits for each domain and the method used to compute the Eddington factor.</td><td>73</td></tr><tr><td>3.2</td><td>Description of the datasets used for training and validation of the neural networks.</td><td>78</td></tr><tr><td>3.3</td><td>Architecture of the neural networks used.</td><td>83</td></tr><tr><td>3.4</td><td>Results of the 10-fold cross-validation for MLP<sub>a</sub>.</td><td>84</td></tr><tr><td>3.5</td><td>Definition of the groups</td><td>85</td></tr><tr><td>3.6</td><td>Comparison of simulation times (CPU time).</td><td>87</td></tr><tr><td>4.1</td><td>Initial hydrodynamic quantities.</td><td>99</td></tr><tr><td>4.2</td><td>Frequency bounds of each group used in the M1-multigroup simulations.</td><td>100</td></tr><tr><td>4.3</td><td>Computational cost of the simulations in CPU seconds, hours, and days.</td><td>107</td></tr><tr><td>5.1</td><td>Initial conditions of the test case used to evaluate the PINNs.</td><td>120</td></tr></tbody></table>
**Adam** Adaptative Moment Estimation

**AI** Artificial Intelligence

**BFGS** Broyden-Fletcher-Goldfarb-Shano

**BNN** Bayesian Neural Networks

**CFL** Courant-Friedrichs-Lewy

**CMB** Cosmological Background

**CNN** Convolutional Neural Networks

**FLD** Flux Limited Diffusion

**GAN** Generative Adversarial Networks

**HADES** Hydrodynamique Adaptée à la Description des Ecoulements Supersoniques

**HLL** Harten-Lax-van Leer

**HLLE** Harten-Lax-van Leer-Einfeld

**ICF** Inertial Confinement Fusion

**L-BFGS** Limited-memory Broyden-Fletcher-Goldfarb-Shano

**LSTM** Long Short-Term Memory

**LTE** Local Thermodynamic Equilibrium

**MLP** Multi-Layer Perceptron
**NeuHPC** Neural High-Performance Computing

**ODE** Ordinary Differential Equation

**PDE** Partial Differential Equations

**PINN** Physics-Informed Neural Networks

**RNN** Recurrent Neural Networks

**SGD** Stochastic Gradient Descent

**SKA** Square Kilometer Array

**VAE** Variational Autoencoders

**VTEF** Variable Tensor Eddington Factor
<figure><img src="image_5.png" /></figure>
CHAPTER 1

INTRODUCTION

Astrophysics constitutes a true natural laboratory in which physical theories can be
tested under extreme conditions, often inaccessible on Earth. It thus enables us to
refine our models, verify the validity of our understanding of the fundamental laws
of physics in such regimes, and deepen our knowledge of how the Universe operates. This
field has, for instance, revealed that time is not an absolute quantity but a relative one: it
depends in particular on the velocity of the observer and on the gravitational field in which
they evolve. By its very nature, astrophysics is interdisciplinary, drawing on all branches of
physics, since astrophysical phenomena involve processes arising from quantum mechanics,
electromagnetism, thermodynamics, nuclear physics, and relativity. However, observations
of astrophysical objects alone are insufficient to capture the full complexity of the physical
mechanisms at play. A rigorous theoretical modelling is required to interpret these observations
and formulate general laws. This constant dialogue between theory and observation is a key
driver of scientific progress in the field. Thus, the increased precision of Mercury's orbital
measurements revealed an anomaly inexplicable by Newtonian mechanics, which led Einstein
to develop the theory of general relativity [1]. Among other things, this theory predicts
the existence of gravitational waves, minute distortions of space-time, whose existence was
experimentally confirmed in 2016 [2].

In particular, the extreme densities and temperatures reached during violent phenomena
such as supernovae or accretion around young stars produce intense radiation that interacts
with the surrounding matter and alters the structure of the observed flows. These complex
interactions require precise modelling, falling within the domain of radiative hydrodynamics.
Understanding such phenomena demands a combination of complementary approaches. Astrophysical observation plays a crucial role because, as mentioned earlier, it gives access to phenomena that cannot be reproduced on Earth. When feasible, laboratory reproduction of these processes also constitutes a fundamental tool, offering a controlled environment in which to isolate the mechanisms at work and verify that the physical laws invoked indeed account for what is observed in the cosmos. To establish a meaningful correspondence between terrestrial experiments and astrophysical situations, it is necessary to define appropriate scaling laws [3, 4, 5, 6]. Finally, and this is the central focus of this thesis, the development of theoretical models constitutes an indispensable pillar for addressing the full richness of the physical processes involved. Such modelling may be analytical, allowing one to isolate simple and well-understood mechanisms, even if this means neglecting nonlinear effects or complex interactions. It may also take the form of numerical simulations, which preserve a more complete formulation of the physical equations and allow one to explore more realistic configurations while maintaining a framework that remains interpretable from a theoretical standpoint.

## Scientific context

In this thesis, we focus on the modelling of violent phenomena in which radiative processes play a predominant role in the observed dynamics. But under what physical conditions do such radiative effects become significant? There is no sharp boundary separating a regime dominated by radiation from one in which radiation can be neglected. However, it is possible to introduce two dimensionless numbers that allow to assess the relative importance of radiative effects in a given configuration. The first is the Mihalas number R, which measures the ratio between the internal energy u of a gas and the radiative energy E_R. In the case of a monoatomic gas in radiative equilibrium<fnref n="1" />, this ratio can be expressed as:

$$R = \frac{u}{E_R} \approx \frac{1}{36} \frac{N}{T^3}.$$

This expression shows that the higher the temperature T of the medium, the more dominant the radiative energy becomes in the energy budget. Conversely, the larger the gas density N, the less significant the radiative effects are. The gas is said to be radiation-dominated when R ≪ 1. The second dimensionless number is the Boltzmann number Bo, which compares the energy flux transported by advection F with the energy flux transported by radiation F_R. It is given by:

$$Bo = \frac{F}{F_R} \approx \frac{1}{9} \frac{v}{c} \frac{N}{T^3} = 4 \frac{v}{c} R.$$

<footnote n="1">¹See definition in Chapter 2. *Radiative hydrodynamics*.</footnote>
<figure><img src="image_6.png" /><figcaption>(a) Cassiopeia A supernova remnant.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) Descriptive diagram of Inertial Confinement
Fusion.</figcaption></figure>

<caption>Figure 1.1: Examples of applications in the field of radiative hydrodynamics.</caption>

In the applications considered here, the fluids are non-relativistic, meaning that their velocity *v* is much smaller than the speed of light *c*. As a result, radiative transport mechanisms become important well before radiative energy becomes dominant in the overall energy balance. To sum up, the phenomena studied in this thesis correspond to hot and sufficiently low-density media in which the transport of radiative flux plays a crucial role, often even before the regime in which radiation dominates energetically is reached.

Such physical conditions arise naturally in astrophysics, a field in which radiative hydro-
dynamics plays a key role in understanding the dynamical and energetic phenomena affecting
stars and other stellar objects. In stars, for example, nuclear fusion produces intense radiation
that continually interacts with the surrounding plasma. This interaction profoundly influences
the internal structure, stability, and evolution of stars. Radiative transfer modelling is espe-
cially crucial for studying mass loss and the chemical composition of red supergiants, whose
extended outer envelopes are subject to complex radiative mechanisms [7, 8]. Supernova
remnants, formed from the cataclysmic explosions of dying stars, are also strongly shaped by
the coupled effects of radiation and hydrodynamics. The interaction between shock waves
and radiation, both during the initial explosion and throughout the subsequent evolution of the
ejecta, modulates their dynamics, can induce instabilities, and largely determines their observ-
able emission [9, 10]. Many questions remain open regarding the observed structures, such as
the filaments present in some supernova remnants (see figure 1.1a); more accurate modelling
that includes radiative effects could provide new insights into these unresolved issues.

The relevance of these phenomena is not purely astrophysical: they also have major
implications in the field of energy production, in particular for the development of clean and sustainable sources. A prominent example is nuclear fusion, envisioned as an energy alternative that produces little waste and emits no carbon dioxide. The most advanced method to date for reaching fusion conditions is the Inertial Confinement Fusion (ICF), which achieved a historic milestone in 2022 with the first net-energy ignition at the National Ignition Facility, where the energy produced by fusion exceeded the energy delivered to the system [11, 12]. In this approach, a capsule containing a mixture of deuterium and tritium is placed in a vacuum chamber and compressed and heated to extremely high temperatures by a radiative shock generated by X-ray radiation produced from powerful laser beams interacting with aohlraum target (see figure 1.1b). This process is extremely sensitive to hydrodynamic instabilities, particularly Rayleigh-Taylor type instabilities, similar to those encountered in supernovae, which can compromise the efficiency of the implosion by causing substantial energy losses [13, 14]. When poorly controlled, these instabilities can prevent the conditions required for fusion from being reached [15]. A detailed understanding of radiation-matter interactions, both in astrophysical settings and in laboratory environments, therefore constitutes a major scientific and technological challenge.

Existing numerical codes

Once the physical situations requiring radiative hydrodynamics modelling have been identified, it becomes necessary to determine how to describe these phenomena realistically while maintaining a reasonable computational cost. This involves choosing appropriate approximations that account for the characteristics of the media under consideration, as well as numerical limitations. Two classical simplifications are often used, especially in astrophysics: the diffusion approximation, suitable for optically thick media, and the free-streaming approximation, valid for optically thin ones. However, these models are not valid in intermediate regimes, where more sophisticated approaches are required to faithfully capture the physics of radiative transfer.

Among these advanced approaches, the first category comprises statistical Monte Carlo methods. These consist of simulating radiative transport through the stochastic emission of photon packets, whose direction and frequency are randomly sampled. The opacity of the medium then determines the probability of interactions (absorption or scattering) between the photons and the matter. These methods have the advantage of being able to describe all opacity regimes, without restrictive assumptions about the radiation-matter coupling. However, they are extremely costly in optically thick media, where photons undergo a lot of interactions before
being absorbed. Examples of codes using this approach include MC3D [16], RADMC [17], CRASH [18], and AREPO-MCRT [19].

A second class of methods relies on the spatial–angular discretisation of the specific intensity<fnref n="2" />, i.e. the fundamental quantity describing radiation. These methods, referred to as angular discretisation schemes, approximate integrals over the solid angle by sums over a finite set of discrete directions. They offer good accuracy at a moderate numerical cost. However, they are subject to two well-known artefacts: ray effects, caused by angular discretisation, may generate artificial structures in free-streaming regions; and numerical diffusion (*false scattering*), induced by spatial discretisation and interpolation, may produce an unphysical smoothing of the solution [20]. These artefacts can be mitigated through finer spatial and angular meshes. A representative example of this approach is the Athena++ code [21, 22, 23].

To overcome some of these difficulties, the long- and short-characteristics methods have been developed. These techniques integrate the radiative transfer equation along the trajectories of light rays. The long-characteristics method performs this integration from each grid point to the domain boundaries, accounting for boundary conditions. Although highly accurate, it is costly because it requires to integrate over long distances across the mesh. Conversely, the short-characteristics method restricts the integration to the neighbouring cell, drastically reducing computational cost. However, the necessary interpolations introduce increased numerical diffusion. The FLASH code illustrates this approach [24].

Another family of methods relies on solving the moment equations of the specific intensity. These methods solve the first angular moments of the radiative transfer equation and then close the system by introducing a relation between higher- and lower-order moments. The simplest model in this class is the Flux Limited Diffusion (FLD) model, which relies on the two lowest moments (radiative energy and flux) and assumes a closure valid only in optically thick media. Within this framework, the radiative flux is always aligned with the gradient of radiative energy, an assumption that becomes incorrect as soon as the medium is even moderately transparent. Codes such as ZEUS2D-FLD [25], COSMOS [26], CO<sub>5</sub>BOLD [27, 28, 29], and a module of FLASH [30] are based on this model. To improve accuracy in intermediate opacity regimes, the Variable Tensor Eddington Factor (VTEF) formalism was developed. This method also relies on the first three moments (radiative energy, flux, and pressure) but determines the closure relation linking pressure to energy via the Eddington tensor, computed by solving for the specific intensity from the full radiative transfer equation. This approach provides a precise estimate of the closure relation, though it is not strictly local with respect to radiative

<footnote n="2">²Definition given in Chapter 2. *Radiative hydrodynamics*.</footnote>
quantities. It captures radiation–matter coupling more effectively, but is numerically costly and complex to implement. Examples of codes implementing this method include ZEUS2D [31], ZEUSMP [32], and TITAN [33].

Finally, the model studied in this thesis is the M1 model, which also relies on the first three moments of the specific intensity. It proposes a local closure derived from an analytical formalism, in which the radiation field is assumed to be symmetric about the direction of the radiative flux. This assumption is equivalent to considering local thermodynamic equilibrium for the radiation, which allows expressing the Eddington tensor as a purely local function of the radiative quantities. The M1 model constitutes a good compromise between physical accuracy and computational cost: it is more general than the FLD model, notably by providing a better description of intermediate transport regimes, while being far less expensive than Monte Carlo, explicit angular discretisation, or VTEF approaches. It is used in several reference codes, such as HERACLES [34, 35], AREPO-RT [36], and, most importantly for this thesis, the HADES code (for Hydrodynamique Adaptée à la Description des Ecoulements Supersoniques) [37, 38, 39].

The HADES code is a massively parallelised Fortran simulation tool based on finite-volume numerical schemes. It enables modelling radiative transfer using two approaches: the M1-gray model, used for global radiation–matter interaction studies, and the M1-multigroup model, which offers a more refined treatment of spectral effects. This code has been widely used to study various stellar physics phenomena, such as supernova remnants [40, 41, 42, 43], accretion shocks in binary systems [44], jets from young stars [45], and, more fundamentally, radiative shocks [46]. Although the M1 model makes radiative hydrodynamics simulations more accessible than more expensive methods (Monte Carlo, angular discretisation, VTEF), these simulations remain extremely demanding in computational resources, particularly when spectral effects are included through the M1-multigroup model. This is due to two main reasons:

• Computation of the Eddington factor<fnref n="3" /> in the M1-multigroup model

In this model, the Eddington factor has no explicit analytical expression. The code therefore relies on numerical search algorithms to estimate it accurately. Although reliable, these methods are particularly costly in computational time.

• The Courant–Friedrichs–Lewy (CFL) condition

To ensure numerical stability when solving the radiative hydrodynamics equations, the

<footnote n="3">³Details provided in Section 2.2.3. *The Eddington factor*.</footnote>
time step $\Delta t$ and spatial step $\Delta x$ must satisfy a CFL type condition: $c\Delta t/\Delta x \le 1$, where c is the speed of light. This constraint forces the use of very small time steps, significantly increasing the computational cost of simulations.

The objective of this thesis is therefore to explore the possibility of reducing the cost of such simulations through approaches derived from Artificial Intelligence (AI), and in particular through the use of neural networks.

## Artificial Intelligence

Neural networks are machine-learning algorithms whose architecture is inspired by the functioning of the human brain. Their origin dates back to 1943, when Walter Pitts and Warren McCulloch proposed a simplified mathematical model to represent thought processes and decision-making in human cognition [47]. This model associated elementary processing units, artificial neurons, with logical functions, constituting an initial attempt to formalise human reasoning. These pioneering works quickly attracted the attention of prominent researchers. Among them, Alan Turing drew inspiration from this approach to develop the famous “Turing test” [48], designed to assess machine intelligence. This test, which consists in determining whether a human interlocutor can distinguish a machine from a human being through a conversational exchange, marked the beginning of the exploration of the capabilities of AI.

Over the following decades, research in AI and machine learning gave rise to a wide variety of neural network architectures, each tailored to specific types of problems. Among the most emblematic are:

* **MLP (Multi-Layer Perceptron) networks:** this is the simplest and most widespread architecture. An MLP consists of several layers of interconnected artificial neurons arranged hierarchically. It can model complex nonlinear relationships between inputs and outputs and remains widely used for classification and regression tasks;

* **Generative Adversarial Networks (GAN, Generative Adversarial Networks):** introduced in the 2010s, GANs rely on a dual architecture composed of a generator and a discriminator trained in competition. The generator attempts to produce data realistic enough to fool the discriminator, while the latter tries to distinguish real data from synthetic data. This competitive dynamic leads the generator to progressively improve the realism of the generated samples, enabling innovative applications in image synthesis, data augmentation, and simulation;
* **Variational Autoencoders (VAE, Variational Autoencoders):** these probabilistic generative models learn a compact and continuous representation of the input data. They are particularly effective for anomaly detection, as they allow the identification of samples that deviate significantly from the learned distribution, thereby revealing atypical behaviours or errors.

These architectures illustrate the diversity and richness of neural network approaches. Each model offers specific advantages depending on the nature of the data and the intended application. For a more comprehensive presentation of these concepts, the reader is referred to standard references such as Bishop (1995) [49] and Wilamowski (2009) [50]. Today, AI techniques are experiencing rapid growth in astrophysics, where their applications can generally be grouped into three main categories:

1. **Instrumentation control and robotics:** in space missions, AI enhances the planning of observations onboard Martian rovers or observational satellites [51, 52, 53, 54]. In ground-based observatories, it eases the automation of optical-instrument calibration, fault diagnosis, and optimisation of image quality [55];

2. **Data processing:** AI plays a central role in analysing observational data. It enables, for instance, the extraction of the Cosmological Background (CMB) signal from astrophysical and instrumental foregrounds [56, 57, 58, 59, 60], the detection of gravitational waves in complex and non-stationary noise backgrounds [61], or the optimisation of massive data-flow management for the Square Kilometer Array (SKA) project by reducing storage needs through intelligent detection of objects and anomalies during acquisition [62, 63];

3. **Simulation and modelling:** more directly related to this thesis, AI is used to improve the accuracy and efficiency of numerical simulations. Applications include super-resolution of physical fields [64] or the generation of adaptive meshes, computed more rapidly than with traditional techniques [65].

A detailed review of AI applications in astrophysics is provided in the article by Smith (2023) [66]. Despite its growing adoption, AI remains a relatively young tool that must be carefully adapted and integrated into the specificities of each domain, particularly in the context of numerical simulations. In this regard, a new research field has recently emerged: Neural High-Performance Computing (NeuHPC) [67], which aims to integrate AI models directly into high-performance computing codes. The present thesis is fully aligned with this approach.
Among the neural network strategies applied to numerical simulation, a particularly notable class is the Physics-Informed Neural Networks (PINN)s. Introduced by Raissi et al. (2019) [68], building upon the pioneering work of Lagaris et al. (1998) [69], PINNs constitute an innovative approach for solving problems governed by PDEs. The core idea is to approximate the solution of such a problem using a neural network trained by minimising a loss function comprising several components: the PDE residual, the initial conditions, and the boundary conditions when available. The spatial and temporal derivatives required to evaluate the residual are computed through automatic differentiation [70], a technique based on decomposing any computer program into a sequence of elementary operations and standard functions. This class of networks is distinctive in that it embeds physical laws directly into the learning process, significantly reducing the amount of required data while ensuring that the predicted solutions remain consistent with the fundamental equations governing the studied problem.

## Contributions

This thesis is organised around three main scientific contributions, each corresponding to a clearly identified research direction:

1. The development of an efficient method for computing the closure relation in the M1-multigroup model, enabling the incorporation of the spectral character of radiation into radiative hydrodynamics simulations (see Chapters 2. *Radiative hydrodynamics* and 3. *Calculation of the Eddington factor with Artificial Intelligence*);

2. The analysis of the impact of a fine spectral modelling of radiation on the structure of radiative shocks (see Chapter 4. *Study of radiative shocks*);

3. A first exploration of the potential of Physics-Informed Neural Networks (PINN) to extrapolate radiative hydrodynamics simulation data to later instants (see Chapter 5. *Extrapolation of radiative shock simulations with Physics-Informed Neural Networks*).

The first two parts of this work have led to the publication of the following articles:

* Radureau, G., Michaut, C. & Comport, A. I., AI-based computation method for the Eddington factor in the M1-multigroup model, *Phys. Rev. E* **111**, 035301 (2025);

* Radureau, G. & Michaut, C., Impact of frequency-dependent radiation on the dynamics and structure of radiative shocks, *Phys. Rev. E* **111**, 045213 (2025).
Beyond its scientific contributions, this manuscript also aims to serve as a useful reference for those wishing to deepen or expand the use of AI in the context of radiative hydrodynamics simulations. With this in mind, I have sought to adopt an approach as pedagogical as possible, both in the presentation of the foundations of radiative hydrodynamics and in the introduction of the AI methods employed. The manuscript is structured around four main chapters, in addition to the present introduction, as well as four technical appendices. These various components introduce the key concepts used in this work, present the results obtained, and detail the methodological aspects necessary for the reproducibility of the proposed approaches.

Chapter 2 is devoted to the physical modelling of radiative hydrodynamics used in this study. The M1-multigroup model, as implemented in the HADES code, is described with an emphasis on the assumptions and simplifications adopted.

Chapter 3 introduces the main AI techniques applied to numerical simulation, and then presents an original method developed during this thesis to efficiently estimate the closure relation of the M1-multigroup model using neural networks.

Chapter 4 analyses the influence of a rigorous spectral modelling of radiation on the internal structure of radiative shocks, relying on the Eddington factor computation method developed in the preceding chapter.

Finally, Chapter 5 presents preliminary attempts to use Physics-Informed Neural Networkss for the extrapolation of radiative shock simulations, thereby opening the way to new perspectives for applying neural networks in radiative hydrodynamics.
RADIATIVE HYDRODYNAMICS

Astrophysical fluids, such as those found in stars or supernovae, interact strongly with the surrounding radiation. Their flow must therefore be described by the equations of fluid mechanics, such as the Euler or Navier-Stokes equations. However, to model these interactions accurately, it is also necessary to account for the propagation of radiation and its exchanges of energy and momentum with matter. The combined study of these phenomena falls within the field of radiative hydrodynamics. In this context, the modelling of matter-radiation interactions depends on the degree of equilibrium reached by the system. Three levels of equilibrium may be considered:

1. **Chemical equilibrium:** this is achieved when chemical kinetics have evolved toward a steady state, meaning that all chemical reactions are balanced and that the concentration of each species remains constant in time;

2. **Thermal equilibrium:** this is achieved when all species constituting the fluid share a single temperature, indicating that internal energy exchanges are sufficiently rapid to homogenise the thermal distribution;

3. **Radiative equilibrium:** this is achieved when matter is in equilibrium with radiation, implying that the radiation temperature is identical to that of the matter.

When a fluid is simultaneously in chemical and thermal equilibrium, it is said to be in Local Thermodynamic Equilibrium (LTE). This concept is fundamental, as it determines how
matter–radiation interactions must be modelled. Indeed, depending on whether LTE and radiative equilibrium are assumed, the complexity of the model can vary significantly.

In this chapter, we will assume that the matter is in LTE, while developing a theoretical framework that does not rely on the assumption of radiative equilibrium, in order to better capture the out-of-equilibrium effects frequently observed in many astrophysical environments. To clearly introduce the modelling assumptions adopted and to lay the groundwork for the approach used in the HADES code, I will begin by presenting the equations of classical hydrodynamics in the absence of interaction with radiation. I will then describe the radiative transport model in a fluid at rest, before detailing how interactions between radiation and a moving fluid are accounted for in the model implemented in HADES.

## 2.1 Hydrodynamics

Let us begin by establishing the general framework for describing a moving fluid, assuming that it does not interact with any external element. In such a context, the fluid evolves autonomously, without exchanging mass, energy, or momentum with its environment. It therefore satisfies the fundamental conservation principles:

1. **Mass conservation:** The fluid neither gains nor loses matter during its evolution. This principle is expressed by the continuity equation, which states that the total mass contained in a moving fluid element remains constant;

2. **Momentum conservation:** In the absence of external forces, the momentum of the fluid is preserved, in accordance with the fundamental principle of dynamics. This is described by the Euler equations (for an inviscid fluid) or by the Navier–Stokes equations (when viscosity is taken into account);

3. **Energy conservation:** The total energy of the fluid (including kinetic, internal, and potential energy) remains constant in the absence of any external interaction. This conservation is expressed by the energy equation, which links variations of internal energy to heat transfer and the work of forces.

These laws form the foundation of fluid dynamics and allow the behaviour of a fluid to be described within a purely hydrodynamic framework. We now examine the modelling of such fluids before introducing the influence of interactions with radiative phenomena.
### 2.1.1 Mass conservation

Let $V$ be a control volume in space and $m$ the total mass of the fluid contained within it. This mass can be expressed as follows:

$$m(t) = \iiint_V \rho(t, \vec{x}) \, dV, \qquad (2.1)$$

where $\rho(t, \vec{x})$ denotes the mass density of the fluid at a point $\vec{x}$ within the control volume $V$ at time $t$. Since the mass of the fluid is conserved over time, we can write:

$$\frac{dm(t)}{dt} = \frac{d}{dt} \iiint_V \rho(t, \vec{x}) \, dV = 0. \qquad (2.2)$$

However, according to the Reynolds transport theorem, the time derivative of the integral can be written as:

$$\frac{d}{dt} \iiint_V \rho \, dV = \iiint_V \partial_t \rho \, dV + \iint_S \rho \vec{v} \cdot \vec{n} \, dS, $$

with $\vec{v}$ the velocity field of the fluid, $S$ the boundary surface of the control volume $V$, and $\vec{n}$ the unit normal vector to the surface element $dS$. Using the Green–Ostrogradsky theorem, the surface integral can be rewritten as a volume integral. Thus, we obtain the expression:

$$\iint_S \rho \vec{v} \cdot \vec{n} \, dS = \iiint_V \vec{v} \cdot (\rho \vec{v}) \, dV .$$

Thus, the mass conservation equation (2.2) can be rewritten in the form:

$$\iiint_V (\partial_t \rho + \vec{v} \cdot (\rho \vec{v})) \, dV = 0 .$$

Considering that this equation is valid for any control volume $V$, we can deduce the continuity equation, also called the local conservation equation of mass:

$$\partial_t \rho + \vec{v} \cdot (\rho \vec{v}) = 0. \qquad (2.3)$$

It is important to emphasise that the transition from the integrated equation to the local equation is valid “almost everywhere”, in the sense of Lebesgue integrals. This means that at the location of a discontinuity, such as a shock, the local equation is no longer valid. This subtlety will be crucial when discussing of Physics-Informed Neural Networks (see Section 5.1. Introduction to PINNs).
## 2.1.2 Momentum conservation

To establish the conservation law for momentum, we rely on Newton's second law of motion. This principle states that the time derivative of the momentum of the particles, $\vec{p}_m$, is equal to the sum of the external forces $\vec{F}_{ext}$ acting on them.

$$ \frac{d\vec{p}_m}{dt} = \sum \vec{F}_{ext} . \qquad (2.4) $$

In the present case, the momentum of a fluid contained in a control volume V corresponds to the sum of the momenta of the particles that compose it. It can therefore be expressed as:

$$ \vec{p}_m = \iiint_V \rho \vec{v} \, dV . $$

For the same reasons as in Section 2.1.1. *Mass conservation*, using the Reynolds transport theorem as well as the Green–Ostrogradsky theorem, and writing $\vec{v} = (v_1, v_2, v_3)$, the time derivative of the momentum can be written as:

$$ \frac{d\vec{p}_m}{dt} = \frac{d}{dt} \iiint_V \rho \vec{v} \, dV = \left[ \iiint_V \left( \partial_t (\rho v_1) + \vec{v} \cdot (\rho v_1 \vec{v}) \right) \right. \\ \left. \qquad \iiint_V \left( \partial_t (\rho v_2) + \vec{v} \cdot (\rho v_2 \vec{v}) \right) \right. \\ \left. \qquad \iiint_V \left( \partial_t (\rho v_3) + \vec{v} \cdot (\rho v_3 \vec{v}) \right) \right] . $$

Considering the coordinates $\vec{x} = (x_1, x_2, x_3)$, the divergence term in the previous equation can be expanded as follows:

$$ \begin{aligned} \begin{bmatrix} \vec{\nabla} \cdot (\rho v_1 \vec{v}) \\ \vec{\nabla} \cdot (\rho v_2 \vec{v}) \\ \vec{\nabla} \cdot (\rho v_3 \vec{v}) \end{bmatrix} &= \begin{bmatrix} \partial_{x_1}(\rho v_1^2) + \partial_{x_2}(\rho v_1 v_2) + \partial_{x_3}(\rho v_1 v_3) \\ \partial_{x_1}(\rho v_2 v_1) + \partial_{x_2}(\rho v_2^2) + \partial_{x_3}(\rho v_2 v_3) \\ \partial_{x_1}(\rho v_3 v_1) + \partial_{x_2}(\rho v_3 v_2) + \partial_{x_3}(\rho v_3^2) \end{bmatrix} = \vec{\nabla} \cdot \begin{pmatrix} \rho \begin{bmatrix} v_1 \\ v_2 \\ v_3 \end{bmatrix} \begin{bmatrix} v_1 & v_2 & v_3 \end{bmatrix} \end{pmatrix} \\ &= \vec{\nabla} \cdot (\rho \vec{v} \otimes \vec{v}) , \end{aligned} $$

where $\otimes$ denotes the dyadic product. Thus, the time derivative of the momentum can be written more concisely as follows:

$$ \frac{d\vec{p}_m}{dt} = \iiint_V \left( \partial_t (\rho \vec{v}) + \vec{v} \cdot (\rho \vec{v} \otimes \vec{v}) \right) dV . $$

We still need to determine the external forces acting on the fluid. We distinguish between *surface forces* $\vec{F}_S$, which act on the boundary of the control volume V, and *body forces* $\vec{F}_V$, which act on every particle of the fluid.
For the surface forces, we may distinguish between *pressure forces* and *viscous stress forces*, which arise from the interaction of the particles at the boundary of the control volume *V* with their environment. In the present case, we consider only the *pressure forces* *p*, which take the following form:

$$ \vec{F}_S = \oint_S -p \vec{n} \, dS . $$

According to the Green-Ostrogradski theorem, this force can be rewritten in the form:

$$ \vec{F}_S = - \iiint_V \vec{\nabla} p \, dV . $$

Regarding the *body forces*, various contributions may be taken into account, such as gravitational forces, electromagnetic forces, and, in particular, the radiative force. For now, we will neglect body forces, and later detail in Section 2.3.2. *Radiative hydrodynamics equations* how to include the radiative force. Finally, the fundamental principle of dynamics can be rewritten in the form:

$$ \iiint_V \left( \partial_t (\rho \vec{v}) + \vec{\nabla} \cdot (\rho \vec{v} \otimes \vec{v}) + p \mathbb{I}_d \right) dV = 0 , $$

where $\mathbb{I}_d$ denotes the identity matrix. By considering that this relation holds for any control volume $V$, we obtain the *momentum conservation equation*:

$$ \partial_t (\rho \vec{v}) + \vec{\nabla} \cdot (\rho \vec{v} \otimes \vec{v}) + p \mathbb{I}_d = 0 . \qquad (2.5) $$

Once again, the transition from the integrated equation to the local equation is valid “almost everywhere”, in the sense of Lebesgue integrals.

### 2.1.3 Total energy conservation

The conservation law for the total energy $E_{tot}$ of a system is based on the first law of thermodynamics, which states that the variation of the total energy of the fluid within a control volume $V$ is equal to the sum of the power of the non-conservative external forces $P_e$ and the heating power exchanged with the fluid $P_c$. This relation can therefore be written in the following form:

$$ \frac{dE_{tot}}{dt} = P_e + P_c . \qquad (2.6) $$

First, the total energy is the sum of the energy contained within a control volume $V$ and can be expressed as:

$$ \frac{dE_{tot}}{dt} = \frac{d}{dt} \iiint_V E \, dV , $$
where E denotes the total energy density of the fluid. Analogously to what was presented in Section 2.1.1. *Mass conservation*, the time derivative of the total energy can be written as follows:

$$ \frac{dE_{tot}}{dt} = \iiint_V \left( \partial_t E + \vec{\nabla} \cdot (E\vec{v}) \right) dV . $$

The power term associated with non-conservative forces may include the power of the gravitational force, the power of the electric field, or the power of pressure forces. In the present case, we consider only the power of the pressure forces, which can be written as follows:

$$ \mathcal{P}_e = \iint_S -p\vec{v} \cdot \vec{n} \, dS , $$

which, according to the Green–Ostrogradsky theorem, can be rewritten in the following form:

$$ \mathcal{P}_e = - \iiint_V \vec{\nabla} \cdot (p \vec{v}) \, dV . $$

Concerning the heating power transferred to the fluid, one may account for thermal conduction as well as radiative energy losses or gains. For the moment, we neglect both phenomena ($\mathcal{P}_c = 0$), but we will later detail in Section 2.3.2. *Radiative hydrodynamics equations* how to include energy exchanges between the fluid and the radiation. Thus, the first law of thermodynamics can be rewritten as:

$$ \iiint_V \left( \partial_t E + \vec{\nabla} \cdot \left( (E + p)\vec{v} \right) \right) dV = 0 . $$

Since this integral is valid for any control volume $V$, we obtain the *total energy conservation equation*:

$$ \partial_t E + \vec{\nabla} \cdot \left( (E + p) \vec{v} \right) = 0 . \qquad (2.7) $$

Once again, the transition from the integrated equation to the local equation is valid “almost everywhere”, in the sense of Lebesgue integrals.

### 2.1.4 Euler's equations

Finally, we have obtained the conservation equations for mass (2.3), momentum (2.5), and total energy (2.7), which govern the evolution of a perfect fluid. These equations neglect the effects of gravity, electromagnetic forces, radiation, and thermal exchanges with the environment. They form a system of coupled equations known as the *Euler equations*, which play a central role in fluid mechanics. These equations depend on four fundamental quantities: the mass
density $\rho$, the velocity field $\vec{v}$, the pressure $p$, and the total energy density $E$, and can be summarised as follows:

$$ \left\{ \begin{aligned} \partial_t \rho + \vec{\nabla} \cdot (\rho \vec{v}) &= 0, \\ \partial_t (\rho \vec{v}) + \vec{\nabla} \cdot (\rho \vec{v} \otimes \vec{v}) + p \mathbb{I}_d &= 0, \\ \partial_t E + \vec{\nabla} \cdot ((E+p)\vec{v}) &= 0. \end{aligned} \right. \qquad (2.8) $$

However, we have only three equations to describe four quantities, which makes the system insufficiently constrained. To complete the system and make the model solvable, it is necessary to introduce a fourth equation that relates the missing variables and allows the system to be closed. To do so, we will express the total energy density of the fluid in terms of the other hydrodynamic quantities. Moreover, since temperature plays a central role in describing the interactions of the fluid with radiation, we will also detail how to compute the equation of state of the fluid, which relates the pressure to the fluid's temperature.

## ❖ Equation of state

Let us begin by establishing the expression of the equation of state, which relates the pressure to the temperature of the fluid. By neglecting interactions between the particles composing the gas, the partial pressure $p_i$ of each particle $i$ can be described using the ideal gas equation of state.

$$ p_i = n_i k_B T_i, \qquad (2.9) $$

where $k_B$ is the Boltzmann constant, $n_i$ is the number density of particle $i$, and $T_i$ is its temperature. The total pressure and total density can then be expressed respectively as follows:

$$ p = \sum_i p_i, $$

$$ \rho = \sum_i m_i n_i, $$

where $m_i$ is the mass of particle $i$. Assuming LTE, which implies a common temperature for all particles ($\forall i$, $T_i = T$), the global equation of state of the fluid becomes:

$$ p = \frac{\rho k_B T}{\mu m_H}, \qquad (2.10) $$

where $m_H$ is the mass of the hydrogen atom and $\mu m_H$ is the mean atomic weight of the medium. The quantity $\mu$ is defined by the relation:

$$ \mu = \frac{\sum_i n_i m_i}{\sum_i n_i m_H}. \qquad (2.11) $$

We may distinguish four different cases for computing $\mu$:
1. **Pure, non-ionised gas:** when a gas consists of particles of mass $m_a$, the mean atomic weight is given by:

$$ \mu = \frac{m_a}{m_H}, \qquad (2.12) $$

2. **Pure, ionised gas:** if the gas consists of neutral particles a, ions i, and electrons e, the mean atomic weight is:

$$ \mu = \frac{m_e n_e + m_i n_i + m_a n_a}{(n_e + n_i + n_a)m_H}, \qquad (2.13) $$

where $n_e$, $n_i$, and $n_a$ are the number densities of electrons, ions, and neutral atoms, respectively. Since the electron mass is very small compared to that of ions and neutral atoms, and since neutral atoms have approximately the same mass as ions ($m_e \ll m_i$, $m_a$ and $m_a \approx m_i$), one may neglect terms involving the electron mass and consider the mass of ions and neutral atoms to be equal. This simplifies the expression for $\mu$:

$$ \mu = \frac{m_i}{(1 + Z)m_H} = \frac{\mu^*}{(1 + Z)}, \qquad (2.14) $$

where $Z = n_e/(n_i + n_a)$ is the ionisation degree and $\mu^* = m_a/m_H$. Under LTE, Z depends only on the density and temperature of the medium. Consequently, the equation of state for an ionised gas is classically written as:

$$ p = (1 + Z) \frac{\rho k_B T}{\mu^* m_H}, \qquad (2.15) $$

3. **Non-reactive, non-ionised mixture:** no particular simplification exists for $\mu$, and $\mu m_H$ represents the average mass of the mixture's particles. The mean atomic weight $\mu$ remains constant over time and is computed using expression (2.11);

4. **Non-reactive, ionised mixture:** the equation of state becomes more complex due to the evolution of the species' ionisation degree, which depends on time and on the fluid state;

5. **Reactive, potentially ionised mixture:** the equation of state becomes even more complex, as it depends not only on time and on the fluid state, potentially accounting for ionisation, but also on ongoing chemical reactions. No simplification is possible for computing the mean atomic weight $\mu$. This case deviates from the LTE approximation and provides a more accurate description of non-equilibrium phenomena.

In this work, although the regimes studied in radiative hydrodynamics generally involve ionised gases, we restrict ourselves to the case of pure, non-ionised gases (case 1), for which the mean atomic weight $\mu$ remains constant throughout the flow. This choice is justified by the
main objective of the thesis, namely to improve the computation of radiation transport and to study its specific influence on the fluid dynamics.

## ❖ Total energy of the fluid

It is necessary to express the total energy of the fluid in order to relate it to the hydrodynamic quantities and to close the Euler equations (2.8). It consists of three main contributions: the macroscopic kinetic energy $E_{k,macro}$, the macroscopic potential energy $E_{p,macro}$, and the internal energy $u$:

$$E = E_{k,macro} + E_{p,macro} + u, \qquad (2.16)$$

Let us detail the physics involved in each of these terms:

* **Macroscopic kinetic energy:** It is associated with the bulk motion of the fluid and can be expressed as follows:

$$E_{k,macro} = \frac{1}{2} \rho ||\vec{v}||^2,$$

where $||\vec{v}||$ denotes the norm of the velocity field.

* **Macroscopic potential energy:** It arises from external forces acting on the fluid, such as gravitational or electromagnetic fields. In the present case, none of these effects are considered, and one has $E_{p,macro} = 0$.

* **Internal energy:** It represents the total energy within the gas and is composed of two contributions:

1. **Internal kinetic energy:** It can be expressed as a function of the gas temperature as follows:

$$E_{k,int} = \sum_{i} E_{k,i},$$

where $E_{k,i}$ is the average kinetic energy associated with particles $i$, and for a perfect gas it is given by:

$$E_{k,i} = \frac{\ell_i n_i k_B T_i}{2},$$

where $\ell_i$ is the number of degrees of freedom of particle $i$ (3 for a monoatomic particle, 5 for a diatomic particle). In the same way as in the equation of state,
assuming that all particles have the same temperature, the total internal energy can be written as:

$$E_{k,int} = \rho \frac{\ell k_B}{2\mu m_H} T = \rho c_v T,$$

where $\ell = (\sum_i \ell_i n_i)/(\sum_i n_i)$ is the average number of degrees of freedom of the particles, $c_v$ is the specific heat capacity at constant volume, and $\mu$ is the same quantity as detailed in the section Equation of state.

2. **Internal potential energy:** It accounts for chemical binding energies, particle state transitions, or nuclear transitions. In the case of a pure ionized gas, this energy corresponds to the ionization energy of the particles. In the present case, particle ionization is not taken into account, so this energy is zero.

Thus, in the case where particle ionization, chemical transformations, or particle state transitions are not taken into account, the internal energy is given by:

$$u = \rho \frac{\ell k_B}{2\mu m_H} T = \rho c_v T.$$

Lets note that the enthalpy of the fluid can thus be expressed as:

$$h = u + p = \rho \frac{(\ell + 2)k_B}{2\mu m_H} T = \rho c_p T,$$

where $c_p$ is the specific heat capacity at constant pressure. The adiabatic index $\gamma$ is defined as the ratio $c_p/c_v$, which can be expressed as a function of the number of degrees of freedom of the particles.

$$\gamma = c_p / c_v = \frac{\ell + 2}{\ell}.$$

The specific heat capacity at constant volume can then be written as:

$$c_v = \frac{k_B}{(\gamma - 1)\mu m_H}.$$

Based on considerations of the internal degrees of freedom of the particles, the adiabatic index takes the value $\gamma = 5/3$ for a monoatomic gas, and $\gamma = 7/5 = 1.4$ for a diatomic gas. By contrast, for fluids composed of more complex particles, mixtures of different species, real or non-ideal gases, or partially or fully ionized plasmas, the value of $\gamma$ may deviate from these two standard values. It may even evolve in time, depending on the local hydrodynamic conditions of the fluid, such as temperature or density, which directly influence the degrees of freedom that are effectively excited.
Finally, the total energy considered here can be written as:

$$E = \frac{1}{2} \rho |\vec{v}|^2 + \rho c_v T, \qquad (2.17)$$

or, equivalently, using the equation of state developed in the previous section and the expression of $c_v$:

$$E = \frac{1}{2} \rho |\vec{v}|^2 + \frac{p}{\gamma - 1}. \qquad (2.18)$$

## 2.2 Radiative transfer

In astrophysics, flows interact with radiation in different ways depending on the properties of the medium. In many cases, radiation propagates freely through space and acts only as an energy loss mechanism for the fluid. Such media are referred to as optically thin. In this situation, cooling functions can be introduced to model these energy losses. This framework is particularly relevant for supernova remnants [71, 72, 73, 10] as well as for accretion processes [74, 75].

Conversely, some media exhibit very rapid interactions between radiation and matter, making the radiative equilibrium approximation relevant. This regime is typical of star formation [76], stellar pulsations [77, 78, 79], and protoplanetary disks [80].

However, these two extreme cases are not sufficient to describe the full complexity of radiation-matter interactions. There exist intermediate situations in which radiation interacts significantly with matter, but too slowly to reach local thermal equilibrium with its surroundings. This phenomenon is encountered in particular at stellar surfaces, where radiation undergoes a transition between an optically thick medium (the stellar atmosphere) and an optically thin medium (the interstellar medium), as well as in protoplanetary disks, which are optically thick in their interior but allow radiation to escape in certain regions [81]. Moreover, these simplified descriptions neglect the dependence of radiative behavior on photon frequency, which limits the accuracy of predictions of observable radiative fluxes. To overcome these limitations, we develop in this section a model that enables a more realistic description of radiation propagation and its interactions with matter, while remaining within the framework of LTE.

Finally, to simplify the mathematical expressions, we explicitly omit the spatio-temporal dependencies on $\vec{x}$ and $t$, although all the quantities considered generally vary in time and space.
<figure><img src="image_8.png" /><figcaption>Figure 2.1: Schematic representation of a radiation beam crossing the surface $d\vec{S}$ and propagating within the solid angle $d\Omega$.</figcaption></figure>

### 2.2.1 Radiative transfer equation

#### ❖ Specific intensity

In order to describe radiation, it is essential to define quantities that characterize the radiation field. The most important one is the specific intensity $I_v(\vec{n})$, which is evaluated at a time $t$, at a position $\vec{x}$, and describes radiation propagating in a direction $\vec{n}$ (see figure 2.1). The amount of radiative energy with frequencies between $\nu$ and $\nu + d\nu$, crossing the surface $d\vec{S}$, and propagating within a solid angle $d\Omega$, at position $\vec{x}$ and time $t$, depends on the specific intensity as follows:

$$dE_v = I_v(\vec{n}) \vec{n} \cdot d\vec{S} d\Omega d\nu dt \quad (2.19)$$

Here, the vector $d\vec{S}$ is oriented along the normal to the considered surface. It is important to note that, in the absence of interactions with matter, the specific intensity remains constant along the direction of propagation of the radiation. By contrast, when interactions with matter occur, radiation may lose or gain energy, leading to a variation of the specific intensity. It is therefore necessary to model these interactions in order to correctly describe the evolution of the radiation intensity within the medium.

#### ❖ Absorption, emission and scattering

When radiation propagates through a material medium, it can interact with matter in three different manners, thereby modifying its behavior: through photon absorption, emission, or
<figure><img src="image_9.png" /><figcaption>Figure 2.2: Types of interactions between light and matter.</figcaption></figure>

scattering.

1. **Photon absorption (figure 2.2a):** photons from the radiation field can be absorbed by matter, transferring their energy to the medium. This can increase the temperature of the material or induce changes in the state of matter, such as molecular dissociation or ionization;

2. **Photon emission (figure 2.2b):** after absorbing energy, matter can re-emit photons. This may occur during energy transitions in the atoms or molecules of the medium, thereby releasing photons;

3. **Photon scattering (figure 2.2c):** photons can also be deflected from their original trajectory as they interact with particles in the medium. Scattering can be elastic, when the photon energy remains unchanged, or inelastic, when it is modified.

These interactions affect both the intensity and the direction of the radiation and must be taken into account to model radiative transport through a medium. Let us consider radiation propagating through a layer of matter of thickness $d\ell$, over which photons are absorbed. The amount of specific intensity lost is:

$$dI_v = -\chi_v(\vec{n})I_v d\ell, \qquad (2.20)$$

where $\chi_v(\vec{n})$ is called the absorption coefficient and is composed of two components: *pure absorption* $\kappa_v(\vec{n})$ and *pure scattering* $\sigma_v(\vec{n})$,

$$\chi_v(\vec{n}) = \kappa_v(\vec{n}) + \sigma_v(\vec{n}), \qquad (2.21)$$
Matter can also emit photons as the light beam traverses a layer of material of thickness $d\ell$. This gain in photons is modeled by the following equation:

$$dI_v = \eta_v(\vec{n}) d\ell, \qquad (2.22)$$

where $\eta_v(\vec{n})$ is called the emission coefficient and is also composed of two parts: thermal emission $\eta_v^t(\vec{n})$ and scattering emission $\eta_v^s(\vec{n})$,

$$\eta_v(\vec{n}) = \eta_v^t(\vec{n}) + \eta_v^s(\vec{n}). \qquad (2.23)$$

All these coefficients provide a general description of the interactions between matter and radiation. Their values are determined by four fundamental physical processes [82, 83]:

1. **Free-free processes (bremsstrahlung):** a free electron, accelerated or decelerated through interaction with an ion or another electron, emits a photon, producing continuous radiation. Conversely, a free electron may absorb a photon in the presence of an ion, thereby increasing its kinetic energy;

2. **Bound-free processes (photoionization):** an atom or ion can absorb a photon with sufficient energy to eject an electron, leading to ionization of the particle;

3. **Free-bound processes (recombination):** conversely, a free electron can be captured by an ion or an atom, resulting in the emission of a photon whose energy equals the difference between the initial and final energy states of the electron;

4. **Bound-bound processes:** an atom or molecule may absorb or emit a photon when one of its electrons changes energy level. These processes are responsible for the characteristic emission and absorption spectra of the elements.

In a medium at rest, the absorption coefficient $\kappa_v$ and the thermal emission coefficient $\eta_v^t$ are isotropic. When a photon is scattered, the energy absorbed through the pure scattering coefficient $\sigma_v$ is re-emitted via the scattering emission coefficient $\eta_v^s$. These two coefficients are therefore related. In general, determining $\eta_v^s$ is not trivial, as this coefficient strongly depends on the microscopic properties of the particles present. However, under certain simplifying assumptions, it is possible to derive an analytical expression. These assumptions are:

1. The pure scattering coefficient $\sigma_v$ is isotropic (which is often the case);

2. Photons are scattered isotropically;
3. Scattering is elastic and does not modify the photon energy,

The analytical expression of this coefficient is then:

$$ \eta_v^s(\vec{n}) = \frac{\sigma_v}{4\pi} \int_{4\pi} I_v(\vec{n}') d\Omega' . $$

Furthermore, if radiation is assumed to be in LTE, Kirchhoff's law allows the thermal emission coefficient to be expressed as follows:

$$ \eta_v'(\vec{n}) = \kappa_v B_v(T) , $$

where $B_v(T)$ is the Planck function, given by:

$$ B_v(T) = \frac{2h\nu^3}{c^2} \left[ e^{-\frac{h\nu}{k_B T}} - 1 \right]^{-1} , $$

where $h$ is Planck's constant, $c$ is the speed of light, $\nu$ is the radiation frequency, and $T$ is the temperature. Finally, the total emission coefficient can be written entirely in terms of the absorption coefficients as:

$$ \eta_v = \eta_v^t + \eta_v^s = \kappa_v B_v(T) + \frac{\sigma_v}{4\pi} \int_{4\pi} I_v(\vec{n}') d\Omega' . \quad (2.24) $$

## ❖ Expression of the radiative transfer equation

The radiative transfer equation describes the evolution of the specific intensity of radiation as it propagates through a material medium. It expresses the balance between emission and absorption processes that modify this intensity. This balance can be written in the following form:

$$ dI_v(\vec{n}) = \eta_v(\vec{n}) d\ell - \chi_v(\vec{n}) I_v(\vec{n}) d\ell , $$

or, equivalently, in Cartesian coordinates:

$$ \left( \frac{1}{c} \partial_t + \vec{n} \cdot \vec{\nabla} \right) I_v(\vec{n}) = \eta_v(\vec{n}) - \chi_v(\vec{n}) I_v(\vec{n}) . \quad (2.25) $$

Taking into account all the assumptions previously mentioned regarding radiation-matter interactions, the radiative transfer equation can also be written as:

$$ \left( \frac{1}{c} \partial_t + \vec{n} \cdot \vec{\nabla} \right) I_v(\vec{n}) = \kappa_v (B_v(T) - I_v(\vec{n})) + \sigma_v \left( \frac{1}{4\pi} \int_{4\pi} I_v(\vec{n}') d\Omega' - I_v(\vec{n}) \right) . \quad (2.26) $$
### 2.2.2 The moment equations

In practice, the specific intensity of radiation is not directly observable. Instead, one has access to integrated quantities, such as the radiative flux, which corresponds to the amount of energy transported by radiation through a given surface, or the radiative pressure, which represents the mechanical force exerted by light on matter. This effect plays a key role in a wide range of phenomena, such as the dynamics of orbiting satellites or the operation of solar sails. In order to reformulate the radiative transfer equation (2.25) in a more directly usable form, we therefore seek to express it in terms of three fundamental quantities: radiative energy, radiative flux, and radiative pressure.

The monochromatic radiative energy represents the energy density of photons of frequency $v$ at a given point in space. It corresponds to the zeroth-order moment of the specific intensity and is expressed as:

$$E_v = \frac{1}{c} \int_{4\pi} I_v(\vec{n}) \, d\Omega . \qquad (2.27)$$

The monochromatic radiative flux quantifies the energy transported by photons of frequency $v$ through a given surface. It corresponds to the first-order moment of the specific intensity and is written as:

$$\vec{F}_v = \int_{4\pi} I_v(\vec{n}) \vec{n} \, d\Omega . \qquad (2.28)$$

Finally, the monochromatic radiative pressure represents the stress exerted by radiation on matter when it interacts with photons of frequency $v$. It corresponds to the second-order moment of the specific intensity and is expressed as:

$$P_v = \int_{4\pi} I_v(\vec{n}) \vec{n} \otimes \vec{n} \, d\Omega . \qquad (2.29)$$

To obtain the monochromatic zeroth-order moment equation (referred to as the *radiation energy equation*), one integrates the radiative transfer equation (2.25) over all solid-angle directions. Using definitions (2.27) and (2.28), this equation can be written as:

$$\partial_t E_v + \vec{\nabla} \cdot \vec{F}_v = -c S_v^0 . \qquad (2.30)$$

To obtain the monochromatic first-order moment equation (referred to as the *radiation momentum equation*), the radiative transfer equation is multiplied by $\vec{n}/c$ and integrated over all solid-angle directions. This equation reads:

$$\partial_t (c^{-2} \vec{F}_v) + \vec{\nabla} \cdot P_v = - \vec{S}_v . \qquad (2.31)$$
The terms $S_v^0$ and $\vec{S}_v$ correspond to source terms modeling the interactions between light and matter and are given by:

$$S_v^0 = -c^{-1} \int_{4\pi} (\eta_v(\vec{n}) - \chi_v(\vec{n})I_v(\vec{n})) d\Omega, \qquad (2.32)$$

$$\vec{S}_v = -c^{-1} \int_{4\pi} (\eta_v(\vec{n}) - \chi_v(\vec{n})I_v(\vec{n})) \vec{n} d\Omega . \qquad (2.33)$$

Working in the comoving frame and adopting the assumptions stated previously regarding radiation-matter interactions, these source terms can be written as:

$$S_v^0 = \kappa_v \left( E_v - \frac{4\pi}{c} B_v(T) \right), \qquad (2.34)$$

$$\vec{S}_v = \chi_v \vec{F}_v / c . \qquad (2.35)$$

The system of equations (2.30) and (2.31) describes the full set of radiative phenomena and is equivalent to the radiative transfer equation (2.25). However, it is underdetermined, as it involves three quantities for only two equations. To close the system, it is therefore necessary to introduce an additional relation, referred to as a *closure relation*. Several models propose different approaches to define this relation.

## 1. The isotropic Flux-Limited Diffusion model [84]:

In this model, radiation is assumed to be close to radiative equilibrium. Only the radiative energy equation (2.30) is retained, while the evolution of the radiative flux is neglected, leading to the approximation $\vec{F}_v = c/\chi_v \vec{\nabla} \cdot \mathbb{P}_v$. Assuming isotropic radiation, the radiative pressure is written as $\mathbb{P}_v = 1/3E_v \mathbb{I}_d$. This model is the fastest to solve numerically, but it does not guarantee the fundamental physical constraint on the radiative flux, namely $||\vec{F}_v|| \le cE_v$. An improvement consists in introducing a flux limiter $\lambda$, which modifies the expression of the radiative flux such that $\vec{F}_v = c\lambda/(3\chi_v) \vec{\nabla} E_v$, where the factor $\lambda$ is chosen to ensure that the physical constraint $||\vec{F}_v|| \le cE_v$ is satisfied.

## 2. The P1 model [85]:

This model still assumes that radiation is close to radiative equilibrium, while retaining both the radiative energy equation (2.30) and the radiative flux equation (2.31). It imposes a closure based on the assumption of isotropic radiation, $\mathbb{P}_v = 1/3E_v \mathbb{I}_d$. Unlike the previous model, it guarantees the inequality $||\vec{F}_v|| \le cE_v$ and allows for more general situations to be treated, but it loses accuracy when the radiation exhibits strong anisotropy (radiative equilibrium not valid anymore).
### 3. The M1 model [86, 87]:

This model generalizes the description of radiative transfer by relaxing the radiative equilibrium assumption. It retains both the radiative energy and flux equations (2.30) and (2.31), while allowing for a certain degree of radiation anisotropy. Assuming that the radiation has an axis of symmetry aligned with the direction of the radiative flux, the radiative pressure can be expressed as:

$$P_v = D_v E_v, \qquad (2.36)$$

where $D_v$ is called the *Eddington tensor* and is defined by the relation [86]:

$$D_v = \frac{1 - \chi_{R,v}}{2} \mathbb{I}_d + \frac{3\chi_{R,v} - 1}{2} \frac{\vec{F}_v \otimes \vec{F}_v}{||\vec{F}_v||^2}, \qquad (2.37)$$

with $\chi_{R,v}$ the *Eddington factor*, whose expression depends on the adopted assumptions. Several authors have proposed different formulations for $\chi_{R,v}$:

- **Kershaw (1976) [88]** proposed a second-order polynomial expression, which is simple but does not account for the physical nature of the radiation field;

- **Minerbo (1978) [89]** proposed an approach based on the maximization of radiative entropy, assuming that $E_v \ll 4\pi h v^3/c^3$;

- **Levermore (1979, 1981) [90, 91]** relied on the isotropic diffusion model;

- **Pomraning (1973) [92]** developed a formulation adapted to relativistic fluids and applicable to radiative quantities integrated over all frequencies.

This model satisfies the inequality $||\vec{F}_v|| \le cE_v$ and efficiently describes situations in which radiation is strongly anisotropic, making it one of the most accurate approaches. However, in 2D or 3D, its assumption of symmetry around the flux axis limits its ability to correctly represent complex configurations, such as two intersecting light beams<fnref n="1" />. Despite this limitation, this model will be used, and its features will be detailed in the following sections.

From a numerical standpoint, it is impossible to represent radiation propagation at all frequencies $v$, as this would require solving an infinite number of equations. An approach that enables the numerical solution of these equations consists in discretizing the electromagnetic spectrum into frequency “groups”. This corresponds to the M1-multigroup model.

<footnote n="1">¹In this scenario, the M1 model predicts that the two light beams “collide”, which is not physically realistic.</footnote>
<figure><img src="image_10.png" /><figcaption>Figure 2.3: Electromagnetic spectrum. Example of discretization into $\mathcal{G}$ groups: the 1st group includes gamma rays, the 2nd includes X-rays, and the last includes long-wavelength radiation.</figcaption></figure>

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

### ❖ The M1-multigroup model

The M1-multigroup model relies on a spectral discretization of radiation into $\mathcal{G}$ frequency groups. This approach makes it possible to capture the spectral variability of radiative transfer by considering several frequency bands separately. Before addressing the detailed formulation of the model, it is necessary to define the structure of the frequency groups. One thus introduces a set of $\mathcal{G}$ frequency intervals $\{[\nu_{g-1/2}, \nu_{g+1/2}]\}_{g \in \llbracket 1, \mathcal{G} \rrbracket}$ covering the entire electromagnetic spectrum, with boundary conditions $\nu_{1/2} = 0$ and $\nu_{\mathcal{G}+1/2} = +\infty$. Each group $g \in \llbracket 1, \mathcal{G} \rrbracket$ is defined by a pair of frequencies $(\nu_{g-1/2}, \nu_{g+1/2})$, which delimit the corresponding frequency range. Figure 2.3 illustrates an example of a partition of the electromagnetic spectrum into several groups: the first group describes the behavior of photons in the gamma-ray domain, the second in the X-ray domain, and so on. The overall behavior of radiation within a given group is then described from the contributions of photons with frequencies lying between these two bounds. This model, developed by Turpault [93, 94], makes it possible to account for the spectral behavior of radiation within a numerical simulation code. It constitutes an improvement over the M1-gray model, which treats radiation as a whole without distinguishing between different spectral components. The latter is in fact a particular case of the M1-multigroup model, in which only a single group covering the entire spectrum is considered, with frequency bounds $\nu_{1/2} = 0$ and $\nu_{3/2} = +\infty$.

Before detailing the equations of the model, we introduce a few fundamental definitions that will be used in this section. Let us begin by defining the Planck and Rosseland mean
opacities associated with each frequency group g, according to the following expressions:

$$ \kappa_{P,g} = \left[ \int_{\nu_{g-1/2}}^{\nu_{g+1/2}} \kappa_{\nu} B_{\nu}(T) \, d\nu \right] \left[ \int_{\nu_{g-1/2}}^{\nu_{g+1/2}} B_{\nu}(T) \, d\nu \right] , \qquad (2.38) $$

$$ \kappa_{R,g} = \left[ \int_{\nu_{g-1/2}}^{\nu_{g+1/2}} \left( \frac{\partial B_{\nu}(T)}{\partial T} \right) d\nu \right] \left[ \int_{\nu_{g-1/2}}^{\nu_{g+1/2}} \chi_{\nu}^{-1} \left( \frac{\partial B_{\nu}(T)}{\partial T} \right) d\nu \right] . \qquad (2.39) $$

Let us also define the operator $< \cdot >_g$ as follows:

$$ < \cdot >_g = c^{-1} \int_{\nu_{g-1/2}}^{\nu_{g+1/2}} \int_{4\pi} \cdot d\Omega \, d\nu . \qquad (2.40) $$

The M1-multigroup model aims to describe the global behavior of radiation within each spectral group g (see figure 2.3). To this end, $\mathcal{G}$ sets of radiative quantities associated with group g are introduced, including the energy $E_g$, the flux $\vec{F}_g$, and the radiative pressure $\vec{P}_g$, defined as follows:

$$ E_g = \int_{\nu_{g-1/2}}^{\nu_{g+1/2}} E_{\nu} \, d\nu , \qquad (2.41) $$

$$ \vec{F}_g = \int_{\nu_{g-1/2}}^{\nu_{g+1/2}} \vec{F}_{\nu} \, d\nu , \qquad (2.42) $$

$$ \vec{P}_g = \int_{\nu_{g-1/2}}^{\nu_{g+1/2}} \vec{P}_{\nu} \, d\nu . \qquad (2.43) $$

Moreover, we define the total radiative energy $E_R$, the total radiative flux $\vec{F}_R$, and the total radiative pressure $\vec{P}_R$ as:

$$ E_R = \int_0^{+\infty} E_\nu \, d\nu = \sum_{g=1}^{\mathcal{G}} E_g , \qquad (2.44) $$

$$ \vec{F}_R = \int_0^{+\infty} \vec{F}_{\nu} \, d\nu = \sum_{g=1}^{\mathcal{G}} \vec{F}_g , \qquad (2.45) $$

$$ \vec{P}_R = \int_0^{+\infty} \vec{P}_{\nu} \, d\nu = \sum_{g=1}^{\mathcal{G}} \vec{P}_g . \qquad (2.46) $$

By integrating the monochromatic radiative transfer equations (2.30) and (2.31) over the frequency interval $[\nu_{g-1/2}; \nu_{g+1/2}]$, one obtains the radiative transfer equations within group g:

$$ \forall g \in [1, \mathcal{G}], \quad \left\{ \begin{array}{l} \partial_t E_g + \vec{\nabla} \cdot \vec{F}_g = -cS_g^0 , \\ \partial_t (c^{-2}\vec{F}_g) + \vec{\nabla} \cdot \vec{P}_g = -\vec{S}_g , \end{array} \right. \qquad (2.47) $$
where $S_g^0$ and $\vec{S}_g$ are the source terms describing interactions between photons of each group g and matter. Their general expressions are:

$$S_g^0 = - \langle (\eta_v(\vec{n}) - \chi_v(\vec{n})I_v(\vec{n})) \rangle_g, \qquad (2.48)$$

$$\vec{S}_g = - \langle (\eta_v(\vec{n}) - \chi_v(\vec{n})I_v(\vec{n})) \vec{n} \rangle_g . \qquad (2.49)$$

Adopting the assumptions previously introduced regarding radiation-matter interactions, these terms can be written as:

$$S_g^0 = \kappa_{E,g} E_g - \kappa_{P,g} a_R \theta_g^4(T), \qquad (2.50)$$

$$\vec{S}_g^i = \chi_{F,g}^i \vec{F}_g^i / c, \qquad (2.51)$$

where $a_R$ is the radiation constant and $\theta_g$ represents, in a sense, the radiative temperature of group g at radiative equilibrium, and is given by:

$$\theta_g(T) = \left( \frac{\langle B_v(T) \rangle_g}{a_R} \right)^{1/4} . \qquad (2.52)$$

The opacities $\kappa_{E,g}$ and $\chi_{F,g}^i$ are energy- and flux-weighted opacities for group g and are defined as:

$$\kappa_{E,g} = \left[ \int_{\nu_{g-1/2}}^{\nu_{g+1/2}} \kappa_v E_v \, d\nu \right] \left[ \int_{\nu_{g-1/2}}^{\nu_{g+1/2}} E_v \, d\nu \right], \qquad (2.53)$$

$$\chi_{F,g}^i = \left[ \int_{\nu_{g-1/2}}^{\nu_{g+1/2}} \chi_v \vec{F}_v^i \, d\nu \right] \left[ \int_{\nu_{g-1/2}}^{\nu_{g+1/2}} \vec{F}_v^i \, d\nu \right] . \qquad (2.54)$$

It is clear that directly computing the mean opacities $\kappa_{E,g}$ and $\chi_{F,g}^i$ is difficult, as their evaluation requires solving the spectral moment equations (2.30) and (2.31), as well as precise knowledge of the spectral opacities of matter. This approach, in addition to being complex, would significantly increase the computational cost associated with solving the M1-multigroup model (2.47). To simplify this issue, Mihalas [95] proposed replacing $\kappa_{E,g}$ and $\chi_{F,g}^i$ with the Planck and Rosseland mean opacities, $\kappa_{P,g}$ and $\kappa_{R,g}$, respectively. In the present framework, we assume that the fluid is in LTE, which implies that these opacities depend only on the temperature and density of the fluid.

This approximation is justified near the radiative equilibrium regime, where the radiative energy tends toward the Planck function ($E_v \propto B_v(T)$), and where the radiative flux can be expressed as $\vec{F}_v \approx c/\chi_v \vec{\nabla} \cdot \mathbb{P}_v \propto \partial B_v(T)/(\chi_v \partial T)$. It is particularly well suited to optically thick media, in which radiation efficiently thermalizes with matter.
However, this approximation becomes less valid in optically thin regimes, where radiation is more directional and departs from radiative equilibrium. Nevertheless, within the M1-multipgroup framework, the finer the spectral groups, the more the mean opacities converge toward the spectral opacities: one has $\kappa_{P,g}$, $\kappa_{E,g} \to \kappa_v$ and $\kappa_{R,g}$, $\chi_{F,g}^i \to \chi_v$, which reduces the impact of this approximation. In this context, the source terms then take the following form:

$$S_g^0 = \kappa_{P,g}(E_g - a_R \theta_g^4(T)), \qquad (2.55)$$

$$\vec{S}_g^i = \kappa_{R,g} \vec{F}_g^i / c, \qquad (2.56)$$

It remains to define the closure relation for each group in order to solve the system. In the present approach, as before, radiation is assumed to be symmetric around the direction of the radiative flux $\vec{F}_g$ within each group. Under this assumption, the closure relation takes a form analogous to that given in (2.36):

$$P_g = D_g E_g, \qquad (2.57)$$

where $D_g$ is the *Eddington tensor* associated with group g and is given by [86]:

$$D_g = \frac{1 - \chi_g}{2} \mathbb{I}_d + \frac{3\chi_g - 1}{2} \frac{\vec{F}_g \otimes \vec{F}_g}{||\vec{F}_g||^2}, \qquad (2.58)$$

where $\chi_g$ denotes the *Eddington factor* associated to the group g. Its determination requires formulating an assumption on the shape of the specific intensity in order to render the problem mathematically tractable. In the following section, we present the developments that enable its numerical computation.

### 2.2.3 The Eddington factor

Within the framework of the M1-multipgroup model, the objective is to determine an expression for the specific intensity that minimizes the radiative entropy, in order to close the system of equations for each group and to derive an expression for the Eddington factor. The radiative entropy density is given by [96]:

$$h_R(I) = \frac{2k_B\nu^2}{c^3} [n \ln(n) - (n+1) \ln(n+1)], \qquad (2.59)$$

where n is the occupation number, defined by $n = \frac{c^2}{2h\nu}I$. The total radiative entropy $H_R(I)$ is then defined as

$$H_R(I) = \sum_{g=1}^{G} \langle h_R(I) \rangle_g, \qquad (2.60)$$
As seen in equations (2.41) and (2.42), the specific intensity is related to the radiative energy and flux in each group. We therefore seek to retain only the specific intensity $I_v$ that minimizes the radiative entropy, i.e. the specific intensity that is the solution of the following minimization problem:

$$H_R(I_v) = \min_I \left\{ H_R(I) \mid \forall g \in [1, G], \langle I \rangle_g = E_R, \langle c \vec{n} I \rangle_g = \vec{F}_g \right\} . \quad (2.61)$$

It can be shown that the intensity satisfying this problem can be written as [97]:

$$I_v = \sum_{g=1}^{G} \mathbb{1}_g(v) I_{v,g}, \quad (2.62)$$

where $\mathbb{1}_g(v) : \mathbb{R}^+ \to \{0, 1\}$ is the characteristic function, equal to 1 for frequencies $v$ in $[\nu_{g-1/2}; \nu_{g+1/2}]$ and 0 otherwise, and where $I_{v,g}$ denotes the specific intensity in group $g$, given by:

$$I_{v,g} = \frac{2h\nu^3}{c^2} \left[ \exp\left(\frac{h\nu}{k_B} \overrightarrow{m}_g \cdot \overrightarrow{\alpha}_g\right) - 1 \right]^{-1}, \quad (2.63)$$

where $\overrightarrow{\alpha}_g = (\alpha_{0,g}, \overrightarrow{\alpha}_{1,g})$ are the Lagrange multipliers of group $g$ associated with problem (2.61), $\overrightarrow{m}_g$ is a four-vector defined as $\overrightarrow{m}_g = (1, \overrightarrow{n}_g)$, and $\overrightarrow{n}_g$ is the propagation direction of light in group $g$. By setting $\overrightarrow{\alpha}_{1,g} \cdot \overrightarrow{n}_g = \alpha_{0,g}\beta_g \cos(\theta_g)$, the specific intensity in group $g$ can also be written as:

$$I_{v,g} = \frac{2h\nu^3}{c^2} \left[ \exp\left(\frac{h\nu}{k_B}\alpha_{0,g}(1 + \beta_g \cos(\theta_g))\right) - 1 \right]^{-1}, \quad (2.64)$$

In this context, $\alpha_{0,g}$ represents the inverse of a temperature and takes values in $\mathbb{R}^+$, while $\beta_g$ measures the degree of anisotropy of the radiation in group $g$, with values in the interval $]-1; 1[$. It is important to note that if $\alpha_{0,g}$ and $\beta_g$ are identical for all groups (i.e. if $\forall g \in [1, G]$, $\alpha_{0,g} = \alpha_0$ and $\beta_g = \beta$), one recovers the specific intensity of the M1-gray model.

However, unlike the M1-gray model, no general analytical formula exists for the Eddington factor $\chi_g$ in this framework. Nevertheless, Turpault showed that, in the M1-multigroup model, for any physically admissible pair $(E_g, ||\vec{F}_g||)$, there exists a unique pair of Lagrange multipliers $(\alpha_{0,g}, \beta_g)$ [97]. Consequently, the Eddington factor can be computed numerically using a root-finding algorithm, such as the bisection-Newton method or a line-search method, described in Appendix B. *Numerical calculation of the radiative quantities*.

In the next section, I will analyze the behavior of the Eddington factor in the M1-multigroup model, in accordance with the study presented in the corresponding article [98]. In addition, in Section 3.2. *Eddington factor in the M1-multigroup model*, I will present a method based on AI, developed during this thesis, that enables an efficient computation of $\chi_g$.
To begin with, let us examine the parameters that influence the Eddington factor in the M1-multigroup model. The radiative energy $E_g$, the norm of the radiative flux $||\vec{F}_g||$, and the radiative pressure component $P_g = \chi_g E_g$ in the direction of the radiative flux can be expressed as follows:

$$E_g = \frac{1}{c} \int_{\nu_{g-1/2}}^{\nu_{g+1/2}} \int_{4\pi} \mathcal{I}_v \, d\Omega \, d\nu, \qquad (2.65)$$

$$||\vec{F}_g|| = \int_{\nu_{g-1/2}}^{\nu_{g+1/2}} \int_{4\pi} \mathcal{I}_v \cos(\theta) \, d\Omega \, d\nu, \qquad (2.66)$$

$$P_g = \frac{1}{c} \int_{\nu_{g-1/2}}^{\nu_{g+1/2}} \int_{4\pi} \mathcal{I}_v \cos^2(\theta) \, d\Omega \, d\nu, \qquad (2.67)$$

or, equivalently, by noting that $d\Omega = \sin(\theta) d\theta d\phi$, setting $\mu = \cos(\theta)$, and using the expression of the specific intensity in group g, one can write:

$$E_g = \frac{4h\pi}{c^3} \int_{\nu_{g-1/2}}^{\nu_{g+1/2}} \int_{-1}^{1} \nu^3 \left[ e^{\frac{h\nu}{k_B}\alpha_{0,g}(1+\beta_g\mu)} - 1 \right]^{-1} d\mu \, d\nu, \qquad (2.68)$$

$$||\vec{F}_g|| = \frac{4h\pi}{c^2} \int_{\nu_{g-1/2}}^{\nu_{g+1/2}} \int_{-1}^{1} \nu^3 \mu \left[ e^{\frac{h\nu}{k_B}\alpha_{0,g}(1+\beta_g\mu)} - 1 \right]^{-1} d\mu \, d\nu, \qquad (2.69)$$

$$P_g = \frac{4h\pi}{c^3} \int_{\nu_{g-1/2}}^{\nu_{g+1/2}} \int_{-1}^{1} \nu^3 \mu^2 \left[ e^{\frac{h\nu}{k_B}\alpha_{0,g}(1+\beta_g\mu)} - 1 \right]^{-1} d\mu \, d\nu. \qquad (2.70)$$

It is then possible to distinguish four particular cases depending on the spectral bounds of the group under consideration:

1. $\nu_{g-1/2} = 0$ and $\nu_{g+1/2} = \infty$: this corresponds to the M1-gray model, in which the Eddington factor $\chi_R$ admits an analytical expression [92]:

$$\chi_R = \frac{3 + 4 f_R^2}{5 + 2\sqrt{4 - 3 f_R^2}}, \qquad (2.71)$$

where $f_R$ is the reduced flux, defined as $f_R = ||\vec{F}_R||/(c E_R)$, with $\vec{F}_R$ and $E_R$ denoting the radiative flux and radiative energy integrated over the entire electromagnetic spectrum, respectively;

2. $\nu_{g-1/2} = 0$ and $\nu_{g+1/2} < \infty$: Eddington factor for the first group;

3. $\nu_{g-1/2} > 0$ and $\nu_{g+1/2} = \infty$: Eddington factor for the last group;

4. $\nu_{g-1/2} > 0$ and $\nu_{g+1/2} < \infty$: general Eddington factor for all remaining groups.
In the following sections, I will analyze the parameters that influence the Eddington factor in each of these cases.

## ❖ Case $v_{g-1/2} > 0$ and $v_{g+1/2} < \infty$

In this case, it is possible to perform the change of variable $\tilde{v} = v/v_{m,g}$, where $v_{m,g}$ is the geometric mean of the frequency bounds of the group ($v_{m,g} = \sqrt{v_{g-1/2} v_{g+1/2}}$), in integrals (2.68), (2.69), and (2.70). The radiative quantities can then be expressed as follows:

$$ 
\begin{aligned}
\frac{k_B^4 E_g}{h^4 v_{m,g}^4} &= a_R \frac{15}{\pi^4} \int_{\delta_g^{1/2}}^{\delta_g^{-1/2}} \int_{-1}^{1} \tilde{v}^3 [e^{\tilde{v}} \tilde{\alpha}_{0,g} (1+\beta_g \mu) - 1]^{-1} d\mu \, d\tilde{v} \, , \\
\frac{k_B^4 ||\vec{F}_g||}{h^4 v_{m,g}^4} &= a_R c \frac{15}{\pi^4} \int_{\delta_g^{1/2}}^{\delta_g^{-1/2}} \int_{-1}^{1} \tilde{v}^3 \mu [e^{\tilde{v}} \tilde{\alpha}_{0,g} (1+\beta_g \mu) - 1]^{-1} d\mu \, d\tilde{v} \, , \\
\frac{k_B^4 P_g}{h^4 v_{m,g}^4} &= a_R \frac{15}{\pi^4} \int_{\delta_g^{1/2}}^{\delta_g^{-1/2}} \int_{-1}^{1} \tilde{v}^3 \mu^2 [e^{\tilde{v}} \tilde{\alpha}_{0,g} (1+\beta_g \mu) - 1]^{-1} d\mu \, d\tilde{v} \, ,
\end{aligned}
$$

where $\delta_g = v_{g-1/2}/v_{g+1/2}$ is a parameter taking values in $]0; 1[$, referred to as the group narrowness, and $\tilde{\alpha}_{0,g} = h v_{m,g} \alpha_{0,g} / k_B$ is the Lagrange multiplier $\alpha_{0,g}$ made dimensionless using the frequency $v_{m,g}$. In other words, by introducing the reduced flux $f_g = ||\vec{F}_g||/(cE_g)$ and the Eddington factor $\chi_g = P_g/E_g$, one can write:

$$ \mathcal{T}_g = \left( \frac{15}{\pi^4} \int_{\delta_g^{1/2}}^{\delta_g^{-1/2}} \int_{-1}^{1} \tilde{v}^3 [e^{\tilde{v}} \tilde{\alpha}_{0,g} (1+\beta_g \mu) - 1]^{-1} d\mu \, d\tilde{v} \right)^{1/4} , \quad (2.72) $$

$$ f_g = \frac{\int_{\delta_g^{1/2}}^{\delta_g^{-1/2}} \int_{-1}^{1} \tilde{v}^3 \mu [e^{\tilde{v}} \tilde{\alpha}_{0,g} (1+\beta_g \mu) - 1]^{-1} d\mu \, d\tilde{v}}{\int_{\delta_g^{1/2}}^{\delta_g^{-1/2}} \int_{-1}^{1} \tilde{v}^3 [e^{\tilde{v}} \tilde{\alpha}_{0,g} (1+\beta_g \mu) - 1]^{-1} d\mu \, d\tilde{v}} , \quad (2.73) $$

$$ \chi_g = \frac{\int_{\delta_g^{1/2}}^{\delta_g^{-1/2}} \int_{-1}^{1} \tilde{v}^3 \mu^2 [e^{\tilde{v}} \tilde{\alpha}_{0,g} (1+\beta_g \mu) - 1]^{-1} d\mu \, d\tilde{v}}{\int_{\delta_g^{1/2}}^{\delta_g^{-1/2}} \int_{-1}^{1} \tilde{v}^3 [e^{\tilde{v}} \tilde{\alpha}_{0,g} (1+\beta_g \mu) - 1]^{-1} d\mu \, d\tilde{v}} , \quad (2.74) $$

where $\mathcal{T}_g = k_B T_g / (h v_{m,g})$ is the dimensionless radiative temperature and $T_g = (E_g / a_R)^{1/4}$ is the radiative temperature of group g. As for any physically admissible pair $(E_g, ||\vec{F}_g||)$, there exists a unique pair of Lagrange multipliers $(\alpha_{0,g}, \beta_g)$. The equations above show that, for any pair $(\mathcal{T}_g, f_g)$, there exists a unique pair of Lagrange multipliers $(\alpha_{0,g}, \beta_g)$, which then allows the Eddington factor to be computed. One therefore concludes that, in the present case, the Eddington factor depends only on three parameters:
<figure><img src="image_12.png" /><figcaption>Figure 2.4: Evolution of the Eddington factor $\chi_g$ at different radiative temperatures $T_g$, using a group narrowness $\delta_g = 10^{-4}$.</figcaption></figure>

1. The dimensionless radiative temperature $T_g$: comparison between the equivalent radiative energy of the group, $k_B T_g$, and the geometric mean of the photon energy levels within the group, $h\nu_{m,g}$;

2. The reduced flux $f_g$: degree of anisotropy of the radiation within group $g$. Its values lie in [0; 1] ($f_g = 0$: isotropic radiation, $f_g \to 1$: highly directional radiation);

3. The group narrowness $\delta_g$: indicator of the separation between the frequency bounds of the group, $\nu_{g-1/2}$ and $\nu_{g+1/2}$, defined by the relation $\delta_g = \nu_{g-1/2}/\nu_{g+1/2}$, with values in the interval $[0; 1]$ ($\delta_g \to 0$: widely separated frequency bounds, $\delta_g \to 1$: very close frequency bounds).

In order to gain insight into the dependence of the Eddington factor on the reduced flux, I plotted its evolution for different values of the dimensionless radiative temperature in figure 2.4. Two common features can be observed across all curves:

1. For any value of $T_g$, $\chi_g$ always varies from 1/3 at $f_g = 0$ to 1 at $f_g = 1$, and has a zero slope at $f_g = 0$;

2. $T_g$ influences the evolution of $\chi_g$ only at intermediate values of $f_g$.

Overall, the Eddington factor depends more strongly on the reduced flux $f_g$ than on the radiative temperature and the group narrowness. However, this dependence does exist, and neglecting it leads to errors in the estimation of the Eddington factor of up to 16% for values
of $f_g = 0.65$, as shown in figure 2.4. Let us now analyze in more detail the dependence on $\mathcal{T}_g$ and $\delta_g$.

The ratio of the Eddington factor $\chi_g$ in the M1-multigroup model to its gray expression is shown as a function of the radiative temperature $\mathcal{T}_g$ for different values of the group narowness $\delta_g$ (see figure 2.5). Qualitatively, when the group is wide ($\delta_g \ll 1$), the presence of a central “plateau” can be observed within domain $I_b$ (see figure 2.5a). This plateau lies approximately in the interval $\mathcal{T}_g \in [\sqrt{\delta_g}; 1/\sqrt{\delta_g}]$, which corresponds to radiative temperatures such that $k_B T_g \in [h\nu_{g-1/2}; h\nu_{g+1/2}]$. As the group becomes narrower, i.e. as $\delta_g$ increases, the size of the “plateau” decreases until it disappears (see figures 2.5b and 2.5c). When the frequencies become nearly identical, two distinct regimes emerge: one at low temperature and one at high temperature. The transition between these regimes occurs at lower radiative temperatures as the frequencies approach one another ($\delta_g \to 1$, see figure 2.5d). Overall, three main regimes can be identified in figure 2.5:

1. **Domain L:** low radiative temperature. In this domain, $\mathcal{T}_g \ll \sqrt{\delta_g}$. $\chi_g$ depends only on $f_g$. This domain corresponds to that described by Minerbo [89];

2. **Domain I:** intermediate radiative temperatures. $\chi_g$ depends strongly on $\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}$. The value of $\chi_g$ depends only on $f_g$.

The intermediate domain I, based on the group narrowness, can be further subdivided into three cases:

1. **Case $\delta_g \to 0$:** wide group. Three intermediate domains can be identified (see figure 2.5a):

- **I_a:** $\chi_g$ depends strongly on $\mathcal{T}_g$ and $f_g$;

- **I_b:** corresponds to the domain in which the Eddington factor of the M1-gray model is valid. This domain is located approximately in the interval $\mathcal{T}_g \in [\sqrt{\delta_g}; 1/\sqrt{\delta_g}]$, and its width is therefore linked to $\delta_g$, i.e. it decreases as $\delta_g$ increases;

- **I_c:** $\chi_g$ depends strongly on $\mathcal{T}_g$ and $f_g$.

2. **Intermediate $\delta_g$:** in domains $I_d$ and $I_e$, $\chi_g$ depends on $\mathcal{T}_g$, $f_g$, and $\delta_g$ (see figures 2.5b and 2.5c);

3. **Case $\delta_g \to 1$:** narrow group. In domain $I_f$, $\chi_g$ is observed to depend only on $\mathcal{T}_g$ and $f_g$.
<figure><img src="image_13.png" /></figure>

(a) $\delta_g = 10^{-4}$. L, H: asymptotic domains; $I_a$: low-to-intermediate radiative temperatures; $I_b$: domain in which $\chi_R$ is valid; $I_c$: high-intermediate radiative temperatures.

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

(c) $\delta_g = 0.3$. L, H: asymptotic domains; $I_d$: low-to-intermediate radiative temperatures; $I_e$: intermediate-to-high radiative temperatures.

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

(b) $\delta_g = 0.07$. L, H: asymptotic domains; $I_d$: low-to-intermediate radiative temperatures; $I_e$: intermediate-to-high radiative temperatures.

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

(d) $\delta_g = 0.99$. L, H: asymptotic domains; $I_f$: intermediate radiative temperatures.

Figure 2.5: Ratio between the Eddington factor of the M1-multigroup model, $\chi_g$, and that of the M1-gray model, $\chi_R$, for a reduced flux $f_g = 0.65$, at different values of the group finesse $\delta_g$.

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

Figure 2.6: Evolution of the ratio between the Eddington factor of the M1-multigroup model, $\chi_g$, and that of the M1-gray model, $\chi_R$, as a function of $\delta_g$, illustrating the shift of the curve as $\delta_g$ approaches 1.
In the case $\delta_g \to 1$, figure 2.6 shows the curve of the Eddington factor $\chi_g$ as a function of the dimensionless radiative temperature $T_g$ for values of the group narrowness $\delta_g$ close to 1. The shape of the curve remains unchanged but is shifted toward lower radiative temperatures. In order to understand this behavior, one can seek an approximate expression for the Eddington factor at very high radiative temperatures $T_g$. Let us assume that the dimensionless Lagrange multiplier $\overline{\alpha}_{0,g}$ satisfies $\overline{\alpha}_{0,g} \ll \sqrt{\delta_g}$ (i.e. $\alpha_{0,g} \ll k_B/(h\nu_{g+1/2})$). The specific intensity in group g can then be written as:

$$I_{\nu,g}(\mu) \approx \frac{2k_B}{c^2} \frac{\nu^2}{\alpha_{0,g}(1 + \beta_g \mu)}, \qquad (2.75)$$

which allows the following expressions for the radiative temperature, the reduced flux, and the Eddington factor to be derived:

$$T_g = \left( \frac{5h^3 \frac{\nu_{g+1/2}^3 - \nu_{g-1/2}^3}{\alpha_{0,g}} \arctanh(\beta_g)}{\pi^4 k_B^3 \beta_g} \right)^{1/4}, \qquad (2.76)$$

$$f_g = - \frac{\arctanh(\beta_g) - \beta_g}{\beta_g \arctanh(\beta_g)}, \qquad (2.77)$$

$$\chi_g = \frac{\arctanh(\beta_g) - \beta_g}{\beta_g^2 \arctanh(\beta_g)}. \qquad (2.78)$$

Equations (2.77) and (2.78) show that, in this case, the Eddington factor depends only on the reduced flux $f_g$ and makes it possible to recover the correct expression for domain H. Since $\alpha_{0,g} \ll k_B/(h\nu_{g+1/2})$, one can derive a condition on $T_g$ that guarantees the validity of this approximation:

$$T_g \gg \frac{h\nu_{g+1/2}}{\pi k_B} \left(5(1 - \delta_g^3) \frac{\arctanh(\beta_g)}{\beta_g}\right)^{1/4}.$$

Since $\arctanh(\beta_g)/\beta_g \ge 1$ for $\beta_g \in ]-1; 1[$, this inequality can be rewritten in terms of $T_g$ as:

$$T_g \gg \frac{(1 - \delta_g^3)^{1/4}}{\sqrt{\delta_g}}. \qquad (2.79)$$

This inequality shows that the expression of the Eddington factor in domain H extends to progressively lower radiative temperatures. This explains the behavior observed in figure 2.6.

These observations apply in the general case where $\nu_{g-1/2} > 0$ and $\nu_{g+1/2} < +\infty$. However, if $\nu_{g-1/2} = 0$ or $\nu_{g+1/2} = +\infty$, the group narrowness $\delta_g$ can no longer be defined, nor can the dimensionless radiative temperature $T_g$. Let us now examine what happens in these two cases.

❖ Case $\nu_{g-1/2} = 0$ and $\nu_{g+1/2} < \infty$
<figure><img src="image_18.png" /><figcaption>Figure 2.7: Ratio between the Eddington factor of the M1-multigroup model, $\chi_g$, and that of the M1-gray model, $\chi_R$, for $f_g = 0.65$ and $\nu_{g-1/2} = 0$. $I_b$: domain in which $\chi_R$ is valid; $I_c$: intermediate-to-high radiative temperatures; H: asymptotically high radiative temperatures.</figcaption></figure>

In the case where $\nu_{g-1/2} = 0$, one can perform the change of variable $\tilde{\nu} = \nu/\nu_{g+1/2}$, which allows the radiative temperature, the reduced flux, and the Eddington factor to be expressed as:

$$ \frac{k_B T_g}{h \nu_{g+1/2}} = \left( \frac{15}{\pi^4} \int_0^1 \int_{-1}^1 \tilde{\nu}^3 [e^{\tilde{\nu}} \overline{\alpha_{0,g}} (1+\beta_g \mu) - 1]^{-1} d\mu \, d\tilde{\nu} \right)^{1/4}, \quad (2.80) $$

$$ f_g = \frac{\int_0^1 \int_{-1}^1 \tilde{\nu}^3 \mu [e^{\tilde{\nu}} \overline{\alpha_{0,g}} (1+\beta_g \mu) - 1]^{-1} d\mu \, d\tilde{\nu}}{\int_0^1 \int_{-1}^1 \tilde{\nu}^3 [e^{\tilde{\nu}} \overline{\alpha_{0,g}} (1+\beta_g \mu) - 1]^{-1} d\mu \, d\tilde{\nu}}, \quad (2.81) $$

$$ \chi_g = \frac{\int_0^1 \int_{-1}^1 \tilde{\nu}^3 \mu^2 [e^{\tilde{\nu}} \overline{\alpha_{0,g}} (1+\beta_g \mu) - 1]^{-1} d\mu \, d\tilde{\nu}}{\int_0^1 \int_{-1}^1 \tilde{\nu}^3 [e^{\tilde{\nu}} \overline{\alpha_{0,g}} (1+\beta_g \mu) - 1]^{-1} d\mu \, d\tilde{\nu}}, \quad (2.82) $$

where $\overline{\alpha_{0,g}} = h\nu_{g+1/2}\alpha_{0,g}/k_B$ is the dimensionless Lagrange multiplier, defined using the frequency $\nu_{g+1/2}$. Using these expressions, one can introduce the dimensionless radiative temperature $\mathcal{T}_g^{\nu_{g+1/2}} = k_B T_g / (h\nu_{g+1/2})$, which must be used in this case. This temperature can also be defined in the general case where $\nu_{g-1/2} > 0$ and $\nu_{g+1/2} < +\infty$, and it is related to the radiative temperature $\mathcal{T}_g$ by the relation $\mathcal{T}_g^{\nu_{g+1/2}} = \mathcal{T}_g \sqrt{\delta_g}$.

In this context, one can deduce that the Eddington factor depends only on two variables: the reduced flux $f_g$ and the dimensionless radiative temperature $\mathcal{T}_g^{\nu_{g+1/2}}$. The evolution of the Eddington factor as a function of the radiative temperature $\mathcal{T}_g^{\nu_{g+1/2}}$ is shown in figure 2.7. In this figure, one can identify domains $I_b$, $I_c$, and H, previously observed in the case of a wide group ($\delta_g \ll 1$). The only difference in the computation of the Eddington factor in this case lies in the choice of the dimensionless radiative temperature used.
<figure><img src="image_19.png" /><figcaption>Figure 2.8: Ratio between the Eddington factor of the M1-multigroup model, $\chi_g$, and that of the M1-gray model, $\chi_R$, for $f_g = 0.65$ and $\nu_{g+1/2} \to \infty$. L: asymptotically low radiative temperatures; $I_a$: low-to-intermediate radiative temperatures; $I_b$: domain in which $\chi_R$ is valid.</figcaption></figure>

### ❖ Case $\nu_{g-1/2} > 0$ and $\nu_{g+1/2} = \infty$

In the case where $\nu_{g+1/2} = +\infty$, one can instead perform the change of variable $\tilde{\nu} = \nu/\nu_{g-1/2}$, which allows the radiative temperature, the reduced flux, and the Eddington factor to be expressed as:

$$ \frac{k_B T_g}{h \nu_{g-1/2}} = \left( \frac{15}{\pi^4} \int_1^{+\infty} \int_{-1}^1 \tilde{\nu}^3 \left[ e^{\tilde{\nu} \widehat{\alpha_{0,g}} (1+\beta_g \mu)} - 1 \right]^{-1} d\mu \, d\tilde{\nu} \right)^{1/4}, \quad (2.83) $$

$$ f_g = \frac{\int_1^{+\infty} \int_{-1}^1 \tilde{\nu}^3 \mu \left[ e^{\tilde{\nu} \widehat{\alpha_{0,g}} (1+\beta_g \mu)} - 1 \right]^{-1} d\mu \, d\tilde{\nu}}{\int_1^{+\infty} \int_{-1}^1 \tilde{\nu}^3 \left[ e^{\tilde{\nu} \widehat{\alpha_{0,g}} (1+\beta_g \mu)} - 1 \right]^{-1} d\mu \, d\tilde{\nu}}, \quad (2.84) $$

$$ \chi_g = \frac{\int_1^{+\infty} \int_{-1}^1 \tilde{\nu}^3 \mu^2 \left[ e^{\tilde{\nu} \widehat{\alpha_{0,g}} (1+\beta_g \mu)} - 1 \right]^{-1} d\mu \, d\tilde{\nu}}{\int_1^{+\infty} \int_{-1}^1 \tilde{\nu}^3 \left[ e^{\tilde{\nu} \widehat{\alpha_{0,g}} (1+\beta_g \mu)} - 1 \right]^{-1} d\mu \, d\tilde{\nu}}, \quad (2.85) $$

where $\widehat{\alpha_{0,g}} = h\nu_{g-1/2}\alpha_{0,g}/k_B$ is the dimensionless Lagrange multiplier, defined using the frequency $\nu_{g-1/2}$. From these expressions, one can introduce the dimensionless radiative temperature $T_g^{\nu_{g-1/2}} = k_B T_g / (h\nu_{g-1/2})$, which must be used in this case. This temperature can also be defined in the general case where $\nu_{g-1/2} > 0$ and $\nu_{g+1/2} < +\infty$, and it is related to the radiative temperature $T_g$ by the relation $T_g^{\nu_{g-1/2}} = T_g / \sqrt{\delta_g}$.

In this case, one can deduce that the Eddington factor depends only on two variables: the reduced flux $f_g$ and the radiative temperature $T_g^{\nu_{g-1/2}}$. The evolution of the Eddington factor as a function of the radiative temperature $T_g^{\nu_{g-1/2}}$ is shown in figure 2.8. In this figure, one finds
domains L, $I_a$, and $I_b$, previously observed in the case of a wide group ($\delta_g \ll 1$). The only
difference lies in the choice of the radiative temperature used.

## 2.3 Coupling

Up to this point, the radiative model has been formulated within the framework of a fluid at rest. However, in radiation hydrodynamics models, light propagates through a moving fluid, which induces Doppler effects that modify the energy and momentum exchanged between the fluid and the radiation. Two choices of reference frame are then possible to describe the radiative quantities:

1. **The comoving frame**, in which the fluid is locally at rest. In this framework, matter–radiation interactions (emission and absorption) are isotropic, so that the source terms of the radiative transfer equations (2.47) remain unchanged. However, transforming the radiative quantities from the comoving frame to the laboratory frame introduces additional terms into the equations. In particular, advection effects appear in the form of terms such as $\vec{\nabla} \cdot (\vec{v}E_g)$ for the radiative energy and $\vec{\nabla} \cdot (\vec{v} \otimes \vec{F}_g)$ for the radiative flux [95, 99, 100]. Moreover, the Doppler effect induces a change in photon frequency between the two frames. Within the M1-multigroup model, this implies transfers of energy and momentum between spectral groups, represented by the terms $\vec{\nabla}\vec{v} : \int_{v_{g-1/2}}^{v_{g+1/2}} \partial_v (vPv) dv$ and $\vec{\nabla}\vec{v} : \int_{v_{g-1/2}}^{v_{g+1/2}} \partial_v (vQ_v) dv$. These terms make the solution of the radiative transfer equations in the M1-multigroup model particularly complex [35];

2. **The laboratory frame**, in which the fluid is in motion. In this frame, the radiative transfer equations (2.47) retain their usual form, without the introduction of advection terms or energy transfers between spectral groups in the M1-multigroup model. By contrast, the source terms become more complex, as they must include Doppler effects and geometric aberrations [101].

In this work, we adopt this second approach by expressing the radiative quantities in the
laboratory frame.

### 2.3.1 Lorentz transformation

Before proceeding, let us recall the expression of the Lorentz transformation. Consider two
reference frames, $\mathcal{R}$ and $\mathcal{R}'$, where $\mathcal{R}'$ moves with velocity $\vec{v}$ relative to $\mathcal{R}$. In special relativity,
the transformation of a four-vector $X^\alpha$ from the reference frame $\mathcal{R}$ to the reference frame $\mathcal{R}'$ is expressed as follows:

$$X'^\beta = \Lambda^\beta_\alpha X^\alpha, \qquad (2.86)$$

where $X'^\beta$ is the four-vector in the reference frame $\mathcal{R}'$, and $\Lambda^\beta_\alpha$ is the Lorentz tensor, given by:

$$\Lambda^\beta_\alpha = \begin{bmatrix} \gamma_u & -\gamma_u \vec{\beta}_u^T \\ -\gamma_u \vec{\beta}_u & \mathbb{I}_d + \frac{\gamma_u^2}{\gamma_u+1} \vec{\beta}_u \otimes \vec{\beta}_u \end{bmatrix}, \qquad (2.87)$$

with:

$$\gamma_u = \frac{1}{\sqrt{1 - \beta_u^2}}, \qquad \vec{\beta}_u = \vec{v}/c, \qquad \beta_u = ||\vec{\beta}_u||.$$

The quantity $\gamma_u$ is called the *Lorentz factor*. Let us first consider the case of the M1-multigroup model. We then define the second-order radiation energy–momentum tensor of group g, denoted $\mathbb{R}_g = (R_g^{\alpha\beta})_{\alpha,\beta}$, as:

$$R_g^{\alpha\beta} = c^{-1} \int_{\nu_{g-1/2}}^{\nu_{g+1/2}} \int_{4\pi} I_v(\vec{n}) n^\alpha n^\beta \, d\Omega \, dv, \qquad (2.88)$$

where $n^\alpha = (n^0, \vec{n}) = (1, \vec{n})$. This tensor can also be written in the form:

$$R_g^{\alpha\beta} = \begin{bmatrix} E_g & \vec{F}_g^T/c \\ \vec{F}_g/c & \mathbb{P}_g \end{bmatrix}. \qquad (2.89)$$

To transform this tensor from the comoving frame to the laboratory frame, it is necessary to apply the Lorentz transformation to each contravariant index. Knowing that the laboratory frame moves with velocity $\vec{v}$ relative to the comoving frame, the radiative quantities in the latter, namely $E_g$, $\vec{F}_g$, and $\mathbb{P}_g$, can be expressed as functions of the radiative quantities in the laboratory frame according to the following relations:

$$E_{g,0} = \gamma_u^2 \left[ E_g - \frac{2\vec{v} \cdot \vec{F}_g}{c^2} + \vec{\beta}_u^T \mathbb{P}_g \vec{\beta}_u \right], \qquad (2.90)$$

$$\vec{F}_{g,0} = \gamma_u \left[ \vec{F}_g - \left\{ \mathbb{P}_g + \gamma_u \mathbb{I}_d \left( E_g - \frac{2\gamma_u + 1}{\gamma_u + 1} \frac{\vec{v} \cdot \vec{F}_g}{c^2} + \frac{\gamma_u \vec{\beta}_u^T \mathbb{P}_g \vec{\beta}_u}{\gamma_u + 1} \right) \right\} \vec{v} \right], \qquad (2.91)$$
$$
\begin{equation}
\begin{split}
\mathbb{P}_{g,0} = \mathbb{P}_g - \gamma_u \frac{\vec{v} \otimes \vec{F}_g + \vec{F}_g \otimes \vec{v}}{c^2} - \gamma_u^2 \left[ \frac{\mathbb{P}_g \vec{\beta}_u \otimes \vec{\beta}_u + \vec{\beta}_u \otimes \vec{\beta}_u \mathbb{P}_g}{\gamma_u + 1} + \right. \\
\left. \qquad \vec{\beta}_u \otimes \vec{\beta}_u \left\{ \mathbb{E}_g - \frac{2\gamma_u}{\gamma_u + 1} \left( \frac{\vec{v} \cdot \vec{F}_g}{c^2} - \frac{\gamma_u}{2} \vec{\beta}_u^T \mathbb{P}_g \vec{\beta}_u \right) \right\} \right] .
\end{split}
\tag{2.92}
\end{equation}
$$

Let us now define the four-vector of source terms of the radiative transfer equation as $S_g^\alpha = (S_g^0, \vec{S}_g)$. By applying the Lorentz transformation, the source terms in the laboratory frame can be expressed in terms of those in the comoving frame through the following relations:

$$
S_g^0 = \gamma_u (S_{g,0}^0 + \vec{\beta}_u \cdot \vec{S}_{g,0}) , \qquad (2.93)
$$

$$
\vec{S}_g = \vec{S}_{g,0} + \gamma_u \vec{\beta}_u \left( S_{g,0}^0 + \frac{\gamma_u \vec{\beta}_u \cdot \vec{S}_{g,0}}{\gamma_u + 1} \right) . \quad (2.94)
$$

As we saw in Section 2.2. Radiative transfer, in the comoving frame the source term can
be written as:

$$
\begin{align*}
S_{g,0}^{0} &= \kappa_{P,g} (\mathbf{E}_{g,0} - a_{R} \theta_{g}^{4}(\mathbf{T})) , \\
\vec{S}_{g,0} &= \kappa_{R,g} \vec{F}_{g,0} / c ,
\end{align*}
$$

where $\kappa_{P,g}$ and $\kappa_{R,g}$ are the Planck and Rosseland opacities of group g in the comoving frame.
Using transformations (2.93) and (2.94), the source terms can therefore be rewritten in the
laboratory frame as:

$$
\begin{align*}
S_g^0 &= \gamma_u \left( \kappa_{P,g} (E_{g,0} - a_R \theta_g^4(T)) + \kappa_{R,g} \frac{\vec{v} \cdot \vec{F}_{g,0}}{c^2} \right) , \\
\vec{S}_g &= \kappa_R \vec{F}_{g,0} / c + \gamma_u \vec{\beta}_u \left( \kappa_{P,g} (E_{g,0} - a_R \theta_g^4(T)) + \frac{\gamma_u \kappa_{R,g}}{\gamma_u + 1} \frac{\vec{v} \cdot \vec{F}_{g,0}}{c^2} \right) .
\end{align*}
$$

Finally, by rewriting this source term in terms of the radiative quantities expressed in the
laboratory frame and using transformations (2.90), (2.91), and (2.92), one obtains the following
expression for the source terms in the laboratory frame [101]:

$$
S_g^0 = \gamma_u \left[ \begin{aligned}
& \left( \gamma_u^2 (\kappa_{P,g} - \beta_u^2 \kappa_{R,g}) \mathbf{E}_{g,0} - \kappa_{P,g} a_R \theta_g^4(\mathbf{T}) \right) + \\
& \qquad \gamma_u^2 \left\{ \left( (1+\beta_u^2) \kappa_{R,g} - 2\kappa_{P,g} \right) \vec{v} \cdot (c^{-2}\vec{F}_g) - (\kappa_{R,g} - \kappa_{P,g}) \vec{\beta}_u^T \mathbb{P}_g \vec{\beta}_u \right\}
\end{aligned} \right] , \quad (2.95)
$$

$$
\begin{equation}
\begin{split}
\vec{S}_g ={}& \gamma_u \left[ \kappa_{R,g} \vec{F}_g / c - \kappa_{R,g} \mathbb{P}_g \vec{\beta}_u - \vec{\beta}_u \left\{ \gamma_u^2 (\kappa_{R,g} - \kappa_{P,g}) \mathbb{E}_g + \kappa_{P,g} a_R \theta_g^4 (\mathbf{T}) \right. \right. \\
& \qquad \left. \left. - \gamma_u^2 (\kappa_{R,g} - \kappa_{P,g}) (2 \vec{v} \cdot (c^{-2} \vec{F}_g) - \vec{\beta}_u^T \mathbb{P}_g \vec{\beta}_u) \right\} \right] .
\end{split}
\tag{2.96}
\end{equation}
$$
In this work, we assume that the fluid is non-relativistic ($\beta_u \ll 1$). Therefore, these source terms can be expanded to first order in $\beta$, noting that $\gamma_u = 1 + O(\beta_u^2)$:

$$S_g^0 = \kappa_{P,g} (E_g - a_R \theta_g^4(T)) + (\kappa_{R,g} - 2\kappa_{P,g}) \vec{v} \cdot (c^{-2}\vec{F}_g) , \quad (2.97)$$

$$\vec{S}_g = \kappa_{R,g} \vec{F}_g / c - \kappa_{R,g} \mathbb{P}_g \vec{\beta}_u - \vec{\beta}_u \{ (\kappa_{R,g} - \kappa_{P,g}) E_g + \kappa_{P,g} a_R \theta_g^4(T) \} . \quad (2.98)$$

## 2.3.2 Radiative hydrodynamics equations

We can now combine the hydrodynamics equations (2.8) and the radiative transfer equations into a coupled system, forming what is known as the *radiation hydrodynamics equations*. These equations describe the evolution of a fluid interacting with radiation.

Let us first present these equations within the framework of the M1-multigroup model for the description of radiation. In this context, the source terms $cS_g^0$ and $\vec{S}_g$ introduced previously (equations (2.97) and (2.98)) represent, respectively, the exchanges of energy and momentum between the fluid and the radiation contained in each frequency group. These exchanges are taken into account by adding the total contribution of these interactions as source terms on the right-hand side of the hydrodynamics equations (2.8) and by incorporating them into the radiative transfer equations (2.47), formulated in the laboratory frame. One thus obtains the complete system of radiation hydrodynamics equations:

$$\left\{
\begin{aligned}
\partial_t \rho + \vec{\nabla} \cdot (\rho \vec{v}) &= 0, \\
\partial_t (\rho \vec{v}) + \vec{\nabla} \cdot (\rho \vec{v} \otimes \vec{v}) + p \mathbb{I}_d &= \sum_{g=1}^{\mathcal{G}} \vec{S}_g, \\
\partial_t E + \vec{\nabla} \cdot ((E + p)\vec{v}) &= \sum_{g=1}^{\mathcal{G}} c S_g^0, \\
\partial_t E_g + \vec{\nabla} \cdot \vec{F}_g &= -c S_g^0, && \forall g \in [1, \mathcal{G}], \\
\partial_t (c^{-2} \vec{F}_g) + \vec{\nabla} \cdot \mathbb{P}_g &= -\vec{S}_g, && \forall g \in [1, \mathcal{G}].
\end{aligned}
\right.
\quad (2.99)$$

In the particular case of the M1-gray model, the radiation hydrodynamics equations take a simpler form:

$$\left\{
\begin{aligned}
\partial_t \rho + \vec{\nabla} \cdot (\rho \vec{v}) &= 0, \\
\partial_t (\rho \vec{v}) + \vec{\nabla} \cdot (\rho \vec{v} \otimes \vec{v}) + p \mathbb{I}_d &= \vec{S}, \\
\partial_t E + \vec{\nabla} \cdot ((E + p)\vec{v}) &= c S^0, \\
\partial_t E_R + \vec{\nabla} \cdot \vec{F}_R &= -c S^0, \\
\partial_t (c^{-2} \vec{F}_R) + \vec{\nabla} \cdot \mathbb{P}_R &= -\vec{S},
\end{aligned}
\right.
\quad (2.100)$$

where $S^0$ and $\vec{S}$ are the source terms of the single group considered, integrated over the entire electromagnetic spectrum.
### 2.3.3 Asymptotic cases

Now that a general radiation hydrodynamics model has been established, it is essential to verify whether it is able to recover the well-known asymptotic limits: the optically thin regime and the optically thick regime. To this end, one can first identify several characteristic lengths defined from the source terms (2.97) and (2.98):

$$ \lambda_{th} = 1/\max_{v}(\kappa_{v}) , $$

$$ \lambda_{tr} = 1/\max_{v}(\chi_{v}) , $$

$$ \lambda_{dyn,1} = 1/\max_{v}(|\chi_{v} - 2\kappa_{v}|) , $$

$$ \lambda_{dyn,2} = 1/\max_{v}(\sigma_{v}) . $$

These characteristic lengths, in turn, allow four associated characteristic times to be defined:

$$ \tau_{th} = \frac{\lambda_{th}}{c} , \qquad (\text{temps de thermalisation}) $$

$$ \tau_{tr} = \frac{\lambda_{tr}}{c} , \qquad (\text{temps de transport}) $$

$$ \tau_{dyn,1} = \frac{\lambda_{dyn,1}}{\beta_{u}c} , \qquad (\text{temps dynamique 1}) $$

$$ \tau_{dyn,2} = \frac{\lambda_{dyn,2}}{\beta_{u}c} . \qquad (\text{temps dynamique 2}) $$

These characteristic times play a key role in the dynamics of radiation-matter interactions:

* $\tau_{th}$ corresponds to the time required for the gas to reach thermal equilibrium with the radiation,

* $\tau_{tr}$ represents the characteristic timescale of radiative transport,

* the times $\tau_{dyn,1}$ and $\tau_{dyn,2}$ quantify the dynamical interaction between radiation and matter in the presence of nonzero fluid velocities.

Three important observations can be made:

#### 1. Hierarchy of timescales:

In the non-relativistic framework ($\beta_{u} \ll 1$), one observes that $\tau_{th}$, $\tau_{tr}$ $\ll \tau_{dyn,1}$, $\tau_{dyn,2}$.

This means that radiation has time to thermalize and to be transported before dynamical
effects associated with fluid motion (Doppler effect, relativistic distortions) become significant. Thus, on short timescales, radiation behaves as if it were in a rest frame;

## 2. Contribution of higher-order terms in $β_u$:

If the expansion of the source terms $S_0$ and $\vec{S}_g$ had not been truncated at first order in $β_u$, additional contributions would appear, associated with characteristic times of the form $τ_n = λ/(β_u^n c)$ with $n ≥ 2$. However, as long as one considers situations in which the characteristic durations are sufficiently short, the impact of these terms remains negligible on the global dynamics of the fluid and the radiation;

## 3. Case where scattering is neglected ($σ_ν ≈ 0$):

By neglecting photon scattering and assuming that the absorption coefficient is independent of photon frequency ($κ_ν ≈ χ_ν = κ$), the characteristic times simplify:

$$τ_{tr} ≈ τ_{th} = τ = \frac{1}{κC}, \quad τ_{dyn,1} ≈ \frac{τ}{β_u}, \quad τ_{dyn,2} → +∞.$$

In this case, the behavior of the system depends on the ratio between the characteristic time $t$ of the situation under consideration and the timescales defined above [102]:

### a) $t << τ$: Optically thin regime (free streaming)

Radiation propagates freely with negligible interaction with matter. This regime corresponds to the optically thin case;

### b) $t >> τ$ and $t << τβ_u$: Static diffusion regime

Radiation is in thermal equilibrium with matter, but its dynamical interaction with the fluid remains negligible;

### c) $t >> τ$, $t >> τβ_u$, and $t << τβ_u^2$: Dynamic diffusion regime

The interaction between radiation and gas becomes dominant, and radiative transport is essentially governed by fluid advection. This regime remains within the validity limits of the model developed here.

Thus, depending on the characteristic timescale considered, the physical situation differs, and different approximations may be employed. Let us now present the two main existing approximations used to describe radiation.

### ❖ Optically thin case
In the interstellar medium, matter is very often optically thin because of the low densities and temperatures that prevail there. This implies weak interactions between radiation and matter, which translates into small opacities $κ_P$ and $κ_R$. Consequently, the characteristic thermalization and transport times, denoted $τ_{th}$ and $τ_{tr}$, respectively, are very long. In other words, radiation takes a long time to exchange energy with the surrounding matter and to propagate through the medium. It is therefore reasonable to assume that, in this regime, the influence of radiation on the fluid dynamics remains negligible.

Let us now consider a situation in which, initially, radiation is in radiative equilibrium with matter, which implies that the radiative flux is nearly zero ($\vec{F}_R \approx 0$) and that the radiative energy satisfies $E_R = a_R T_i^4$, where $T_i$ is the initial temperature of the fluid. Suppose that, at some time, the fluid is heated, reaching a temperature $T_f$ such that $T_f^4 \gg T_i^4$. In this case, the source terms describing matter–radiation interactions, given by expressions (2.97) and (2.98), can be approximated by:

$$S^0 \approx -\kappa_P(\rho, T)a_R T^4 = -\Lambda(\rho, T), \qquad (2.101)$$

$$\vec{S} \approx 0, \qquad (2.102)$$

where $κ_P$ is the Planck mean opacity, defined by equation (2.38), taking $ν_{g-1/2} = 0$ and $ν_{g+1/2} = +\infty$. One recognizes here the typical form of the cooling function $\Lambda(\rho, T)$, which depends only on the gas density $\rho$ and its temperature $T$. Assuming that the Planck opacity can be described by a power law in $\rho$ and $T$, this cooling function can be written in the general form:

$$\Lambda(\rho, T) \approx \Lambda_0 \rho^\alpha T^\beta, \qquad (2.103)$$

where $Λ_0$ is a constant that depends in particular on the radiation constant $a_R$, while $\alpha$ and $\beta$ are exponents characterizing the physical processes involved. For instance, to model bremsstrahlung cooling in a hot interstellar medium ($T > 10^7$ K), one would choose $\alpha = 2$ and $\beta = 0.5$. By contrast, for cooling via inverse Compton scattering of electrons, the commonly used values are $\alpha = 1$ and $\beta = 1$.

Since we assume that the gas follows the ideal-gas equation of state (2.10), it is also possible to express this cooling function in terms of the thermal pressure $p$ and the density $\rho$ as follows:

$$\Lambda(\rho, p) \approx \Lambda'_0 \rho^\epsilon p^\zeta, \qquad (2.104)$$

where $Λ'_0$ is a constant related to $Λ_0$ through $Λ'_0 = Λ_0μm_H/k_B$. The exponents $\epsilon$ and $\zeta$ are related to $\alpha$ and $\beta$ by $\epsilon = \alpha - \beta$ and $\zeta = \beta$. In this framework, the simplified radiation
hydrodynamics equations take the form:

$$ 
\left\{ 
\begin{aligned}
\partial_t \rho + \vec{\nabla} \cdot (\rho \vec{v}) &= 0, \\
\partial_t (\rho \vec{v}) + \vec{\nabla} \cdot (\rho \vec{v} \otimes \vec{v} + p) &= 0, \\
\partial_t E + \vec{\nabla} \cdot ((E + p) \vec{v}) &= -\Lambda(\rho, p).
\end{aligned}
\right. 
\qquad (2.105) $$

This formulation is widely used because of its numerical simplicity: by avoiding the explicit resolution of the radiative transfer equations, it significantly reduces both the complexity of the calculations and the computational cost.

### ❖ Optically thick case

Conversely, in dense regions such as stellar interiors, the medium is highly opaque and interacts strongly with radiation. In this regime, radiation is in near-equilibrium with matter, since the characteristic thermalization and transport times are very short compared with the dynamical timescales under consideration ($t \gg \tau_{th}, \tau_{tr}$). One can then assume that radiative equilibrium is reached instantaneously and express the radiative pressure as well as the radiative flux as follows:

$$ \mathbb{P}_{R,0} = \frac{1}{3} \mathbb{I}_d E_{R,0} , \qquad \vec{F}_{R,0} = - \frac{c}{3\kappa_R} \vec{\nabla} E_{R,0} , $$

where $\kappa_R$ is the Rosseland opacity defined by (2.39), with $\nu_{g-1/2} = 0$ and $\nu_{g+1/2} = +\infty$. These expressions reflect the fact that, in an optically thick medium, radiation interacts strongly with matter and propagates by diffusion. The radiative pressure is isotropic, while the radiative flux follows a diffusion law, proportional to the gradient of the radiative energy, with a diffusion coefficient given by $c/(3\kappa_R)$.

Using these approximations, it becomes possible to retain only the radiative energy equation to describe radiation. Keeping only terms of order $O(v/c)$ as well as certain terms of order $O((v/c)^2)$, and expressing the radiative quantities in the comoving frame, the radiation hydrodynamics equations can be rewritten in the form [42]:

$$ 
\left\{ 
\begin{aligned}
\partial_t \rho + \vec{\nabla} \cdot (\rho \vec{v}) &= 0, \\
\partial_t (\rho \vec{v}) + \vec{\nabla} \cdot (\rho \vec{v} \otimes \vec{v} + p) &= -\frac{1}{3} \vec{\nabla} E_{R,0}, \\
\partial_t E + \vec{\nabla} \cdot ((E + p) \vec{v}) &= c \kappa_P (E_{R,0} - a_R T^4) + \frac{1}{3} \vec{v} \cdot \vec{\nabla} E_{R,0}, \\
\partial_t E_{R,0} + \vec{\nabla} \cdot \left( \frac{4}{3} \vec{v} E_{R,0} \right) &= \vec{v} \cdot \left( \frac{c}{3\kappa_R} \vec{\nabla} E_{R,0} \right) - c \kappa_P (E_{R,0} - a_R T^4) - \frac{1}{3} \vec{v} \cdot \vec{\nabla} E_{R,0}.
\end{aligned}
\right. 
\qquad (2.106) $$
Compared with the complete system (2.100), this simplified model contains only four equations instead of five, which greatly facilitates its numerical solution. It thus avoids the explicit treatment of detailed radiative transfer, while still capturing the essential interactions between matter and radiation in an optically thick medium. This formalism describes both the static diffusion regime and the dynamic diffusion regime, as introduced by Krumholz [102]. However, it does not allow one to treat regions where radiative equilibrium is locally broken, as is the case in radiative shocks (see Section 4.1. *Introduction to radiative shocks*).

## 2.4 Synthesis

In this chapter, I have presented the model used to simulate radiation hydrodynamics. This model relies on several simplifying assumptions:

* the gas is ideal, inviscid, and without thermal diffusion,

* it is assumed to be non-ionized,

* gravitational effects are neglected,

* the gas and the radiation are in LTE,

* photon scattering is assumed to be elastic and isotropic, with an isotropic scattering coefficient,

* the fluid is treated within the non-relativistic framework.

These assumptions make it possible to derive the radiation hydrodynamics equations based on the M1-multigroup model for the description of radiative transfer. Since the Eddington factor associated with this model does not admit an analytical expression but can be determined numerically, I have studied its dependencies and variations in order to identify its main characteristics. The next chapter introduces AI and describes the work I have carried out to leverage it with the aim of improving the efficiency of radiation hydrodynamics simulations.
## CALCULATION OF THE EDDINGTON FACTOR WITH ARTIFICIAL INTELLIGENCE

The HADES code [37, 38, 39], developed to solve the radiation hydrodynamics equations, models the interaction between radiation and matter outside radiative equilibrium. It was designed to limit approximations in the treatment of the radiation field as much as possible. Two radiative transfer models are implemented in the code: the M1-gray model and the M1-multigroup model, which make it possible to solve the equations introduced in the previous chapter using classical spatial and temporal discretization methods. We first present the numerical methods used in this code and rewrite the radiation hydrodynamics equations in the following form:

$$ \begin{cases} \partial_t \mathbf{q}_h + \nabla \cdot \mathbf{F}(\mathbf{q}_h) = \mathbf{s}_h & (\text{hydrodynamics}) \\ \partial_t \mathbf{q}_r + \nabla \cdot \mathbf{F}(\mathbf{q}_r) = \mathbf{s}_r & (\text{radiative transfer}) \end{cases} $$

The algorithm used in HADES to solve the radiation hydrodynamics equations consists of two main steps [37]:

1. Update of the radiative quantities $\mathbf{q}_r$ over two time steps $\Delta t$,

2. Update of the hydrodynamic quantities $\mathbf{q}_h$ over two time steps $\Delta t$.

Each of these steps is further decomposed into three substeps:

1. Resolution of the Ordinary Differential Equation (ODE) $\partial_t \mathbf{q} = \mathbf{s}$ over a time step $\Delta t/2$ (explicit scheme for hydrodynamics, implicit scheme for radiation),
2. Resolution of the Partial Differential Equations (PDE) $\partial_t q + \nabla \cdot \mathbf{F}(q) = 0$ over a time step $\Delta t$ (for instance, an Harten-Lax-van Leer-Einfeld (HLLE) Riemann solver for hydrodynamics and an Harten-Lax-van Leer (HLL) Riemann solver for radiation),

3. Resolution of the ODE $\partial_t q = s$ over a time step $\Delta t/2$ (explicit scheme for hydrodynamics, implicit scheme for radiation).

This algorithm guarantees second-order accuracy in the solution of the radiation hydrodynamics equations. However, two major limitations remain:

1. **Computation of the Eddington factor in the M1-multigroup model**

In this model, the Eddington factor does not admit an explicit analytical expression. The code therefore relies on numerical root-finding algorithms to estimate it accurately. Although these methods provide precise values of the Eddington factor, they are particularly expensive in terms of computational cost;

2. **Courant–Friedrichs–Lewy (CFL) condition**

To ensure numerical stability when solving the radiation hydrodynamics equations, the time step $\Delta t$ and the spatial step $\Delta x$ must satisfy a CFL-type condition: $c\Delta t/\Delta x \le 1$, where c denotes the speed of light. This constraint forces the use of very small time steps, which significantly increases the computational cost of the simulations.

In order to overcome these limitations and accelerate the simulations, I explored the potential contribution of AI. This technology has attracted growing interest and is experiencing rapid development in fields as diverse as medicine, robotics, and space exploration. More recently, it has been integrated into numerical simulation workflows, giving rise to an emerging research field known as NeuHPC [67]. This field develops solutions along two main directions:

1. **Inverse modeling:** The goal is to discover the analytical or differential equations underlying the data. This approach includes symbolic regression<fnref n="1" /> [103], as well as the use of PINNs to infer the differential equations governing the observed phenomena;

2. **Simulation:** AI is used to solve known differential equations, either by fully integrating them using machine-learning techniques or by combining traditional numerical methods with AI-based approaches.

<footnote n="1">¹Python library: https://github.com/MilesCranmer/PySR, which aims at identifying explicit analytical expressions; Julia library: https://github.com/MilesCranmer/SymbolicRegression.jl</footnote>
In this work, I focus more specifically on the second direction, with the aim of developing more efficient methods to solve the radiation hydrodynamics equations. Three broad categories of neural-network applications in simulation can be identified:

1. **AI solver:** These methods consist in using AI to solve differential equations end to end. For instance, PINNs approximate the solution by treating the neural network as a function that approximates the solution of the underlying differential equations [68, 104, 105]. Another strategy relies on neural operators, which learn to map the solution at time $t$ to that at time $t + \Delta t$, using architectures such as Graph Neural Networks (GNNs) [106, 107, 108, 109] or Fourier Neural Operators (FNOs) [110, 111];

2. **Closure relations:** In this case, neural networks are used to determine or improve the closure relations of the equations. This approach can be applied within simplified models, for example for turbulence modeled by the RANS or LES equations [112, 113, 114], or to establish an unknown closure relation [115];

3. **Error correction:** This approach consists in training a neural network on data generated by very high-resolution simulations and then using it to correct the errors made in low-resolution simulations. This procedure makes it possible to drastically reduce computational cost while maintaining high accuracy in the results [116].

In this work, I explored two complementary approaches to reduce the computational cost of radiation hydrodynamics simulations. The first consists in developing an AI-based method to efficiently estimate the Eddington factor involved in the closure relation of the M1-multigroup model (see Section 3.2. *Eddington factor in the M1-multigroup model*). The second aims at investigating the feasibility of a neural solver, based on PINNs, capable of extrapolating a simulation from the initial data of an existing simulation, while relying directly on the radiation hydrodynamics equations (see Chapter 5. *Extrapolation of radiative shock simulations with Physics-Informed Neural Networks*).

In the present chapter, I will first present the main AI techniques used for training neural networks, before detailing the proposed method for accelerating the evaluation of the Eddington factor.

### 3.1 Steps to generate a functional neural network

When addressing a problem in AI, several fundamental questions must be posed and resolved during the development of the method.
* Which data should be used?

This involves determining the nature, quality, format, and quantity of the data required to train the network.

* How should the model error be measured (cost function)?

The choice of the loss function (evaluation of the neural network's performance) is essential to guide the learning process.

* Which optimization algorithm should be used?

The optimizer adjusts the parameters of the neural network in order to minimize the loss function.

* Which neural network architecture should be adopted?

A neural network can be viewed as a more or less complex function, characterized by a large number of adjustable parameters. The choice of its architecture—that is, the way in which the computational units are organized and connected—determines the form of this function. Among the common architectures in AI, one finds in particular Multi-Layer Perceptron (MLP)s, Generative Adversarial Networks (GAN)s, or Long Short-Term Memory (LSTM)s. This decision depends on the type of relationship to be modeled, the complexity of the problem, as well as constraints in terms of computational resources and available data.

* Should a particular strategy be adopted to promote generalization?

This involves helping the neural network produce accurate predictions on configurations it has not seen during training.

These choices are not independent and must be adjusted iteratively. An initial combination of assumptions is tested and then modified in light of the results obtained, until convergence toward a solution deemed satisfactory. This process is illustrated in figure 3.1 in the form of a decision loop. Indeed, answering one question (for example, *which data should be used?*) makes it possible to address the next one (for example, *which architecture should be adopted?*), and so on, up to the evaluation of the effectiveness of the choices made. This cycle is repeated several times with different combinations before retaining the most successful trial.

In this case, since the objective is to reduce the computational cost of numerical simulations, I imposed a strong constraint on the architecture: I opted for MLPs, one of the simplest architectures, yet sufficiently expressive to approximate complex nonlinear relationships while remaining computationally inexpensive.
<figure><img src="image_20.png" /><figcaption>Figure 3.1: Diagram illustrating the reasoning process involved in developing an artificial intelligence strategy.</figcaption></figure>

In this section, I will first introduce MLPs, before detailing the possible choices regarding the loss function, the optimization strategy, as well as techniques aimed at improving generalization capability.

### 3.1.1 Multi-Layer Perceptrons (MLP)

An MLP can be viewed as a complex mathematical function, having a large number of adjustable parameters depending on the problem being addressed. It is composed of elementary units called *perceptrons* (see figure 3.2), whose formulation was first introduced by Rosenblatt (1958) [117]. Each perceptron receives $n$ inputs $x_1, x_2, \dots, x_n$, associated with *weights* $w_1, w_2, \dots, w_n$, as well as a *bias* $b$. It produces an output $y$, defined by the following expression:

$$y = \sigma \left( \sum_{i=1}^{n} w_i x_i + b \right),$$

where $\sigma$ denotes an *activation function*, which may be linear or nonlinear. The introduction of nonlinearity allows the perceptron to capture complex behaviors and to model nonlinear phenomena. Several activation functions are commonly used:

* the identity function (see figure 3.3a)

$$\sigma(x) = x \in \mathbb{R},$$
<figure><img src="image_21.png" /><figcaption>Figure 3.2: Diagram of a perceptron</figcaption></figure>

* the sigmoid function (see figure 3.3b)

$$ \sigma(x) = \frac{1}{1 + e^{-x}} \in [0; 1] , $$

* the hyperbolic tangent function (see figure 3.3c)

$$ \sigma(x) = \tanh(x) = \frac{1 - e^{-x}}{1 + e^{-x}} \in [-1; 1] , $$

* the arctangent function (see figure 3.3d)

$$ \sigma(x) = \operatorname{Arctan}(x) \in [-\pi/2; \pi/2] , $$

* the Rectified Linear Unit (ReLU) (see figure 3.3e)

$$ \sigma(x) = \begin{cases} 0 & \text{if } x < 0 \\ x & \text{if } x \ge 0 \end{cases} \in \mathbb{R}^{+} , $$

* the soft rectified linear unit (SoftPlus) (see figure 3.3f)

$$ \sigma(x) = \ln(1 + e^{-x}) \in \mathbb{R}^{+} , $$

Although most of these activation functions are nonlinear, they generally do not allow the representation of discontinuities. Some alternative variants have nevertheless been proposed, notably by DellaSanta (2023) [118], which explicitly introduce discontinuities into the activation function in order to endow the network with specific properties. The adoption of such functions, however, requires particular care, both from a theoretical and a numerical standpoint, due to the consequences they may have on training stability and the local differentiability of the model.
<figure><img src="image_22.png" /><figcaption>Figure 3.3: Examples of activation functions that can be used in neural networks.</figcaption></figure>

However, a single perceptron is a relatively simple function with limited capacity to repre-
sent complex functions. To overcome this limitation, multiple perceptrons are grouped together
to form a *layer*. By subsequently stacking several layers of perceptrons, one obtains what is
known as a neural network of type MLP. Three types of layers are generally distinguished:

1. **The input layer**, which receives the input data and transmits it to the network;

2. **The hidden layers**, located between the input layer and the output layer, which extract increasingly abstract and complex features from the data. More than one hidden layer may be used;

3. **The output layer**, which produces the final results of the network, depending on the task (classification, regression, etc.).

Figure 3.4 illustrates an example of an MLP neural network, consisting of an input layer with five neurons, an output layer with three neurons, and two hidden layers each containing four neurons. An MLP neural network has a large number of parameters (weights and biases of the perceptrons) that must be adjusted. Hornik (1989) [119] showed that this type of network is a universal function approximator, meaning that it can theoretically approximate any function, provided that a sufficiently large architecture and an appropriate parameter configuration are available. The output produced by this network, referred to as the *prediction*, depends directly
<figure><img src="image_23.png" /><figcaption>Figure 3.4: Example of an MLP with five inputs and three outputs. The network consists of an input layer, two hidden layers, and an output layer. The red circles represent the neurons (or perceptrons), and $\sigma$ denotes the activation function applied to each of them. The same activation function may be used for all neurons, or different functions may be employed depending on the specific needs of the model.</figcaption></figure>

on these parameters. The process of optimizing the parameters to achieve the best possible performance on a given task constitutes the essence of neural network *training*.

### 3.1.2 Training

Training consists in optimizing the weights and biases of the perceptrons in the neural network in order to minimize the error made by the network. First, one must select the data used to train the neural network. In general, these data are divided into three datasets:

1. **The training set**, which contains labeled examples used to train the neural network;

2. **The validation set**, which is used to evaluate the generalization capability of the model and to tune the hyperparameters in order to avoid overfitting<fnref n="2" />;

3. **The test set**, which is completely independent of the previous datasets and is used to assess the final performance of the model on unseen data.

To evaluate the overall quality of the predictions produced by the neural network, a *cost function* is introduced. This function measures the global error made by the model. Let us assume that the dataset contains *N* samples, that $y^{(i)}$ denotes the expected output of the neural

<footnote n="2">²See Section 3.1.3. *Underfitting and overfitting* for the definition of overfitting.</footnote>
network for the input $x^{(i)}$, and that $h_{\theta}(x^{(i)})$ represents the prediction produced by the network for the same input.

Several cost functions are commonly used in artificial intelligence; some of them are presented below.

### Mean Squared Error (MSE)

MSE is the most commonly used option. It strongly penalizes large errors, which makes it sensitive to outliers.

$$ \mathcal{L}(Y, h_{\theta}(X)) = \frac{1}{N} \sum_{i=1}^{N} [h_{\theta}(x^{(i)}) - y^{(i)}]^{2} . $$

### Mean Absolute Error (MAE)

Unlike MSE, MAE penalizes errors linearly, which makes it more robust to outliers. It is frequently used in regression tasks.

$$ \mathcal{L}(Y, h_{\theta}(X)) = \frac{1}{N} \sum_{i=1}^{N} \left| h_{\theta}(x^{(i)}) - y^{(i)} \right| . $$

### Logarithmic Mean Squared Error (LMSE)

This cost function is particularly well suited when the data are strictly positive and when minimizing relative error is more important than minimizing absolute error. It is often used in problems where target values span several orders of magnitude.

$$ \mathcal{L}(Y, h_{\theta}(X)) = \frac{1}{N} \sum_{i=1}^{N} \left[ \ln(h_{\theta}(x^{(i)})) - \ln(y^{(i)}) \right]^{2} . $$

### Binary Cross-Entropy

Used in binary classification tasks, this function measures the divergence between the predicted probabilities and the true labels $y^{(i)}$, which are assumed to belong to {0, 1}. It requires a network whose outputs lie in the interval [0; 1].

$$ \mathcal{L}(Y, h_{\theta}(X)) = -\frac{1}{N} \sum_{i=1}^{N} \left[ y^{(i)} \ln(h_{\theta}(x^{(i)})) + (1 - y^{(i)}) \ln(1 - h_{\theta}(x^{(i)})) \right] . $$

Next, in order to train the neural network, it is essential to choose an appropriate optimization algorithm capable of determining the weights and biases that minimize the value of the cost function defined above. These algorithms adjust the network parameters to achieve good performance on the training data while preserving generalization capability on the validation and test sets. Some commonly used optimizers are presented below.
### Gradient Descent

This classical method updates the parameters using the gradient of the cost function computed over the entire training dataset. The parameter update is written as:

$$ \theta_{t+1} = \theta_t - \eta \nabla_{\theta} \mathcal{L}(Y, h_{\theta_t}(X)) , $$

where $\theta_t$ denotes the neural network parameters at iteration $t$, $\eta$ is the learning rate, which must be tuned, $\nabla_{\theta}$ is the gradient with respect to the network parameters $\theta$, and $\mathcal{L}$ is the cost function evaluated over the entire training dataset. Each update corresponds to one epoch. At each epoch, the value of the cost function $\mathcal{L}$ decreases. However, this method is sensitive to noisy data, which can degrade the quality of the final neural network.

### Stochastic Gradient Descent (SGD)

Rather than computing the gradient over the full dataset, SGD updates the neural network parameters using each individual training example. For a training set containing $N$ samples, the update corresponding to the $i$-th example is:

$$ \theta_{t+1} = \theta_t - \eta \nabla_{\theta} \mathcal{L}(y^{(i)}, h_{\theta_t}(x^{(i)})) , $$

One epoch therefore corresponds to $N$ updates. This method is, however, more computationally expensive and difficult to parallelize.

### Mini-batch Gradient Descent

This approach lies between classical gradient descent and SGD. It divides the training dataset into small batches (mini-batches) of size $K$ and performs one update per mini-batch. If $\mathcal{B}_K$ denotes the set of all mini-batches, the updates are written as:

$$ \forall (X_K, Y_K) \in \mathcal{B}_K, \theta_{t+1} = \theta_t - \eta \nabla_{\theta} \mathcal{L}(Y_K, h_{\theta_t}(X_K)) , $$

Thus, $N/K$ updates are performed per epoch. This method reduces computation time compared to SGD, is well suited to parallelization, and retains a stochastic character, making training less sensitive to noisy data than classical gradient descent.

### Adaptive Moment Estimation (Adam)

Proposed by Kingma in 2014 [120], Adam adapts the learning rate for each parameter by combining mini-batch gradient descent with adjustments based on the moments of the gradient. Two variables are defined:

$$ v_t = \beta_1 v_{t-1} + (1 - \beta_1) g_t , $$

$$ s_t = \beta_2 s_{t-1} + (1 - \beta_2) g_t^2 , $$
where $g_t$ is the gradient of the cost function evaluated on the mini-batch data, and $\beta_1$ and $\beta_2$ are hyperparameters (usually set to $\beta_1 = 0.9$ and $\beta_2 = 0.999$). In the Adam algorithm, normalized values of $v_t$ and $s_t$ are used to avoid bias:

$$ \hat{v}_t = \frac{v_t}{1 - \beta_1^t}, \qquad \hat{s}_t = \frac{s_t}{1 - \beta_2^t}. $$

The update of the neural network parameters is then given by:

$$ \theta_{t+1} = \theta_t - \eta_t \hat{v}_t, $$

with

$$ \eta_t = \frac{\eta}{\sqrt{\hat{s}_t + \epsilon}}, $$

where $\epsilon$ is a small constant (e.g. $10^{-6}$) introduced to avoid division by zero. Adam is widely used due to its fast convergence.

## Limited-memory Broyden-Fletcher-Goldfarb-Shano (L-BFGS)

This quasi-Newton method uses an approximation of the Hessian matrix to accelerate convergence while limiting memory usage. In Newton's method, the update would be computed as:

$$ \theta_{t+1} = \theta_t - [\vec{\nabla}_{\theta}^2 \mathcal{L}(Y, h_{\theta_t}(X))]^{-1} \vec{\nabla}_{\theta} \mathcal{L}(Y, h_{\theta_t}(X)), $$

where $[\vec{\nabla}_{\theta}^2 \mathcal{L}(Y, h_{\theta_t}(X))]^{-1}$ and $\vec{\nabla}_{\theta} \mathcal{L}(Y, h_{\theta_t}(X))$ denote respectively the inverse Hessian and the gradient of the cost function with respect to the neural network parameters, evaluated on the training dataset. L-BFGS replaces the Hessian with a matrix $B_t$ (generally initialized as the identity matrix), which approximates the Hessian. The algorithm then proceeds through the following steps:

1. Computation of the descent direction:

$$ p_t = -B_t^{-1} \vec{\nabla}_{\theta} \mathcal{L}(Y, h_{\theta_t}(X)), $$

2. Line search to determine the step size $\alpha_t$ such that $\alpha_t = \argmin(\mathcal{L}(\theta_t + \alpha_t p_t))$,

3. Update of the neural network parameters:

$$ \theta_{t+1} = \theta_t + \alpha_t p_t, $$

4. Computation of the vectors $s_t$ and $y_t$:

$$ s_t = \theta_{t+1} - \theta_t, \qquad y_t = \vec{\nabla}_{\theta} \mathcal{L}(Y, h_{\theta_{t+1}}(X)) - \vec{\nabla}_{\theta} \mathcal{L}(Y, h_{\theta_t}(X)), $$
5. Update of the inverse Hessian approximation:

$$ \mathbb{B}_{t+1}^{-1} = \mathbb{B}_{t}^{-1} + \frac{(s_{t}^{T} y_{t} + y_{t}^{T} \mathbb{B}_{t}^{-1} y_{t}) s_{t} s_{t}^{T}}{(s_{t}^{T} y_{t})^{2}} - \frac{\mathbb{B}_{t}^{-1} y_{t} s_{t}^{T} + s_{t} y_{t}^{T} \mathbb{B}_{t}^{-1}}{s_{t}^{T} y_{t}} $$

L-BFGS converges much faster than the previous algorithms, but may suffer from the same sensitivity to noisy data as gradient descent. It is often used to refine the parameters after an initial optimization phase using the Adam optimizer.

### 3.1.3 Underfitting and overfitting

At the end of the training of a neural network, three different situations may arise:

* **Underfitting:** The neural network fails to correctly solve the problem, including on the training set itself. This phenomenon may result from an overly simple architecture, which is unable to capture the complexity of the problem, or from an early termination of the training process;

* **Overfitting:** The neural network fits the training data perfectly but fails to provide reliable predictions on unseen cases. This may occur if the network architecture is too complex or if the training is continued for too long. In this case, the network is said to fail to *generalize*;

* **Satisfactory learning:** The neural network produces predictions with a similar level of error on the training set and on new data. The network is then said to *generalize* correctly.

Figure 3.5 illustrates three typical configurations encountered in the classification of points into two sets: underfitting, satisfactory learning, and overfitting. In the case of underfitting, the function used is too simple, which prevents the model from properly separating the two groups of points. Conversely, in the case of overfitting, the model adapts excessively to the training data, perfectly classifying all points, but at the expense of poor generalization: it also classifies as belonging to the sets certain points that are clearly outliers. Finally, the best compromise is achieved by a model capable of correctly classifying the majority of points while discerning those that do not belong to the main sets.

To limit underfitting, the network should be trained for a sufficient number of epochs, until the solution converges toward a stable error during training. It is also important to use
<figure><img src="image_24.png" /><figcaption>Figure 3.5: Different possible outcomes after the training of a neural network.</figcaption></figure>

an architecture that is sufficiently expressive to capture the dynamics of the problem, while avoiding excessive complexity.

In the case of overfitting, several strategies can be implemented. The first consists in using a stochastic optimization algorithm, such as Adam, which is known to improve generalization capabilities. Two additional techniques are also commonly employed:

1. A common approach consists in adding a regularization term $\mathcal{L}(\theta)$ to the cost function, which explicitly depends on the parameters of the neural network. This term aims to encourage the weights and biases to remain close to zero, thereby improving the generalization capability. The modified cost function is then expressed:

$$ \mathcal{L}^R(Y, h_\theta(X)) = \mathcal{L}(Y, h_\theta(X)) + \mathcal{L}^*(\theta) . $$

Two common forms of regularization are used:

a) **L1 regularization:** penalizes the sum of the absolute values of the parameters:

$$ \mathcal{L}^*(\theta) = \lambda \sum_i |\theta_i| . $$

It promotes a more parsimonious model by inducing sparsity through the cancellation of some parameters, which may improve the interpretability of the neural network;

b) **L2 regularization:** penalizes the sum of the squares of the weights:

$$ \mathcal{L}^*(\theta) = \lambda \sum_i \theta_i^2 . $$

It reduces the magnitude of the weights without setting them exactly to zero, making the network more stable with respect to variations in the input data.
In both cases, $\lambda$ is a hyperparameter that must be tuned depending on the problem being addressed, and the $\theta_i$ denote the parameters of the neural network;

2. **Dropout:** This technique consists in randomly deactivating a certain percentage of perceptrons at each training iteration. In practice, during the parameter update step, a subset of neurons is temporarily ignored, which prevents the network from relying too strongly on relationships that are specific to the training data. This approach helps limit overfitting and improves the robustness of the model. For more details on dropout, see Baldi (2014) [121].

## 3.2 Eddington factor in the M1-multigroup model

Let us briefly recall the definition of the Eddington factor.
Definition of the Eddington factor $\chi_g$

Within the framework of the M1-multigroup model for radiative transfer, the Eddington factor $\chi_g$ associated with the spectral group g, defined between the frequencies $\nu_1$ and $\nu_2$, is a scalar function that relates the radiative pressure tensor $P_g$ to the radiative energy $E_g$ of the group. The closure relation reads:

$$P_g = \left[ \frac{1 - \chi_g}{2} \mathbb{I}_d + \frac{3\chi_g - 1}{2} \frac{\vec{F}_g \otimes \vec{F}_g}{||\vec{F}_g||^2} \right] E_g,$$

where $\vec{F}_g$ is the radiative flux vector of the group and $\mathbb{I}_d$ denotes the identity tensor.

This expression makes it possible to express the radiative pressure knowing lower-order moments and thereby to close the system of the moment equations. The Eddington factor $\chi_g$ depends on three parameters:

* the reduced flux $f_g = ||\vec{F}_g||/(c E_g)$,

* the dimensionless radiative temperature $\mathcal{T}_g = k_B E_g^{1/4}/(a_R^{1/4} h \tilde{\nu})$,

* the group narrowness $\delta_g = \nu_1/\nu_2$, when $\nu_1 > 0$ and $\nu_2 < +\infty$.

The reference frequency $\tilde{\nu}$ may be chosen as $\nu_1$, $\nu_2$, or $\sqrt{\nu_1\nu_2}$, depending on the values of $\nu_1$ and $\nu_2$.

The limiting values of $\chi_g$ correspond to well-identified physical regimes:

$$\chi_g \xrightarrow{f_g \to 0} \frac{1}{3} \quad \text{(isotropic regime) et} \quad \chi_g \xrightarrow{f_g \to 1} 1 \quad \text{(free-streaming regime).}$$

As indicated in Section 2.2.2. *The moment equations*, there is no general analytical formula for computing the Eddington factor within the framework of the M1-multigroup model. However, it can be computed numerically (see Appendix C. *Computing the Eddington factor using search algorithms*). In this section, I present several methods used in the literature to perform this computation, followed by the innovative method based on AI that I developed during my PhD. Finally, I compare the performance of this new approach with that of existing methods.
### 3.2.1 “Historical” methods

Historically, several approaches have been developed to estimate the Eddington factor in the M1-multigroup model. These methods aim to circumvent the absence of a general analytical formula by exploiting different computational and approximation strategies. Three main families of techniques have thus emerged:

1. search algorithms;

2. interpolation of precomputed values;

3. the use of the Eddington factor expression from the M1-gray model.

Each of these methods present advantages and limitations, which affect the accuracy and efficiency of radiative simulations, as will be shown in the following sections.

#### ❖ Search algorithms

As discussed in Section 2.2.3. *The Eddington factor*, in the general case the radiative energy, the reduced flux, and the Eddington factor all depend on the Lagrange multipliers $\alpha_{0,g}$ and $\beta_g$, as well as on the frequency bounds of the group $v_{g-1/2}$ and $v_{g+1/2}$:

$$E_g = e(\alpha_{0,g}, \beta_g; v_{g-1/2}, v_{g+1/2}),$$

$$f_g = f(\alpha_{0,g}, \beta_g; v_{g-1/2}, v_{g+1/2}),$$

$$\chi_g = c(\alpha_{0,g}, \beta_g; v_{g-1/2}, v_{g+1/2}).$$

Turpault [93] showed that, for each physically admissible pair of radiative energy and reduced flux ($E_g, f_g$), there exists a unique pair of Lagrange multipliers ($\alpha_{0,g}, \beta_g$). A classical approach therefore consists in using a search algorithm to determine these multipliers and then deducing the associated Eddington factor $\chi_g$. Two search algorithms have been developed to solve this problem:

1. **The Dichotomy-Newton algorithm**, introduced in the PhD thesis of Hung Chinh Nguyen (2011) [37], which combines a bisection search to estimate $\beta_g$ within the interval $]-1; 1[$, followed by the Newton method to determine $\alpha_{0,g}$ in $\mathbb{R}^+$;

2. **The line-search algorithm**, developed in the course of this thesis, which simultaneously searches $\beta_g$ and $\alpha_{0,g}$ in the space $]-1; 1[$ $\times \mathbb{R}^+$ in order to accelerate convergence and reduce computational cost.
Additional details on the implementation of these search algorithms are provided in Appendix C. *Computing the Eddington factor using search algorithms*. Although this approach yields highly accurate estimates of the Eddington factor, with a relative error on the order of $10^{-5}$ %, it remains computationally expensive, even when using the line-search algorithm, which converges faster than the Dichotomy-Newton algorithm. Moreover, in some cases the algorithm may fail to converge, leading to simulation failure.

### ❖ Interpolation

Due to the prohibitive computational cost of search algorithms, an alternative approach based on the precomputation and interpolation of Eddington factors was proposed by Turpault in 2003, and later adopted and refined within the HADES code by Nguyen in 2011 [97, 37]. This approach aims to reconcile accuracy and computational efficiency by replacing systematic numerical resolution with interpolation over a grid of precomputed values.

In the HADES code, interpolation is performed using cubic splines, based on a grid of precomputed Eddington factor values for a given group. These values are obtained using the search algorithms presented previously and are stored for regularly spaced values of the logarithm of the radiative energy and the reduced flux (that is, $\ln(E_g)$ and $f_g$). This method enables a rapid estimation of the Eddington factor while maintaining good accuracy, provided that the grid is sufficiently dense and properly calibrated.

However, this approach has a major limitation: it relies on a good a priori knowledge of the orders of magnitude of the radiative quantities encountered during a simulation. In the context of the M1-multigroup model, these quantities are difficult to estimate, since they depend on spectral interactions with matter, which vary from one group to another and evolve over time. Insufficient coverage of the interpolation grid can then lead to significant errors in the estimation of the Eddington factor, thereby compromising the accuracy of the simulations.

### ❖ Analytical expression of the M1-gray model

In the literature, some authors, such as Vaytet in his work on the implementation of the M1-multigroup model in the HERACLES code [35, 122], adopt a simplified approach by neglecting the dependence of the Eddington factor on the dimensionless radiative temperature $T_g$ and on the group narrowness $\delta_g$. Instead, they use the analytical expression derived from
the M1-gray model (see equation (2.71)):

$$ \chi_g \approx \frac{3 + 4f_g^2}{5 + 2\sqrt{4 - 3f_g^2}}, $$

This method has the advantage of being extremely fast, since it is fully analytical. However, the expression used is not generalizable within the M1-multigroup framework and can introduce significant errors, reaching up to 16 % relative error in the computation of the Eddington factor. While this approach may remain reasonably accurate for broad frequency groups, its precision rapidly degrades as the spectral resolution increases, leading to increasingly inaccurate predictions for narrow groups.

### 3.2.2 Description of the new method based on Artificial Intelligence

The objective of this study is to develop a new method based on AI to accurately estimate the Eddington factor while maintaining a reasonable computational cost. Before implementing this approach, it is essential to select the AI methodology best suited to the problem. Three types of models were considered:

1. **Symbolic regression [123]** aims at identifying mathematical expressions that fit the data. This approach has the advantage of producing interpretable models that are optimized for numerical evaluation. However, my experiments showed that it is difficult to simultaneously capture all the dependencies of the Eddington factor on the different parameters of the problem;

2. **Gaussian processes [124]** are probabilistic models that provide predictions together with an uncertainty quantification. They are particularly effective trained on small datasets, but their computational cost becomes prohibitive using larger datasets. In the present case, the prediction time associated with Gaussian processes proved to be too high;

3. **Neural networks (introduced previously)** are well suited to complex nonlinear relationships. They offer excellent accuracy and good generalization capabilities, at the cost of reduced interpretability.

In light of these considerations, the use of neural networks, and more specifically a MLP-type architecture, was selected as the preferred approach (see Radureau (2025) [98] for further details). Two main strategies can be considered in this context. The first consists in training
<figure><img src="image_25.png" /><figcaption>Figure 3.6: Schematic overview of the developed method, for the case where $v_{g-1/2} > 0$ and $v_{g+1/2} < +\infty$. $\chi_{g,L}$ and $\chi_{g,H}$ correspond respectively to the polynomials defined by equations (3.6) and (3.5), used to compute $\chi_g$ in domains L and H, while $\chi_R$ denotes the Eddington factor of the M1-gray model (see equation (2.71)), used in domain $I_b$. The MLPs are used to approximate the closure in the intermediate domains $I_a$, $I_c$, $I_d$, $I_e$, and $I_f$. The boundaries of domains L, H, and $I_a$ to $I_f$ are defined in table 3.1. The brown lines delineate the domains of applicability of the different approaches.</figcaption></figure>

a single large network capable of predicting the Eddington factor over the entire domain, in the spirit of general-purpose approaches such as ChatGPT. However, this solution does not meet the objectives set here, as it involves large models whose predictions are computationally expensive. The second strategy relies on the use of more compact networks, specialized on specific subdomains, whose task is essentially limited to interpolation. With the goal of achieving numerical efficiency while maintaining a high level of accuracy, this second approach was therefore favored.

Figure 3.6 illustrates the organization of the method, distinguishing the different strategies adopted depending on the domain. In domains L and H, the Eddington factor is evaluated using polynomials. In domain $I_b$, the analytical expression from the M1-gray model is used. The
<figure><img src="image_26.png" /></figure>

(a) $\delta_g = 10^{-4}$. L, H, I<sub>b</sub>: Domains in which $\chi_g = \mathcal{F}(f_g)$; I<sub>a</sub>: Domain in which $\chi_g = \mathcal{F}(f_g, \mathcal{T}_g^{\nu_{g-1/2}})$; I<sub>c</sub>: Domain in which $\chi_g = \mathcal{F}(f_g, \mathcal{T}_g^{\nu_{g+1/2}})$.

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

(c) $\delta_g = 0.3$. L, H: Domains in which $\chi_g = \mathcal{F}(f_g)$; I<sub>d</sub>: Domain in which $\chi_g = \mathcal{F}(f_g, \mathcal{T}_g^{\nu_{g-1/2}})$; I<sub>e</sub>: Domain in which $\chi_g = \mathcal{F}(f_g, \mathcal{T}_g^{\nu_{g-1/2}}, \delta_g)$.

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

(b) $\delta_g = 0.07$. L, H: Domains in which $\chi_g = \mathcal{F}(f_g)$; I<sub>d</sub>: Domain in which $\chi_g = \mathcal{F}(f_g, \mathcal{T}_g^{\nu_{g-1/2}})$; I<sub>e</sub>: Domain in which $\chi_g = \mathcal{F}(f_g, \mathcal{T}_g^{\nu_{g-1/2}}, \delta_g)$.

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

(d) $\delta_g = 0.99$. L, H: Domains in which $\chi_g = \mathcal{F}(f_g)$; I<sub>f</sub>: Domain in which $\chi_g = \mathcal{F}(f_g, \mathcal{T}_g^s)$.

Figure 3.7: Ratio between the Eddington factor of the M1-multigroup model, $\chi_g$, and that of the M1-gray model, $\chi_R$, for a reduced flux $f_g = 0.65$, in the case where $\nu_{g-1/2} > 0$ and $\nu_{g+1/2} < +\infty$, for different values of the group narrowness $\delta_g$.

remaining domains are handled using five distinct MLPs, each with two or three inputs: one network specialized for domains I<sub>a</sub> and I<sub>d</sub>; one for domain I<sub>c</sub>; one for I<sub>f</sub>; and two for domain I<sub>e</sub>: one for cases where $\delta_g \le 0.03$, and the other for cases where $\delta_g > 0.03$. The remainder of this section is devoted to a detailed description of this approach.

To design this method, I first analyzed the dependence of the Eddington factor on the influential parameters, as detailed in Section 2.2.3. *The Eddington factor*. The main results of this analysis are summarized below:

* When $\nu_{g-1/2} > 0$ and $\nu_{g+1/2} < +\infty$ (figure 3.7), the Eddington factor depends on the dimensionless radiative temperature $\mathcal{T}_g$ (or, equivalently, $\mathcal{T}_g^{\nu_{g+1/2}}$ or $\mathcal{T}_g^{\nu_{g-1/2}}$), on the
<figure><img src="image_30.png" /></figure>

(a) Case where $v_{g+1/2} = +\infty$. L, I<sub>b</sub>: Domains in which $\chi_g = \mathcal{F}(f_g)$: I<sub>a</sub>: Domain in which $\chi_g = \mathcal{F}(f_g, \mathcal{T}_g^{v_{g-1/2}})$.

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

(b) Case where $v_{g-1/2} = 0$. H, I<sub>b</sub>: Domains in which $\chi_g = \mathcal{F}(f_g)$: I<sub>c</sub>: Domain in which $\chi_g = \mathcal{F}(f_g, \mathcal{T}_g^{v_{g+1/2}})$.

Figure 3.8: Ratio between the Eddington factor of the M1-multigroup model, $\chi_g$, and that of the M1-gray model, $\chi_R$, for a reduced flux $f_g = 0.65$, in the cases where $v_{g-1/2} = 0$ and where $v_{g+1/2} = +\infty$.

reduced flux $f_g$, and on the group narrowness $\delta_g$;

* When $v_{g-1/2} = 0$ (figure 3.8b), the Eddington factor depends only on the dimensionless radiative temperature $\mathcal{T}_g^{v_{g+1/2}}$ and on the reduced flux $f_g$;

* When $v_{g+1/2} = +\infty$ (figure 3.8a), it depends on the dimensionless radiative temperature $\mathcal{T}_g^{v_{g-1/2}}$ and on the reduced flux $f_g$.

On this basis, I identified several domains in which the Eddington factor exhibits distinct dependencies:

* **domains L and H (figures 3.7 and 3.8):** asymptotic cases in which the Eddington factor depends only on the reduced flux $f_g$,

* **domain I<sub>b</sub> (figures 3.7a and 3.8):** the Eddington factor depends only on the reduced flux $f_g$, and the analytical expression of the Eddington factor from the M1-gray model is valid,

* **domain I<sub>a</sub> (figures 3.7a and 3.8a):** the Eddington factor depends only on the reduced flux $f_g$ and on the dimensionless radiative temperature $\mathcal{T}_g^{v_{g-1/2}}$,

* **domain I<sub>c</sub> (figures 3.7a and 3.8b):** the Eddington factor depends only on the reduced flux $f_g$ and on the dimensionless radiative temperature $\mathcal{T}_g^{v_{g+1/2}}$,

* **domain I<sub>d</sub> (figures 3.7b and 3.7c):** the Eddington factor depends only on the reduced flux $f_g$, on the dimensionless radiative temperature $\mathcal{T}_g^{v_{g-1/2}}$, and very weakly on the group
narrowness $\delta_g$. I neglected this last dependence and assumed that the Eddington factor values in this domain are the same as in domain $I_a$, which is accurate within a relative error of $10^{-3}$ %,

* **domain $I_e$ (figures 3.7b and 3.7c):** the Eddington factor depends strongly on the reduced flux $f_g$, on the dimensionless radiative temperature $T_g^{\nu_g+1/2}$, and on the group narrowness $\delta_g$,

* **domain $I_f$ (figure 3.7d):** the Eddington factor depends on the reduced flux $f_g$ and on a dimensionless radiative temperature $T_g^s$, which will be referred to as the *shifted temperature*, and which accounts for the shift described in equation (2.79), Section 2.2.3. *The Eddington factor*.

First, let us define the temperature $T_g^*$, related to the dimensionless temperature $T_g$ by:

$$T_g^* = \frac{T_g \sqrt{\delta_g}}{(1 - \delta_g^{3/2})^{1/4}} = \frac{T_g^{\nu_g - 1/2}}{(1 - \delta_g^{3/2})^{1/4}}$$

The shifted temperature can then be defined from $T_g^*$ as follows<fnref n="3" />:

$$T_g^s = T_g^* \frac{(1 - 0.99^3)^{1/4}}{\sqrt{0.99}} = \frac{T_g \sqrt{\delta_g / 0.99}}{\{(1 - \delta_g^3) / (1 - 0.99^3)\}^{1/4}} \quad (3.1)$$

This expression ensures that $T_g^s = T_g$ when $\delta_g = 0.99$.

The boundaries between the different domains H, L, etc. were determined empirically, so that the relative errors induced by the dependency simplifications remain below $10^{-3}$ %. For domains L and H, polynomial expressions are used; for domain $I_b$, the analytical expression from the M1-gray model; and for all other domains, neural networks are used. These choices are summarized in table 3.1.

Before detailing the method, it is crucial to examine the impact of an error in the estimation of the Eddington factor on the components of the radiative pressure tensor. In the context of 2D simulations performed with the HADES code, since the radiative pressure tensor is symmetric, it is sufficient to compute three of its components to fully determine the tensor: the components $P_{g,xx}$, $P_{g,yy}$, and $P_{g,xy}$. Using equation (2.57), each component of the radiative pressure tensor can be expressed as:

$$P_{g,xx} = \left\{ \frac{1 - \chi_g}{2} + \frac{3\chi_g - 1}{2} \frac{F_{g,x}^2}{||F_g||^2} \right\} E_g,$$

<footnote n="3">³This definition differs slightly from the one presented in the article [98], which relied on linear and quadratic adjustments. This difference does not affect the results.</footnote>
<table><thead><tr><th rowspan="2"></th><th colspan="2">Range of δ<sub>g</sub></th><th colspan="2">Range of T<sub>g</sub></th><th rowspan="2">Method **</th></tr><tr><th>Min</th><th>Max</th><th>Min</th><th>Max</th></tr></thead><tbody><tr><td rowspan="2">L</td><td>0</td><td>0.9</td><td>0</td><td>10<sup>-5</sup>√δ<sub>g</sub></td><td>P</td></tr><tr><td>0.9</td><td>1</td><td>0*</td><td>10<sup>-3</sup>/√δ<sub>g</sub>*</td><td>Input: f<sub>g</sub></td></tr><tr><td rowspan="2">H</td><td>0</td><td>0.9</td><td>10<sup>0.6</sup>/√δ<sub>g</sub></td><td>+∞</td><td>P</td></tr><tr><td>0.9</td><td>1</td><td>10<sup>0.6</sup>/√δ<sub>g</sub>*</td><td>+∞*</td><td>Input: f<sub>g</sub></td></tr><tr><td rowspan="2">I<sub>a</sub></td><td>0</td><td>10<sup>-2.2</sup></td><td>10<sup>-5</sup>√δ<sub>g</sub></td><td>10<sup>1.1</sup>√δ<sub>g</sub></td><td>MLP</td></tr><tr><td>10<sup>-2.2</sup></td><td>0.03</td><td>10<sup>-5</sup>√δ<sub>g</sub></td><td>1</td><td>Input: f<sub>g</sub>, T<sub>g</sub><sup>Yg-1/2</sup></td></tr><tr><td rowspan="2">I<sub>b</sub></td><td>0</td><td>10<sup>-3</sup></td><td>10<sup>1.1</sup>√δ<sub>g</sub></td><td>10<sup>-1.5</sup>/√δ<sub>g</sub></td><td>χ<sub>R</sub></td></tr><tr><td>10<sup>-3</sup></td><td>10<sup>-2.2</sup></td><td>10<sup>1.1</sup>√δ<sub>g</sub></td><td>1</td><td></td></tr><tr><td rowspan="2">I<sub>c</sub></td><td>0</td><td>10<sup>-3</sup></td><td>10<sup>-1.5</sup>/√δ<sub>g</sub></td><td>10<sup>0.6</sup>/√δ<sub>g</sub></td><td>MLP</td></tr><tr><td>10<sup>-3</sup></td><td>0.03</td><td>1</td><td>10<sup>0.6</sup>/√δ<sub>g</sub></td><td>Input: f<sub>g</sub>, T<sub>g</sub><sup>Yg+1/2</sup></td></tr><tr><td rowspan="2">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</td></tr><tr><td></td><td></td><td></td><td></td><td>Input: f<sub>g</sub>, T<sub>g</sub><sup>Yg+1/2</sup></td></tr><tr><td rowspan="2">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</td></tr><tr><td></td><td></td><td></td><td></td><td>Input: f<sub>g</sub>, δ<sub>g</sub>, T<sub>g</sub><sup>Yg+1/2</sup></td></tr><tr><td rowspan="2">I<sub>f</sub></td><td>0.9</td><td>1</td><td>10<sup>-3</sup>/√δ<sub>g</sub>*</td><td>10<sup>0.6</sup>/√δ<sub>g</sub>*</td><td>MLP</td></tr><tr><td></td><td></td><td></td><td></td><td>Input: f<sub>g</sub>, T<sub>g</sub><sup>s</sup></td></tr></tbody></table>

<caption>Table 3.1: Limits for each domain and the method used to compute the Eddington factor.</caption>

$$ \mathbb{P}_{g,yy} = \left\{ \frac{1 - \chi_g}{2} + \frac{3\chi_g - 1}{2} \frac{F_{g,y}^2}{||F_g||^2} \right\} E_g, $$

$$ \mathbb{P}_{g,xy} = \frac{3\chi_g - 1}{2} \frac{F_{g,x}F_{g,y}}{||F_g||^2} E_g, $$

When considering an error δχ<sub>g</sub> in the estimation of the Eddington factor, the relative error on the components of the radiative pressure can be expressed as:

$$ \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, $$

$$ \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, $$

<footnote>*Min/max values of the shifted temperature T<sub>g</sub><sup>s</sup> (see equation 3.1).</footnote>

<footnote>**P: polynomial, MLP: neural network, χ<sub>R</sub>: Eddington factor of the M1-gray model. The variables listed after “Input:” correspond to the inputs used for each method.</footnote>
$$ \left| \frac{\delta \mathbb{P}_{g,xy}}{\mathbb{P}_{g,xy}} \right| = \frac{\delta \chi_g}{\chi_g - 1/3}, $$

where $x = F_{g,x}/||F_g||$ and $y = F_{g,y}/||F_g||$. The quantities $|\delta \mathbb{P}_{g,xx}/\mathbb{P}_{g,xx}|$ and $|\delta \mathbb{P}_{g,yy}/\mathbb{P}_{g,yy}|$ reach their maximum relative error when $x = 0$ or $1$ and $y = 0$ or $1$, respectively. In this case, the maximum relative errors can be written as:

$$ \left\{ \begin{array}{l} \left| \frac{\delta \mathbb{P}_{g,xx}}{\mathbb{P}_{g,xx}} \right|_{\max} = \frac{\delta \chi_g}{\min \{ \chi_g, 1 - \chi_g \}}, \\ \\ \left| \frac{\delta \mathbb{P}_{g,yy}}{\mathbb{P}_{g,yy}} \right|_{\max} = \frac{\delta \chi_g}{\min \{ \chi_g, 1 - \chi_g \}}, \\ \\ \left| \frac{\delta \mathbb{P}_{g,xy}}{\mathbb{P}_{g,xy}} \right|_{\max} = \frac{\delta \chi_g}{\chi_g - 1/3}. \end{array} \right. \qquad (3.2) $$

These expressions highlight the importance of an accurate estimation of the Eddington factor, particularly when the reduced flux $f_g$ is close to 0 or 1, in order to minimize relative errors in the components of the radiative pressure tensor. When $f_g$ is close to 0, the component $\mathbb{P}_{g,xy}$ also tends toward zero. Moreover, if $f_g$ is close to 0 and the radiative flux in the $x$ or $y$ direction is weak, the components $\mathbb{P}_{g,xx}$ or $\mathbb{P}_{g,yy}$ also approach zero, which explains the large relative errors observed in such situations. The maximum error on the components of the radiative pressure tensor will therefore be evaluated using the expressions derived above. In the following, I will detail the methods I have developed to compute the Eddington factor in the different domains, taking this impact into account.

❖ Asymptotically high temperatures (domain H, $k_B T_g \gg h \nu_{g+1/2}$)

In this case, equations (2.77) and (2.78), presented in Section 2.2.3. *The Eddington factor*, simplify in such a way that the reduced flux and the Eddington factor can be expressed as follows:

$$ f_g = - \frac{\operatorname{arctanh}(\beta_g) - \beta_g}{\beta_g \operatorname{arctanh}(\beta_g)}, \qquad (3.3) $$

$$ \chi_g = \frac{\operatorname{arctanh}(\beta_g) - \beta_g}{\beta_g^2 \operatorname{arctanh}(\beta_g)}. \qquad (3.4) $$

Here, it is possible to use a simple bisection algorithm to determine the Lagrange multiplier $\beta_g$ within the interval $]-1; 1[$, and then to compute the associated Eddington factor. However, performing this search during a simulation is not optimal in terms of computational cost.
<figure><img src="image_32.png" /><figcaption>Figure 3.9: Comparison of the Eddington factor computed using the polynomial expression $\chi_{g,H}$, defined by equation (3.5), with the reference data obtained by solving equations (3.3) and (3.4) using a bisection algorithm.</figcaption></figure>

Therefore, I chose to use a 12th-order polynomial, whose expression is:

$$ \chi_{g,H}(f_g) = 1/3 + f_g^2 - f_g^3/3 + 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). \quad (3.5) $$

In this formulation, the Eddington factor was constrained to take the value 1/3 at $f_g = 0$ and 1 at $f_g = 1$. In addition, I imposed the conditions $\chi'_{g,H}(0) = 0$ and $\chi'_{g,H}(1) = 1$ in order to ensure appropriate regularity at the boundaries of the domain. This expression guarantees excellent accuracy, with a relative error below $10^{-3}$ % in the estimation of the Eddington factor. The maximum observed error is of the order of $10^{-2}$ % for the $P_{g,xx}$ and $P_{g,yy}$ components of the radiative pressure tensor, and of $10^{-1}$ % for the $P_{g,xy}$ component. Figure 3.9 shows a comparison between the results obtained by numerically solving equations (3.3) (green curve) and (3.4) using a bisection algorithm, and those obtained from the polynomial approximation (3.5) (red curve). Very good agreement between the two curves is observed.
❖ Asymptotically low temperatures (domain L, $k_B T_g \ll h\nu_{g-1/2}$)

When the radiative temperature is very low, there is unfortunately no simplified analytical expression that allows the Eddington factor $\chi_g$ to be computed directly as a function of the reduced flux $f_g$ alone, unlike in the previous case. This means that, in this regime, the radiative temperature still influences $\chi_g$. However, an analysis of the curves shown in figures 3.7 and 3.8 indicates that this dependence becomes very weak in this limit.

In order to ensure an efficient computation of the Eddington factor in this domain, I therefore chose to use an 11th-order polynomial, neglecting the residual dependence of $\chi_g$ on the radiative temperature. The retained expression is the following:

$$ \chi_{g,L}(f_g) = 1/3 + 2f_g^2/3 + \begin{aligned}[t] & 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 (3.6) $$

This polynomial expression was constructed so as to satisfy the following conditions: $\chi_{g,L}(0) = 1/3$, $\chi_{g,L}(1) = 1$, and $\chi'_{g,L}(0) = 0$. The bounds of domain L were chosen in order to guarantee a relative error below $10^{-3}$ % in the estimation of $\chi_g$. The maximum error on the components $P_{g,xx}$ and $P_{g,yy}$ of the radiative pressure tensor is of the order of $10^{-1}$ %, while that on the component $P_{g,xy}$ can reach up to 1 %.

Figure 3.10 compares the estimates obtained from this polynomial approximation (purple curve) with the values of $\chi_g$ computed using the line-search algorithm (green curve). An excellent agreement is observed between the two, confirming the validity of the approximation. It is also worth noting that the expression proposed by Minerbo (1978) [89] (red curve) provides values of the Eddington factor that are very close to those obtained numerically, with a maximum relative error of $10^{-3}$ % (see figure 3.10). This model therefore constitutes a relevant approximation in this specific case, although it cannot be exploited within the adopted
<figure><img src="image_33.png" /></figure>

Figure 3.10: Comparison of the computation of the Eddington factor between the polynomial expression $\chi_{g,L}$, defined by equation (3.6), and reference data obtained by solving equations (2.72), (2.73), and (2.74) using a line-search algorithm, for $\mathcal{T}_g = 10^{-8}$ and $\delta_g = 10^{-4}$. This comparison also includes results obtained by solving the equations proposed by Minerbo (1978) using a Newton algorithm.

approach, since it also requires a search algorithm.

## ❖ Intermediate domains I

For all intermediate domains I, except for domain $I_b$ (see figures 3.7 and 3.8), I used different neural networks specialized for each domain:

* **domain I_a:** neural network MLP_a with two inputs: $f_g$ and $\mathcal{T}_g^{\nu_g}$,

* **domain I_c:** neural network MLP_c with two inputs: $f_g$ and $\mathcal{T}_g^{\nu_{g+1}}$,

* **domain I_d:** I used the neural network MLP_a, which predicts accurate values of the Eddington factor in this domain,

* **domain I_e:** I use two neural networks, MLP_{e,1} and MLP_{e,2}, in order to reduce the required network size and, consequently, the prediction time. MLP_{e,1} is used for $\delta_g \in ]0.03 ; 0.1]$ and MLP_{e,2} for $\delta_g \in ]0.1 ; 0.9]$. These two neural networks have three inputs: $f_g$, $\delta_g$, and $\mathcal{T}_g^{\nu_{g+1}}$,

* **domain I_f:** neural network MLP_f with two inputs: $f_g$ and $\mathcal{T}_g^s$ (defined by equation (3.1)).
The remainder of this section is devoted to detailing the answers I provided to the questions listed in Section 3.1. *Steps to generate a functional neural network*, with the aim of obtaining reliable neural networks.

## Which data should be used?

<table><thead><tr><th></th><th></th><th>Training dataset</th><th>Validation dataset</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>v<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>v<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="4">MLP<sub>e,1</sub></td><td>log<sub>10</sub>(T<sub>g</sub><sup>v<sub>g+1</sub></sup>)</td><td>100</td><td>150</td><td>[-3; 0.6]</td></tr><tr><td>f<sub>g</sub></td><td>100</td><td>100*</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># data</td><td>590 160</td><td>1 467 561</td><td>-</td></tr><tr><td rowspan="4">MLP<sub>e,2</sub></td><td>log<sub>10</sub>(T<sub>g</sub><sup>v<sub>g+1</sub></sup>)</td><td>100</td><td>150</td><td>[-3; 0.6]</td></tr><tr><td>f<sub>g</sub></td><td>100</td><td>100*</td><td>[0; 1]</td></tr><tr><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 rowspan="3">MLP<sub>f</sub></td><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>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.2: Description of the datasets used for training and validation of the neural networks.</caption>

To train the neural networks, I generated training and validation datasets using the line-search algorithm that I developed (Appendix C. *Computing the Eddington factor using search algorithms*). Table 3.2 summarizes the parameters used to generate these datasets. A linear spacing was used for all quantities. I draw the reader's attention to the fact that the expected amount of data does not correspond to the amount of data actually obtained, since the line-search algorithm does not always converge (see Appendix C.2. *Comparison of the search*

<footnote>*The interval $f_g \in [10^{-5}; 0.99999]$ was used in the validation datasets of MLP<sub>e,1</sub> and MLP<sub>e,2</sub> to vary the values of f<sub>g</sub>.</footnote>