algorithms). Here, I used between 2 and 35 times more data in the validation set than in the training set in order to evaluate the performance of the neural networks at intermediate values and to assess their interpolation capabilities. Finally, to build a test set, a reference Marshak wave simulation using this method will be performed in Section 3.2.3. Performance of the proposed method, in order to evaluate both the accuracy of the method and the prediction error.

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

Given the variable impact of errors in the Eddington factor on the components of the radiative pressure tensor (see equation 3.2), I use a modified mean squared error loss function, giving priority to the regions where the Eddington factor has a stronger influence on the relative error of the radiative pressure components. To this end, I defined two loss functions:

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

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

where $\chi_{g,\text{data}}$ and $\chi_{g,\text{pred}}$ are, respectively, the Eddington factor known from the data and the one predicted by the neural network; $x_1$, $x_2$, and $x_3$ are the neural network inputs, whose expressions are given by equations (3.12)–(3.15); and $\alpha_f$, $\alpha_T$, and $\alpha_v$ are factors defined as:

$$ \alpha_T = \min(x_1, 1 - x_1) + \epsilon_T, \quad (3.9) $$

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

$$ \alpha_v = \min(x_3, 1 - x_3) + \epsilon_v, \quad (3.11) $$

where $\epsilon_T$, $\epsilon_f$, and $\epsilon_v$ are constants introduced to prevent these factors from vanishing. I set these constants to 0.1, $10^{-3}$, and 0.1, respectively. These factors serve two purposes:

* $\alpha_T$ and $\alpha_v$ increase the weight of regions located at the interfaces between the domains defined in table 3.1 and illustrated in figure 3.6. They therefore help mitigate the discontinuities induced by the use of multiple specialized neural networks;

* $\alpha_f$ increases the importance of regions where an error in the Eddington factor induces a significant relative error in the radiative pressure (see equation 3.2).

The loss function $\mathcal{L}_1$ is used to train the neural networks MLP$_a$, MLP$_c$, and MLP$_f$, while $\mathcal{L}_2$ is used for MLP$_{e,1}$ and MLP$_{e,2}$.
## Which optimization algorithm should be used?

To train the neural networks, I used the Flux library in Julia, which is well known for its versatility in neural network modeling. The L-BFGS optimizer, provided by Julia's Optim library, is employed because of its fast convergence and its efficient minimization of the loss function. Although L-BFGS does not provide the stochastic properties that are often beneficial for promoting generalization, comparisons performed with the classical Adam optimizer show that networks trained with L-BFGS retain good generalization capabilities in the scenarios considered. I will discuss the generalization performance in more detail in the section *Should a specific strategy be adopted to promote generalization?* The models are trained until the neural network parameters exhibit variations smaller than $10^{-32}$, that is, until the predicted values vary by less than $10^{-32}$, or until the gradient norm falls below $10^{-8}$.

## Which neural network architecture should be used?

All these neural networks do not take the quantities mentioned above directly as inputs, but rather normalized quantities, in order to help them avoid giving too much importance to some inputs relative to others.

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

$$ \left\{ \begin{aligned} x_1 &= \frac{\log_{10}(\mathcal{T}_g^{\nu_{g-1/2}}) - \log_{10}(\mathcal{T}_{g,\min}^{\nu_{g-1/2}})}{\log_{10}(\mathcal{T}_{g,\max}^{\nu_{g-1/2}}) - \log_{10}(\mathcal{T}_{g,\min}^{\nu_{g-1/2}})}, \\ x_2 &= f_g . \end{aligned} \right. \qquad (3.12) $$

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

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

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

$$ \left\{ \begin{aligned} x_1 &= \frac{\log_{10}(\mathcal{T}_g^{\nu_{g+1/2}}) - \log_{10}(\mathcal{T}_{g,\min}^{\nu_{g+1/2}})}{\log_{10}(\mathcal{T}_{g,\max}^{\nu_{g+1/2}}) - \log_{10}(\mathcal{T}_{g,\min}^{\nu_{g+1/2}})}, \\ x_2 &= f_g , \\ x_3 &= \frac{\delta_g - \delta_{g,\min}}{\delta_{g,\max} - \delta_{g,\min}} . \end{aligned} \right. \qquad (3.14) $$
<figure><img src="image_34.png" /><figcaption>Figure 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$.</figcaption></figure>

* Inputs for MLP$_f$:

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

The quantities $X_{min}$ and $X_{max}$ represent the minimum and maximum values over the considered domain, and $x_i$ denotes an input of the neural networks.

Moreover, these neural networks do not directly predict the Eddington factor $\chi_g$, but rather a quantity related to it, according to the following expression:

$$ \chi_g = \frac{1}{3} + \frac{2}{3} f_g^2 \left( 1 + \frac{(1 - f_g)(3y - 2)}{2} \right). \qquad (3.16) $$
<figure><img src="image_35.png" /><figcaption>Figure 3.12: Learning curves of the different architectures tested for MLP<sub>a</sub>. *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.</figcaption></figure>

Here, y denotes the output of the neural network, which has been empirically verified to always lie within the interval [0; 1] in the data. This construction ensures that the predicted Eddington factor satisfies the expected physical conditions: it tends toward 1/3 when f<sub>g</sub> = 0, toward 1 when f<sub>g</sub> = 1, and exhibits a zero slope at f<sub>g</sub> = 0, in accordance with its expected behavior (see figure 2.4 and the definition recalled in Section 3.2. *Eddington factor in the M1-multigroup model*). Consequently, the network output y primarily affects the intermediate values of the Eddington factor, where its dependence on f<sub>g</sub> is the strongest.

For all these neural networks, I considered various activation functions for the hidden and output layers and found that they did not have a significant impact. Therefore, I ultimately used hyperbolic tangent activation functions for the hidden layers, a linear activation function for the input layer, and a sigmoid activation function for the output layer (see the summary figure 3.11). Only two hyperparameters remain to be tuned: the number of hidden layers and the number of neurons per hidden layer. These two parameters are selected via best-performing model selection (see table 3.3 for the final adopted meta-parameters).

For each neural network MLP<sub>a</sub>, MLP<sub>b</sub>, etc., several architectures were generated by varying the number of hidden layers and the number of neurons per hidden layer, while keeping the models sufficiently simple in order to obtain fast predictions (less then 150 parameters to train). Training was performed on the training set specified in table 3.2. To evaluate the performance of each tested architecture, learning curves were plotted, showing the evolution of the cost
function computed on the training set over the epochs. These curves provide two key pieces of information:

1. They make it possible to verify that training has been fully completed (absence of underfitting), which is reflected by a plateau at the end of training;

2. They allow identification of the architecture that reaches the minimum value of the cost function at the end of training.

In the case of the neural network MLP<sub>a</sub> (see figure 3.12), all tested architectures exhibit complete training, and the one consisting of three hidden layers with 15 neurons each proves to be the most accurate, reaching the lowest value of the cost function at convergence. The same methodology was applied to select the architectures of the other neural networks. The final selected hyperparameters are summarized in table 3.3.

<table><caption>Table 3.3: Architecture of the neural networks used.</caption><thead><tr><td>Neural network</td><td>Total # of layers</td><td>Total # of neurons</td><td># of parameters</td></tr></thead><tbody><tr><td>MLP<sub>a</sub></td><td>3</td><td>15</td><td>77</td></tr><tr><td>MLP<sub>c</sub></td><td>3</td><td>15</td><td>77</td></tr><tr><td>MLP<sub>e,1</sub></td><td>4</td><td>24</td><td>146</td></tr><tr><td>MLP<sub>e,2</sub></td><td>3</td><td>19</td><td>116</td></tr><tr><td>MLP<sub>f</sub></td><td>3</td><td>15</td><td>77</td></tr></tbody></table>

### Should a specific strategy be adopted to promote generalization?

To evaluate the generalization capability of the selected neural networks, I performed a 10-fold cross-validation. This procedure consists in dividing the chosen dataset (here, the data from the training set, as described in table 3.2) into 10 subsets of equal size. At each iteration, one of these subsets is used as the validation set, while the remaining nine are used for training. This process is repeated 10 times, so that each subset is used once as the validation set. For each tested architecture, I computed the mean and the standard deviation of a representative performance score in order to assess the robustness of the generalization.

The selected score is the mean absolute relative error on the estimation of the Eddington factor, defined as $|\Delta\chi_g|/\chi_{g,\text{true}}$. The results of this cross-validation for the network MLP<sub>a</sub> are presented in table 3.4. On average, the network exhibits a good generalization capability, as evidenced by a mean score on the validation sets comparable to that observed on the training
<table><caption>Table 3.4: Results of the 10-fold cross-validation for MLP<sub>a</sub>.</caption><thead><tr><th>Block number</th><th>Training set score</th><th>Validation set score</th></tr></thead><tbody><tr><td>1</td><td>1.2 × 10<sup>-4</sup></td><td>2.7 × 10<sup>-4</sup></td></tr><tr><td>2</td><td>4.9 × 10<sup>-5</sup></td><td>4.8 × 10<sup>-5</sup></td></tr><tr><td>3</td><td>3.9 × 10<sup>-3</sup></td><td>1.1 × 10<sup>-3</sup></td></tr><tr><td>4</td><td>1.5 × 10<sup>-4</sup></td><td>1.9 × 10<sup>-4</sup></td></tr><tr><td>5</td><td>2.1 × 10<sup>-4</sup></td><td>4.7 × 10<sup>-4</sup></td></tr><tr><td>6</td><td>2.1 × 10<sup>-4</sup></td><td>5.6 × 10<sup>-4</sup></td></tr><tr><td>7</td><td>8.0 × 10<sup>-5</sup></td><td>5.3 × 10<sup>-4</sup></td></tr><tr><td>8</td><td>6.8 × 10<sup>-4</sup></td><td>1.0 × 10<sup>-3</sup></td></tr><tr><td>9</td><td>2.1 × 10<sup>-2</sup></td><td>1.1 × 10<sup>-2</sup></td></tr><tr><td>10</td><td>8.1 × 10<sup>-5</sup></td><td>2.4 × 10<sup>-4</sup></td></tr><tr><td>Mean score</td><td>2.7 × 10<sup>-3</sup></td><td>1.6 × 10<sup>-3</sup></td></tr><tr><td>Standard deviation</td><td>6.7 × 10<sup>-3</sup></td><td>3.4 × 10<sup>-3</sup></td></tr></tbody></table>

<figure><img src="image_36.png" /><figcaption>Figure 3.13: Comparison of the distribution of the relative error of the MLP<sub>a</sub> prediction for the training set and the validation set.</figcaption></figure>

sets, and below 0.01 %, which is satisfactory. However, some folds (notably folds 7, 8, and 9) show significant discrepancies, indicating either localized overfitting or a difficulty of the model in adapting to certain specific sub-distributions of the data. In addition, the relatively large standard deviation observed in the validation scores reflects a non-negligible variability in performance across the different data splits, highlighting the sensitivity of the model to the distribution of the training set, but which generally leads to a neural network with good generalization properties. Similar observations were made for the other neural networks, with the exception of the network MLP<sub>c</sub>, for which the addition of L2 regularization to the loss function was necessary to ensure good generalization.

To more finely assess the interpolation performance of the final network, I analyzed its behavior on the training and validation sets described in table 3.2. The validation set, approximately twenty times larger than the training set, makes it possible to test the network on
intermediate points that are absent from the training data; these points represent the majority of cases encountered in practice. I examined the distributions of the relative errors on both sets for each trained network in order to verify the absence of overfitting.

In the case of the network MLP<sub>a</sub> (see figure 3.13), the distributions are very similar and centered around 0, which indicates a good generalization capability of the final neural network. Similar results were observed for the other architectures. Finally, the maximum errors measured on the validation sets remain sufficiently small for the requirements of the simulations.

### 3.2.3 Performance of the proposed method

In order to evaluate the performance of this new AI-based method, I carried out a classical one-dimensional Marshak wave simulation using the HADES code, in which a radiative flux is emitted from a heated region (see figure 3.14). Since the code is strictly two-dimensional, the configuration was set up in plane-parallel geometry in order to reproduce an effectively one-dimensional case.

We then compare the simulation time as well as the error on the radiative pressure obtained using four different approaches for computing the Eddington factor:

<table><caption>Table 3.5: Definition of the groups</caption><thead><tr><th>Number of groups</th><th></th><th>ν<sub>g-1/2</sub> [eV]</th><th>ν<sub>g+1/2</sub> [eV]</th><th>δ<sub>g</sub></th></tr></thead><tbody><tr><td rowspan="2">2 groups</td><td>- group 1</td><td>0</td><td>1.0 × 10<sup>-1</sup></td><td>-</td></tr><tr><td>- group 2</td><td>1.0 × 10<sup>-1</sup></td><td>+∞</td><td>-</td></tr><tr><td rowspan="3">3 groups</td><td>- group 1</td><td>0</td><td>1.0 × 10<sup>-3</sup></td><td>-</td></tr><tr><td>- group 2</td><td>1.0 × 10<sup>-3</sup></td><td>1.0 × 10<sup>-1</sup></td><td>1.0 × 10<sup>-2</sup></td></tr><tr><td>- group 3</td><td>1.0 × 10<sup>-1</sup></td><td>+∞</td><td>-</td></tr><tr><td rowspan="4">4 groups</td><td>- group 1</td><td>0</td><td>1.0 × 10<sup>-3</sup></td><td>-</td></tr><tr><td>- group 2</td><td>1.0 × 10<sup>-3</sup></td><td>2.0 × 10<sup>-1</sup></td><td>5.0 × 10<sup>-2</sup></td></tr><tr><td>- group 3</td><td>2.0 × 10<sup>-1</sup></td><td>4.0 × 10<sup>-1</sup></td><td>5.0 × 10<sup>-1</sup></td></tr><tr><td>- group 4</td><td>4.0 × 10<sup>-1</sup></td><td>+∞</td><td>-</td></tr><tr><td rowspan="5">5 groups</td><td>- group 1</td><td>0</td><td>1.0 × 10<sup>-3</sup></td><td>-</td></tr><tr><td>- group 2</td><td>1.0 × 10<sup>-3</sup></td><td>1.0 × 10<sup>0</sup></td><td>1.0 × 10<sup>-3</sup></td></tr><tr><td>- group 3</td><td>1.0 × 10<sup>0</sup></td><td>1.1 × 10<sup>0</sup></td><td>9.1 × 10<sup>-2</sup></td></tr><tr><td>- group 4</td><td>1.1 × 10<sup>0</sup></td><td>1.0 × 10<sup>1</sup></td><td>1.1 × 10<sup>-2</sup></td></tr><tr><td>- group 5</td><td>1.0 × 10<sup>1</sup></td><td>+∞</td><td>-</td></tr></tbody></table>
<figure><img src="image_37.png" /><figcaption>Figure 3.14: Initial and boundary conditions for the Marshak wave simulation.</figcaption></figure>

1. The line search algorithm described in Appendix C. *Computing the Eddington factor using search algorithms*;

2. Interpolation of the Eddington factor from precomputed values [37];

3. The analytical expression of the Eddington factor for the M1-gray model (equation 2.71);

4. The AI-based method that I described in this section.

In this simulation, we consider a domain of size $L_x = 1$ m, discretized into 1 000 cells in the x direction, and of size $L_y = 3$ mm, discretized into 3 cells in the y direction. The domain contains a fluid initially at rest, with a temperature of $T = 300$ K and a density of $\rho = 8.2 \times 10^{-9}$ kg/m³. On its left boundary, the fluid is heated from an imposed temperature of $T_l = 5780$ K, accompanied by an incoming radiative flux governed by the following expression:

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

where $G$ denotes the total number of groups. At the initial time, the radiation is in equilibrium with the gas. The latter is characterized by an atomic mass $\mu m_H = 1$ u and an adiabatic index $\gamma = 5/3$. The simulation covers a physical duration of 13.3 ns, during which different group configurations, presented in table 3.5, are used. Figure 3.14 summarizes the initial and boundary conditions.

## ❖ Computational time
Table 3.6 presents a comprehensive comparison of the simulation times for the different multigroup configurations described in Section 3.2.1. “Historical” methods. The AI-based method demonstrates remarkable efficiency, outperforming the line-search algorithm by a factor of 1 000 to 3 000 in terms of computational time. This approach is, however, slightly slower than the method using the analytical expression of the M1-gray model, which remains the fastest reference. Nevertheless, its computational cost remains comparable, exceeding that of the analytical M1-gray expression by only a factor of 1 to 3.

<table><caption>Table 3.6: Comparison of simulation times (CPU time).</caption><thead><tr><td>Number of groups</td><td>M1-gray model expression</td><td>Interpolation</td><td>AI method</td><td>Line-search algorithm</td></tr></thead><tbody><tr><td>2 groups</td><td>3.11 min</td><td>3.70 min</td><td>4.18 min</td><td>3.58 hours</td></tr><tr><td>3 groups</td><td>3.92 min</td><td>4.69 min</td><td>4.92 min</td><td>7.05 hours</td></tr><tr><td>4 groups</td><td>4.55 min</td><td>5.65 min</td><td>10.37 min</td><td>22.6 hours</td></tr><tr><td>5 groups</td><td>5.53 min</td><td>6.88 min</td><td>7.72 min</td><td>16.6 hours</td></tr></tbody></table>

It should be noted that the simulation times are higher for the 4-group configuration than for the others. This additional cost is mainly due to the use of the neural networks MLP<sub>e,1</sub> and MLP<sub>e,2</sub>, whose more complex architecture (see table 3.3) leads to a longer prediction time. As a result, the method I developed proves to be particularly efficient in terms of speed for a small number of groups, a situation in which the use of the neural networks MLP<sub>e,1</sub> and MLP<sub>e,2</sub> is not required.

### ❖ Precision of the method

In figure 3.15, the profile of the total radiative pressure, defined as $P_{R,xx} = \sum_g P_{g,xx}$, is shown as a function of the spatial position $x$ at time $t = 0.133$ ns, as obtained using the four methods for computing the Eddington factor (top panels). The figure also displays the relative error on the radiative pressure, computed with respect to that obtained using the line-search algorithm, for simulations employing the analytical expression of the M1-gray model, interpolation, and the AI method developed here (bottom panels).

The analysis of the curves shows that the AI-based method I developed exhibits a relative error below $10^{-3}$ % on the Eddington factor. In comparison, the use of the M1-gray model expression can lead to errors of up to 9 %, notably in the five-group simulation at the position $x = 25$ mm, where the reduced flux reaches its maximum value. The error can reach 0.1 %
<figure><img src="image_38.png" /><figcaption>Figure 3.15: Component $P_{R,xx}$ of the total radiative pressure at time $t = 0.133$ ns. The relative error is given by the following formula: $(P_{xx,method} - P_{xx,line search})/P_{xx,line search} \times 100$.</figcaption></figure>

in the case of the Eddington factor interpolation method, in particular for the five-group configuration at the same location.

These results highlight the significant improvement in accuracy provided by this new approach. Moreover, one can observe that the error associated with the M1-gray model expression as well as with the interpolation method tends to increase with the number of groups, unlike the AI-based method, whose accuracy remains stable.

❖ Conclusions on the method

In conclusion, I have developed a neural-network-based approach capable of predicting the Eddington factor within the M1-multigroup model with a mean error between $10^{-3}$ % and $10^{-2}$ %. This method is approximately 3 000 times faster than the line search algorithm, while remaining only 1 to 3 times slower than using the Eddington factor derived from the M1-gray model. An important advantage of this approach is that it does not require any prior knowledge of the radiative energy or the reduced flux, unlike methods based on the interpolation of precomputed Eddington factors, making it one of the most efficient solutions currently available.
However, the simulations carried out so far, including those with a large number of groups or in two dimensions, have not revealed significant differences in the hydrodynamic quantities between the different closure methods. This can be explained by the fact that simulations in which radiative pressure dominates, outside radiative equilibrium, remain too computationally expensive with the HADES code. Nevertheless, this new method constitutes a reference for future work relying on the M1-multigroup model, as it enables an accurate and efficient estimation of the Eddington factor, thereby paving the way for simulations with an unprecedented level of accuracy. In the following chapter, I will apply this method to the study of radiative shocks in order to concretely assess its impact on the structure and dynamics of these phenomena.

## 3.3 Synthesis

In this chapter, I first briefly introduced neural networks, presenting their operating principles, the fundamentals of training, as well as the key steps involved in their design and selection. I then detailed the method I developed during my thesis, based on AI, and aimed at estimating the Eddington factor of the M1-multigroup model in a general manner.

This approach achieves a mean prediction error between $10^{-3}$ % and $10^{-2}$ %, while being approximately 3 000 times faster than a line search algorithm, and only 1 to 3 times slower than using the Eddington factor from the M1-gray model. It also has the advantage of requiring no prior knowledge of the radiative energy or the reduced flux, unlike methods based on the interpolation of precomputed data, which makes it one of the most efficient solutions currently available. However, the simulations conducted so far have not revealed any notable differences in the hydrodynamic quantities, regardless of the number of groups or the dimensionality of the simulation (1D or 2D). This lack of discrepancy can be explained by the current inability to explore configurations out of radiative equilibrium dominated by radiative pressure, whose simulation remains particularly computationally expensive with the HADES code. Despite this, the proposed method constitutes a solid reference for future work relying on the M1-multigroup formalism, by providing an accurate and fast estimation of the Eddington factor.
<figure><img src="image_39.png" /></figure>
STUDY OF RADIATIVE SHOCKS

Radiative shocks constitute a particular category of shock waves in a fluid, in which radiation plays a central role in the dynamics and structure of the shock. Unlike purely hydrodynamic shocks, where energy and momentum are essentially conserved within the fluid, radiative shocks convert a significant portion of them into radiation. Radiative energy transport profoundly influences their behavior: the energy deposited in the shocked region can excite the rotational-vibrational states of molecules [125], and may even ionize atoms in extreme cases. This intense radiation emission in turn leads to additional heating and ionization of the surrounding medium, establishing a complex coupling between fluid dynamics, radiative transfer, and thermodynamic processes. These interactions significantly modify the structure of the shock, affecting the temperature, pressure, and density profiles [95, 126, 127]. Radiative shocks are ubiquitous in many astrophysical environments, making their study crucial for understanding a wide range of phenomena. They play a key role in accretion and ejection processes, particularly during gravitational collapse or supernovae, where high fluid velocities generate shocks accompanied by strong radiation. This radiation is often difficult to observe directly, which makes numerical studies essential for their interpretation. Such shocks occur in contexts as diverse as supernova remnants [128, 10], stellar atmospheres [129, 7, 8, 130], star formation through accretion processes [131, 76, 75], astrophysical jets [132, 133, 134], stellar mass-loss mechanisms [135, 136], and symbiotic stars [137, 74].

Their study relies simultaneously on analytical, experimental, and numerical approaches. Analytical models have been proposed to describe characteristic structures of the shock, such as the cooling layer or the radiative precursor [138, 139, 140, 141]. However, these models rest
on simplifying assumptions intended to handle the complexity of the underlying equations, and therefore cannot capture the full impact and interplay of all the physics involved. Advances in high-energy laser facilities have made it possible to reproduce radiative shocks experimentally in the laboratory [142, 143, 144, 145]. Nevertheless, these experiments remain limited due to the presence of instabilities, observed experimentally [146] and interpreted theoretically as axial [147] or lateral instabilities [148]. In this context, numerical simulations play a central role by providing predictive tools that are invaluable for interpreting observations and exploring fundamental physical mechanisms.

Thanks to the improvements made in the HADES code concerning the calculation of the Eddington factor, I undertook a more detailed analysis of the influence of an accurate treatment of radiation on the structure and dynamics of radiative shocks. In this chapter, I will begin by presenting the current state of knowledge on radiative shocks, before discussing the simulations I have carried out in order to assess the impact of accurate radiation modeling on the evolution and internal structure of radiative shocks.

## 4.1 Introduction to radiative shocks

The behavior of radiative shocks depends strongly on the interaction between radiation and gas, an interaction that varies with the opacity of the medium. In optically thick media, this interaction is dominant: radiation efficiently exchanges energy and momentum with matter, thereby deeply modifying the shock dynamics. Conversely, in optically thin media, its influence becomes negligible and the shock behaves in an essentially hydrodynamic manner. Five parameters can be used to characterize one-dimensional radiative shocks [149]:

### 1. Mach number M:

The Mach number characterizes the ratio between the velocity of a fluid and the speed of sound in that same medium. In the context of shock waves, it quantifies the intensity of the shock by comparing the relative velocity between the fluid and the shock front with the propagation speed of acoustic waves. It is written:

$$M = \frac{|v_{up} - v_{shock}|}{c_s}, \qquad (4.1)$$

where $v_{up}$ is the fluid velocity upstream of the shock, $v_{shock}$ is the shock-front velocity, and $c_s$ is the local upstream sound speed. In an ideal gas, one has $c_s = \sqrt{\gamma p_{up}/\rho_{up}}$. The Mach number determines the strength of the shock: the larger M is, the stronger the shock. By definition, a shock can occur only if $M > 1$.
## 2. Cooling parameter $\chi$:

This parameter compares the ability of the medium to lose its thermal energy through radiation to the dynamical evolution of the fluid. It is defined as the ratio between the radiative cooling time $t_{\text{cool}}$ and the characteristic dynamical time $t_{\text{dyn}}$ of the downstream fluid. In an optically thin medium:

$$ \chi = \frac{t_{\text{cool}}}{t_{\text{dyn}}} = \frac{1}{\gamma - 1} \frac{k_B}{\mu m_H} \frac{\rho_{\text{down}} v_{\text{down}} T_{\text{down}}}{\Lambda(\rho_{\text{down}}, T_{\text{down}}) L}, \quad (4.2) $$

and in an optically thick medium:

$$ \chi = \frac{1}{\gamma - 1} \frac{k_B}{\mu m_H} \frac{\rho_{\text{down}} v_{\text{down}} T_{\text{down}}}{F_{R,\text{down}}}, \quad (4.3) $$

where $\rho_{\text{down}}$, $v_{\text{down}}$, $T_{\text{down}}$, and $F_{R,\text{down}}$ denote respectively the density, velocity, temperature, and radiative flux of the downstream fluid in the shock frame, and $L = v_{\text{down}} t_{\text{dyn}}$ is a characteristic length. The function $\Lambda(\rho, T)$ represents the cooling function used in optically thin modeling (see equation 2.105). The constants $\gamma$, $k_B$, and $\mu m_H$ correspond respectively to the adiabatic index, Boltzmann's constant, and the mean molecular mass of the fluid. The smaller this parameter is, the more important the effects of radiative cooling are downstream.

## 3. Optical depth $\tau$ (upstream/downstream):

This quantity measures how transparent a medium is to radiation. It is defined as the integral of a mean opacity over the entire spectrum, $\kappa_{\text{ref}}$, along a path $s$:

$$ \tau = \int_{s_0}^{s} \kappa_{\text{ref}}(s') \, ds' . \quad (4.4) $$

One may consider either the Planck mean opacity $\kappa_P$ or the Rosseland mean opacity $\kappa_R$, depending on the radiative regime. For a homogeneous medium over a characteristic length $L$, this expression simplifies to:

$$ \tau = \kappa_{\text{ref}} L . \quad (4.5) $$

A medium with $\tau \ll 1$ is called optically thin, whereas a medium with $\tau \gg 1$ is optically thick.

## 4. Mihalas number $R$ (upstream/downstream):

The Mihalas number compares the internal energy of the gas $u$ to the radiation energy $E_R$. It is expressed as:

$$ R = \frac{u}{E_R} = \frac{1}{\gamma - 1} \frac{k_B}{a_R \mu m_H} \frac{\rho}{T^3 g(\tau)}, \quad (4.6) $$
where $\rho$, $v$, and $T$ are the hydrodynamic quantities upstream or downstream of the shock, $a_R$ is the radiation constant, and $g(\tau)$ is the ratio between the actual radiation energy $E_R$ and that of a blackbody $a_R T^4$. One has $g(\tau) \to 1$ when the medium is optically thick, and $g(\tau) \sim \tau$ in the optically thin regime. The smaller $R$ is, the more significant the radiative effects are.

**Boltzmann number Bo (upstream/downstream):**

The Boltzmann number is the ratio between the internal energy flux transported by the fluid and the energy flux transported by radiation. It is given by:

$$ Bo = \frac{(u + p)v}{F_R} = \frac{\gamma}{\gamma - 1} \frac{k_B}{\sigma \mu m_H} \frac{\rho v}{T^3 f(\tau)} , \quad (4.7) $$

where $\sigma$ is the Stefan-Boltzmann constant and $f(\tau)$ corresponds to the ratio between the actual radiative flux $F_R$ and the blackbody flux $\sigma T^4$. One has $f(\tau) \to 1$ when the medium is optically thick, and $f(\tau) \sim \tau$ in the optically thin regime. The smaller this parameter is, the more important radiative energy transport becomes.

Depending on the values taken by these five parameters, the following classes of one-dimensional radiative shocks can be distinguished [138, 140, 149]:

1. **Purely hydrodynamic shock:** Radiation is negligible ($\chi \to +\infty$, $Bo \to +\infty$, $R \to +\infty$, $\tau = 0$). The shock structure reduces to a discontinuity separating the upstream and downstream regions, and the Euler equations are sufficient to describe the situation. The temperature and density profiles in this case therefore display a discontinuity and are shown in figure 4.1a;

2. **Optically thin shock:** The medium is optically thin, which allows photons emitted downstream of the shock front to escape freely. This radiative emission then acts as an efficient mechanism for dissipating energy. As a result, a cooling region forms immediately downstream of the shock. This regime is typically characterized by $\chi \ll 1$, $Bo \gg 1$, and $R \gg 1$. The associated density and temperature profiles are illustrated in figure 4.1b. In this case, the loss of energy through photon emission leads to cooling of the fluid, manifested by a temperature peak followed by a rapid drop in the cooling region.

3. **Hybrid shock:** This type of shock occurs when the medium is optically thick on one side and optically thin on the other. In the optically thick region, radiation interacts strongly with the fluid ($\chi \ll 1$, $Bo \ll 1$), whereas it escapes freely in the optically thin region,
<figure><img src="image_40.png" /><figcaption>(a) Hydrodynamic shock.</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>(c) Optically thick shock, for Mach numbers
$M < 2.4 M_{rad}$.</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>(b) Optically thin shock.</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>(d) Optically thick shock, for Mach numbers
$M > 2.4 M_{rad}$.</figcaption></figure>

Figure 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.

leading to significant losses of energy and momentum. The density and temperature profiles are then similar to those observed in the optically thin case (see figure 4.1b), but the downstream temperature is much higher than the upstream one ($T_{down} \gg T_{up}$);

4. **Optically thick radiative shock dominated by radiative flux:** The medium is optically thick on both sides of the shock, and energy transfer is mainly carried by the radiative flux ($B_0 \ll 1, R \gg 1$). A *radiative precursor* forms upstream, resulting from fluid–radiation interactions, without any significant contribution from radiation energy or pressure. This type of shock may occur at relatively low Mach numbers, typically when $M < M_{rad} \approx 0.25 (\mu m_H)^{-1/6} \gamma^{-1/2} T_{amont}^{-1/2} \rho_{amont}^{1/6}$. The corresponding temperature and density profiles are shown in figure 4.1c. In this regime, photons interact strongly with matter and heat it upstream of the shock, creating a “bulge” in the temperature profile,
visible in the figure, corresponding to the radiative precursor;

5. **Optically thick radiative shock dominated by radiation:** This regime is similar to the previous one, but here it is the radiation energy and pressure that dominate the shock dynamics (Bo ≪ 1, R ≪ 1). This type of shock appears for Mach numbers greater than $M_{rad}$. When $M_{rad} < M < 2.4 M_{rad}$, the density profile exhibits a discontinuity, and the density and temperature profiles remain close to those observed in the flux-dominated case (see figure 4.1c). Conversely, for $M > 2.4 M_{rad}$, the density profile becomes continuous. The density and temperature profiles are then modified: as illustrated in figure 4.1d, the density varies continuously, while the temperature still shows a “bulge” upstream of the shock, a signature of the radiative precursor.

In the case of optically thick radiative shocks, a radiative precursor develops upstream of the shock front as a result of the intense radiation emitted by the compressed material. According to Paul Drake's book "High-Energy-Density Physics: Fundamentals, Inertial Fusion and Experimental Astrophysics" [127], this precursor exhibits a structure composed of two physically distinct regions: the diffusive precursor and the transmissive precursor. The diffusive precursor corresponds to the zone in which radiation interacts strongly with matter, leading to a state close to radiative equilibrium. Radiative energy transport there occurs mainly by diffusion, a mechanism favored by the high opacity of the medium, and the gas temperature increases progressively due to radiation absorption. Conversely, the transmissive precursor corresponds to a region in which the material is not in radiative equilibrium. Radiation propagates more rapidly than the ability of the matter to adjust thermally, which leads to thermal decoupling: the radiation temperature is already high, while the gas temperature remains close to the upstream value. This asymmetry results in a ratio $T_R/T$ significantly different from unity, indicating a nonequilibrium regime. Depending on the structure of the precursor and the intensity of the shock, two types of radiative shocks are distinguished: *subcritical shocks* and *supercritical shocks*.

* In a subcritical shock, the temperature $T_-$ just upstream of the front is lower than the downstream temperature $T_{down}$. The precursor is then essentially transmissive, and no region is in radiative equilibrium (see figure 4.2a). This type of shock occurs for Mach numbers below a critical value, denoted $M_{crit}$. To date, no reliable analytical expression allowing one to compute this critical value $M_{crit}$ has been established in the literature.

* In a supercritical shock, the temperature $T_-$ is equal to $T_{down}$, which indicates radiative heating strong enough to establish local radiative equilibrium. The radiative precursor
<figure><img src="image_44.png" /><figcaption>Figure 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.</figcaption></figure>

then consists of a diffusive precursor close to the shock front, followed by a more distant transmissive precursor (see figure 4.2b). These shocks occur for Mach numbers greater than $M_{crit}$.

### 4.2 Influence of the spectral nature of light

Despite existing classifications, radiative shocks remain imperfectly understood, largely because of the complexity of the physical interactions they involve, particularly when it comes to realistically modeling the couplings between radiation and matter. To date, no model is able to accurately predict the size and structure of the different regions of the radiative precursor, especially when two- or three-dimensional configurations are considered. Only approximations valid in one-dimensional cases exist for estimating the jump relations, and the evolution of the shock propagation speed is still poorly understood.

Furthermore, previous studies have highlighted the decisive influence of spectral effects on the shape of the precursor and on the cooling region [122]. However, such studies most often fall within the framework of a stationary shock, in which the simulation is initialized with a discontinuity whose upstream and downstream conditions are determined from approximate models. This approach therefore introduces, from the outset, prior assumptions about the jump relations, which may bias the analysis.

In addition, this type of configuration does not allow one to examine the actual dynamics of a propagating shock, particularly the effects induced by light-matter interactions, such as radiative braking. These approaches also frequently rely on numerical approximations intended to reduce computation time, which can impair the fidelity of opacity modeling and the accuracy
<figure><img src="image_45.png" /><figcaption>Figure 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 $ρ_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.</figcaption></figure>

of the closure relation within the M1-multigroup model.

In this context, I devoted part of my PhD work to studying in greater detail the impact of accurate spectral modeling of radiation on radiative shocks, by computing the closure relation of the M1-multigroup model with a level of precision never previously achieved. To do so, I used the method I developed (see Section 3.2), without introducing any simplifications on the opacities, and relying on a simulation capable of capturing both the dynamics and the structure of the shock.

This section summarizes the main results of this study, published in the article “*Impact of frequency-dependent radiation on the dynamics and structure of radiative shocks*” [150].

### 4.2.1 Simulation description

We performed simulations with the HADES code of the evolution of a one-dimensional radiative shock generated by a fluid moving toward a rigid wall. Figure 4.3 illustrates this configuration. The initial conditions are homogeneous for the hydrodynamic quantities. A symmetric boundary condition is imposed on the left boundary in order to represent the presence of the wall (zero velocity and zero radiative flux at the wall), while a zero-order extrapolation condition is applied at the right boundary (no imposed condition). At later times, the interaction between the incoming fluid and the wall generates a shock front that progressively propagates away from the wall. This type of configuration has been widely used to study optically thin radiative shocks [147, 151], and constitutes a reference case for
validating radiative-hydrodynamics codes [152].

Two shock configurations were considered: “Mach 4” and “Mach 8”. The initial hydrodynamic parameters are given in table 4.1. For the “Mach 4” simulation, a 10 cm domain was discretized into 2 000 cells, whereas the “Mach 8” simulation was carried out on a 20 cm domain with the same number of cells.

Because HADES is a two-dimensional code, a second direction must be defined. Its size was set to 0.15 mm for the “Mach 4” configuration and to 0.3 mm for the “Mach 8” configuration, and it was discretized on 3 cells in both cases. Periodic boundary conditions were imposed in this direction, thus allowing a plane-parallel simulation. The calculations were continued until the shock reached a steady state.

<table><caption>Table 4.1: Initial hydrodynamic quantities.</caption><thead><tr><th>Case</th><th>Density [kg/m<sup>3</sup>]</th><th>Velocity [m/s]</th><th>Temperature [K]</th><th>Temperature [eV]</th></tr></thead><tbody><tr><td>“Mach 4”</td><td>1</td><td>-5 622.82</td><td>11 604.52</td><td>1</td></tr><tr><td>“Mach 8”</td><td>1</td><td>-11 824.97</td><td>11 604.52</td><td>1</td></tr></tbody></table>

The medium considered in this study consists of pure argon ($\mu m_H = 39.948$ u). The effects of ionization on the mean atomic mass and on the adiabatic index were neglected. Although this simplification tends to overestimate the downstream temperature, it remains justified in the context of the objective of this study, which is to analyze the impact of the spectral behavior of radiation on the physics of radiative shocks.

In all simulations based on this model, a constant group narrowness $\delta_g^1$ was adopted for each spectral group. This choice makes it possible to quantitatively evaluate the influence of using narrower groups on the structure of the radiative shock, which would not be possible with groups of different sizes. Only the first and last groups have different narrownesses: the first is smaller, because the TOPS opacity tables did not contain data below 1.47 eV, while the last is wider in order to capture residual radiation at high frequencies, even though this radiation does not play a predominant role in the phenomena studied. Additional details regarding the group partitioning are presented in table 4.2. The following radiation models were studied:

1. M1-gray model;

2. M1-multigroup model with a constant group narrowness of $\delta_g = 0.13$ (4 groups);

<footnote n="1">¹The group narrowness $\delta_g$ is defined in Section 2.2.3. *The Eddington factor.*</footnote>
3. M1-multigroup model with a constant group narrowness of $\delta_g = 0.50$ (9 groups);

4. M1-multigroup model with a constant group narrowness of $\delta_g = 0.68$ (16 groups).

For the “Mach 8” simulation, we tested an additional frequency spacing:

5. M1-multigroup model with a constant group narrowness of $\delta_g = 0.78$ (23 groups).

<table><caption>Table 4.2: Frequency bounds of each group used in the M1-multigroup simulations.</caption><thead><tr><th>Simulation</th><th># groups</th><th>Group</th><th>Group frequency bounds</th></tr></thead><tbody><tr><td rowspan="3">δ<sub>g</sub> = 0.13</td><td rowspan="3">4</td><td>Group 1</td><td>ν ∈ [1.47 eV; 7.94 eV]</td></tr><tr><td>Groups 2-3</td><td>Group narrowness δ<sub>g</sub> = 0.13</td></tr><tr><td>Group 4</td><td>ν ∈ [500 eV; 16 keV]</td></tr><tr><td rowspan="3">δ<sub>g</sub> = 0.50</td><td rowspan="3">9</td><td>Group 1</td><td>ν ∈ [1.47eV; 3.98eV]</td></tr><tr><td>Groups 2-8</td><td>Group narrowness δ<sub>g</sub> = 0.50</td></tr><tr><td>Group 9</td><td>ν ∈ [500 eV; 16 keV]</td></tr><tr><td rowspan="3">δ<sub>g</sub> = 0.68</td><td rowspan="3">16</td><td>Group 1</td><td>ν ∈ [1.47eV; 2.17eV]</td></tr><tr><td>Groups 2-15</td><td>Group narrowness δ<sub>g</sub> = 0.68</td></tr><tr><td>Group 16</td><td>ν ∈ [500 eV; 16 keV]</td></tr><tr><td rowspan="3">δ<sub>g</sub> = 0.78</td><td rowspan="3">23</td><td>Group 1</td><td>ν ∈ [1.47 eV; 7.94 eV]</td></tr><tr><td>Groups 2-22</td><td>Group narrowness δ<sub>g</sub> = 0.78</td></tr><tr><td>Group 23</td><td>ν ∈ [300 eV; 16 keV]</td></tr></tbody></table>

This last simulation enabled us to further investigate the impact of narrower groups on the dynamics of the “Mach 8” shock. We used the TOPS opacity tables<fnref n="2" /> for argon, provided by Los Alamos National Laboratory. These tables supply the Rosseland and Planck mean opacities for the M1-gray and M1-multigroup models as functions of gas density and temperature, under the LTE assumption.

### 4.2.2 Results

In these simulations, the upstream and downstream media are optically thick (with an optical depth of $\tau \sim 10^4$ downstream and $\tau \sim 1$ to $10^2$ upstream, for both cases considered). From the simulated data, one can compute the three dimensionless parameters characterizing radiative shocks — the cooling parameter $\chi$, the Boltzmann number Bo, and the Mihalas parameter R — using equations (4.3), (4.7), and (4.6), in order to gain insight into the nature of the shocks under

<footnote n="2">²https://aphysics2.lanl.gov/apps/</footnote>
study. For the “Mach 4” shock, the values of these parameters are $\chi \approx 2\,000$, $Bo \approx 0.2$, and $R \approx 3\,000$, whereas for the “Mach 8” shock, the values are $\chi \approx 100$, $Bo \approx 0.01$, and $R \approx 200$. These results show that both configurations correspond to optically thick radiative shocks dominated by radiative flux, with the “Mach 8” case representing a particularly strong example. The simulations, whose final profiles of hydrodynamic and radiative quantities are shown in figures 4.5 and 4.6, demonstrate that a detailed description of the spectral nature of radiation, made possible by using narrow groups in the M1-multigroup model, significantly modifies three essential aspects of the shock: its propagation speed, the hydrodynamic properties of the downstream medium, and the structure of the radiative precursor.

### ❖ Influence on the shock speed

<figure><img src="image_46.png" /><figcaption>Figure 4.4: Evolution of the shock position as a function of simulated time.</figcaption></figure>

First, we determined the position of the shock by identifying it as the point where the density gradient reaches its maximum. From this estimate, we deduced the shock velocity by plotting, in figure 4.4, the evolution of its position as a function of physical time for the “Mach 4” and “Mach 8” configurations. One observes that, at the beginning of the simulations, the shock decelerates (smaller slope on the curves), and then, beyond a certain time, the position of the shock front evolves linearly with respect to time. To avoid bias associated with this initial deceleration, we fitted a linear function to the last 50% of the time data. This regression makes it possible to estimate the final shock velocity, denoted $v_{shock}$.

We then plotted in figure 4.7 the evolution of the shock velocity $v_{shock}$ obtained as a function
<figure><img src="image_47.png" /><figcaption>Figure 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.</figcaption></figure>
<figure><img src="image_48.png" /><figcaption>(a) Gas temperature.</figcaption></figure>

<figure><img src="image_49.png" /><figcaption>(c) Radiative temperature.</figcaption></figure>

<figure><img src="image_50.png" /><figcaption>(e) Ratio of radiative temperature to gas temperature.</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>(b) Density.</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>(d) Velocity.</figcaption></figure>

<figure><img src="image_53.png" /><figcaption>(f) Radiative flux.</figcaption></figure>

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

Figure 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.
of the group narrowness $\delta_g$ in the different simulations. This analysis reveals that radiation affects the shocks differently depending on the configuration: in the “Mach 4” simulation, the shock is accelerated relative to the hydrodynamic case, whereas in the “Mach 8” simulation, it is decelerated. However, reducing the width of the groups in the M1-multigroup model (increasing the group narrowness $\delta_g$) systematically leads to a reduction in shock velocity in both cases, which results in a decrease in the Mach number, since the sound speed of the upstream medium does not change. The only exception occurs between the simulations with $\delta_g = 0.68$ and $\delta_g = 0.78$, where an acceleration of the shock is observed. However, the shock velocity in the M1-multigroup simulation with $\delta_g = 0.78$ remains lower than in all the other simulations.

<figure><img src="image_55.png" /><figcaption>Figure 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.</figcaption></figure>

❖ Influence on downstream hydrodynamic quantities

Since this framework provides the downstream hydrodynamic quantities without bias from the initial conditions, I compared the temperature and density obtained in the simulations with the predictions of an approximate analytical model based on the diffusion regime, which is applicable to optically thick radiative shocks [138]. This model is widely used to provide a first estimate of the downstream hydrodynamic quantities of a radiative shock, or to initialize stationary-shock simulations (see figure 4.8).

Because the shock velocity (and therefore the associated Mach number) decreases as the group width increases (that is, as $\delta_g$ decreases), this has a direct impact on the hydrodynamic quantities predicted by the analytical diffusion regime model, which explains the variations observed in the blue curves in figure 4.8.

However, despite these variations, one finds that this approximate model introduces devia-
<figure><img src="image_56.png" /><figcaption>(a) Downstream temperature in the “Mach 4” simulation.</figcaption></figure>

<figure><img src="image_57.png" /><figcaption>(c) Downstream temperature in the “Mach 8” simulation.</figcaption></figure>

<figure><img src="image_58.png" /><figcaption>(b) Downstream density in the “Mach 4” simulation.</figcaption></figure>

<figure><img src="image_59.png" /><figcaption>(d) Downstream density in the “Mach 8” simulation.</figcaption></figure>

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

Figure 4.8: Comparison of the downstream hydrodynamic quantities obtained from the simulations performed here with those predicted by the analytical model based on the diffusion regime approximation. The simulation with $\delta_g = 0$ corresponds to the simulation using the M1-gray model. The predicted value from the diffusion regime approximation changes because the shock velocity, and therefore the Mach number, evolves as a function of the group narrowness $\delta_g$.

tions ranging from 1 to 8% in the gas temperature and from 3 to 15% in the density, compared with the values obtained from the simulations. For the “Mach 4” configuration, the downstream medium is systematically hotter and less dense than predicted by the diffusion model, and the discrepancies decrease as the group narrowness $\delta_g$ increases. Conversely, in the “Mach 8” configuration, the downstream medium is systematically colder and denser, and the deviation with the approximate model increases as the group narrowness increases.

This highlights the challenges involved in accurately modeling the jump conditions of a radiative shock, and emphasizes the need for caution when relying on analytical expressions derived from the diffusion regime, particularly for simulation initialization or boundary-condition tuning.

### ❖ Influence on the radiative precursor
<figure><img src="image_61.png" /><figcaption>Figure 4.9: Size of the radiative precursor as a function of the group narrowness $\delta_g$. The simulation with $\delta_g = 0$ corresponds to the simulation using the M1-gray model.</figcaption></figure>

We define the diffusive precursor as the region upstream of the shock where the radiation and gas temperatures are approximately equal, that is, where $|T_R/T - 1| < 10^{-3}$. The transmissive precursor, located upstream of the diffusive precursor, is characterized by a significant difference between the radiation and gas temperatures, that is, $|T_R/T - 1| \ge 10^{-3}$ (see figures 4.5e and 4.6e). Figure 4.9 illustrates the evolution of the size of the different radiative precursors as a function of the spectral narrowness of the groups, characterized by the parameter $\delta_g$.

In both simulated cases, one observes that the transmissive precursor expands when radiative effects are modeled with greater accuracy, at the expense of the diffusive precursor and the upstream medium. This suggests that the M1-gray model overestimates the equilibrium between gas and radiation in the radiative precursor. Moreover, it is observed that the radiative flux in the transmissive precursor becomes positive, whereas it remains negative in other regions (see figures 4.5f and 4.6f). This behavior appears to be a general trend that is not well captured by the M1-gray model, particularly in the “Mach 8” simulation.

## ❖ Computational time

In order to study in more detail the influence of the spectral nature of radiation on radiative shocks, it is necessary to perform a larger number of simulations, varying in particular the shock strength (that is, the Mach number), the opacity, and the different classes of radiative shocks presented in Section 4.1. Although optimizing the computation of the Eddington factor in the M1-multigroup model has made such simulations feasible, they remain extremely expensive in terms of computational time, as shown in table 4.3. In this work, we parallelized the calculations, using up to 400 processors for the most demanding cases. Despite this, a
<table><caption>Table 4.3: Computational cost of the simulations in CPU seconds, hours, and days.</caption><thead><tr><th rowspan="2">Simulation</th><th colspan="2">"Mach 4"</th><th colspan="2">"Mach 8"</th></tr><tr><th></th><th>s</th><th></th><th>s</th></tr></thead><tbody><tr><td><strong>Hydrodynamics</strong></td><td>7.6</td><td>s</td><td>17</td><td>s</td></tr><tr><td><strong>M1-gray</strong></td><td>749.5</td><td>h</td><td>1 821.8</td><td>h</td></tr><tr><td></td><td>31.2</td><td>days</td><td>75.9</td><td>days</td></tr><tr><td><strong>M1-multigroup</strong></td><td>5 298.9</td><td>h</td><td>5 102.5</td><td>h</td></tr><tr><td>δ<sub>g</sub> = 0.14 - 4 groups</td><td>220.8</td><td>days</td><td>212.6</td><td>days</td></tr><tr><td><strong>M1-multigroup</strong></td><td>52 847.5</td><td>h</td><td>44 414.0</td><td>h</td></tr><tr><td>δ<sub>g</sub> = 0.48 - 9 groups</td><td>2 202.0</td><td>days</td><td>1 850.6</td><td>days</td></tr><tr><td><strong>M1-multigroup</strong></td><td>121 099.6</td><td>h</td><td>113 203.4</td><td>h</td></tr><tr><td>δ<sub>g</sub> = 0.68 - 16 groups</td><td>5 045.8</td><td>days</td><td>4 716.8</td><td>days</td></tr><tr><td><strong>M1-multigroup</strong></td><td>–</td><td></td><td>191 132.5</td><td>h</td></tr><tr><td>δ<sub>g</sub> = 0.78 - 23 groups</td><td>–</td><td></td><td>7 963.9</td><td>days</td></tr></tbody></table>

more systematic and extensive exploration of the spectral effects of radiation on the structure of radiative shocks remains out of reach for a code such as HADES. This limitation becomes even more pronounced when one is interested in multidimensional (2D/3D) radiative effects. This constraint stems from a fundamental limitation of the explicit schemes used in HADES, namely the CFL stability condition, which imposes a time step defined by:

$$ \Delta t \le \Delta x/v \qquad (4.8) $$

In purely hydrodynamic simulations, the velocity $v$ corresponds to the fluid velocity $|\vec{v}|$, whereas in radiative hydrodynamics it is replaced by the speed of light $c$, which is much larger ($c \gg |\vec{v}|$). As a consequence, the time step must be drastically reduced, which significantly increases the duration of the simulations, as illustrated in table 4.3. Thus, the simulation that describes radiative transfer with the greatest accuracy required 7 963.9 CPU days. After parallelization on 400 processors, it still lasted about 20 days, which remains particularly long for this type of simulation.

It therefore becomes essential to explore other approaches that allow these radiative shocks to be simulated more efficiently. In this perspective, we undertook to investigate the use of PINNs to perform radiative-hydrodynamics simulations. The following chapter presents my first attempts at shock simulation using this neural-network architecture.
## 4.3 Synthesis

In this chapter, I presented my simulation results on the structure of a radiative shock, modeling radiative transfer using the M1-multigroup model, which makes it possible to account finely for the spectral nature of light. This approach led to three main results:

1. The propagation speed of the shock decreases as the radiation spectrum is modeled with greater accuracy;

2. The predictions obtained from the diffusion regime approximation for the jump conditions show notable discrepancies: errors of 1–8% in the gas temperature and 3–15% in the density;

3. The more refined the spectral modeling, the larger the radiative precursor becomes, in particular the nonequilibrium region, which benefits the most from this extension.

These results highlight the crucial importance of accurate radiation-matter interaction modeling in simulations of radiative shocks. They also shed light on the limitations of the diffusion-regime approximation often used to initialize such simulations — an approximation that can introduce significant errors in the upstream conditions and lead to nonphysical propagation of the shock inside the computational domain.

However, simulations carried out with the HADES code remain particularly expensive in computational time, notably because of the Courant–Friedrichs–Lewy condition, which must be satisfied to ensure the stability of the explicit schemes in the code. This condition forces the use of a very small time step to maintain numerical stability when solving the radiative-hydrodynamics equations.
EXTRAPOLATION OF RADIATIVE SHOCK SIMULATIONS WITH
PHYSICS-INFORMED NEURAL NETWORKS

Among recent deep-learning approaches developed to reproduce or extrapolate physical simulations, some stand out for their ability either to explicitly incorporate physical constraints or to directly learn the operators governing the system dynamics. Physics-Informed Neural Networks (PINNs) [68] rely on the explicit incorporation of the governing equations into the loss function: the neural network is therefore trained to provide an approximate solution to these equations. Neural Operators (such as DeepONet or the Fourier Neural Operator) [110, 153, 154] adopt a different strategy by learning the integral operator that maps initial conditions and physical parameters to the solution at a later time. These methods offer remarkable generalization capability to new regimes or geometries, though at the cost of requiring accurate reference data from simulations for training. Finally, hybrid or reduced-order approaches (Hybrid Surrogates / ROM-AI) [155, 156, 157] consist in projecting simulated fields into a lower-dimensional space obtained through a decomposition method (such as Proper Orthogonal Decomposition (POD) or an autoencoder), and then training a neural network to reproduce the temporal evolution of these latent variables. This strategy makes it possible to reconstruct full physical fields at reduced cost while preserving the essential physical coherence and structure of the system. In this work, we chose to explore the PINN technique to extrapolate simulations of radiative shocks. Its formulation, directly based on the governing equations, has the advantage of not relying on any prior database of numerical simulations, making it a particularly suitable approach for our objective.
PINNs constitute an innovative AI-based approach for solving problems governed by PDEs. Introduced by Raissi et al. (2019) [68], it extends the pioneering work of Lagaris et al. (1998) [69]. Although this method can be used to solve inverse problems, the focus here is restricted to forward simulation, which is our primary objective. This approach has already been successfully applied to the modeling of shocks in pure hydrodynamics [158, 159, 160] and, more recently, to radiative transfer [161], but never yet to radiative hydrodynamics. In this chapter, we therefore investigate the potential of this new approach to facilitate and accelerate the modeling of radiative shocks.

## 5.1 Introduction to PINNs

The principle of PINNs is to approximate the solution of a PDE using a **neural network** (typically an MLP), whose training relies on minimizing a loss function composed of several terms: the **PDE residual**, the initial conditions, and the boundary conditions when they are specified. The spatial and temporal derivatives required to evaluate the residuals of the PDE are computed using **automatic differentiation** [70], a procedure based on the fact that any computer program, no matter how complex, can be decomposed into a sequence of elementary arithmetic operations and standard functions. By recursively applying the chain rule to this sequence, it is possible to obtain derivatives of arbitrary order with machine precision. Training of the neural network proceeds as long as the **total loss function** remains above a predefined threshold, or as long as its decrease remains significant. The operation of PINNs is summarized in figure 5.1.

PINNs are characterized by their great flexibility: they can be implemented with various types of neural network architectures, ranging from MLPs (see Section 3.1.1. *Multi-Layer Perceptrons (MLP)*) to more advanced structures such as Convolutional Neural Networks (CNN), Recurrent Neural Networks (RNN), LSTM, GAN, Bayesian Neural Networks (BNN), or even transformers [162]. For a detailed overview of their applications, the reader may refer to the review by Cuomo et al. (2022) [105].

In the context of this work, the focus is placed on the solution of hyperbolic PDE systems, which can be expressed in the following general form :

$$ egin{cases} 
\forall (t, \vec{x}) \in [t_0; t_f] \times \Omega, & \partial_t U(t, \vec{x}) + \vec{\nabla} \cdot F[U(t, \vec{x})] = S[U(t, \vec{x})], \\
\forall \vec{x} \in \Omega, & U(t_0, \vec{x}) = U_0(\vec{x}), \\
\forall (t, \vec{x}) \in [t_0; t_f] \times \partial\Omega, & B[U(t, \vec{x})] = B(t, \vec{x}).
\end{cases} 
\quad (5.1) $$

Here, $U$ denotes the function representing the state of the physical system, $\Omega$ is the spatial
<figure><img src="image_62.png" /><figcaption>Figure 5.1: Illustrative diagram of Physics-Informed Neural Networks.</figcaption></figure>

domain, and $[t_0; t_f]$ is the time interval. The terms $F[U]$ and $S[U]$ correspond respectively to the fluxes and to the source term (which may represent a gain or a loss) associated with the quantity $U$. The functions $U_0(\vec{x})$ and $B(t, \vec{x})$ denote respectively the initial conditions and the boundary conditions prescribed on the boundary of the domain $\partial\Omega$. Moreover, $\mathcal{B}$ represents an operator applied on the boundary function, which is typically a derivative for Neumann boundary conditions or the identity operator for Dirichlet conditions. In the case of the Euler equations (2.8), these terms take the following form:

$$ U = \begin{pmatrix} \rho \\ \rho \vec{v} \\ E \end{pmatrix} \quad ; \quad F[U] = \begin{pmatrix} \rho \vec{v} \\ \rho \vec{v} \otimes \vec{v} + p \mathbb{I}_d \\ (E + p) \vec{v} \end{pmatrix} \quad ; \quad S[U] = 0 . $$

In the case of radiative hydrodynamics, using the M1-gray model for the description of radiation, these terms become (2.100):

$$ U = \begin{pmatrix} \rho \\ \rho \vec{v} \\ E \\ E_R \\ c^{-2} \vec{F}_R \end{pmatrix} \quad ; \quad F[U] = \begin{pmatrix} \rho \vec{v} \\ \rho \vec{v} \otimes \vec{v} + p \mathbb{I}_d \\ (E + p) \vec{v} \\ \vec{F}_R \\ \mathbb{P}_R \end{pmatrix} \quad ; \quad S[U] = \begin{pmatrix} 0 \\ \vec{S} \\ cS^0 \\ -cS^0 \\ -\vec{S} \end{pmatrix} . $$
Unlike traditional numerical methods, which rely on an explicit discretization of space and time, PINNs adopt a so-called mesh-free approach: the solution is approximated continuously by a neural network, without an underlying mesh. This makes it possible to avoid certain classical limitations such as discretization errors, numerical stiffness, or stability constraints on the time step.

The originality of PINNs lies mainly in two key elements: the definition of their loss function and the strategy used to select training points in the space-time domain. These aspects will be discussed in greater detail in the following sections. For a detailed overview of methodological variants, limitations, and possible improvements, the reader may refer to the Master's thesis of Wagenaar (2023) [163].

### 5.1.1 Cost function

In PINNs, the loss function plays a central role. It guides the learning of the neural network by penalizing deviations between the solution approximated by the network and the constraints imposed by the physical problem — differential equations, initial conditions, or even boundary conditions when these are used.

For a system of PDEs of the form (5.1), we denote by $U_{\theta}(t, \vec{x})$ the approximation of the solution by a neural network parameterized by a vector $\theta$ (the weights and biases). The total loss function $\mathcal{L}$ can then be expressed as a weighted sum of several contributions:

$$ \mathcal{L} = \omega_{IC}\mathcal{L}_{IC} + \omega_{BC}\mathcal{L}_{BC} + \omega_{eq}\mathcal{L}_{eq}, \qquad (5.2) $$

where:

* $\mathcal{L}_{IC}$ is the loss on the initial conditions,

* $\mathcal{L}_{BC}$ is the loss on the boundary conditions,

* $\mathcal{L}_{eq}$ measures the residual of the differential equations (PDE),

* $\omega_{IC}$, $\omega_{BC}$, and $\omega_{eq}$ are weights adjusting the relative importance of each term.

In the literature, it is common to impose $\omega_{IC} \gg \omega_{eq}$ in order to force the network to properly satisfy the initial conditions before minimizing the residual of the physical equations [160]. A typical value is $\omega_{IC} = 100 \omega_{eq}$. As for $\omega_{BC}$, few authors specify its optimal value, but it is generally chosen such that $\omega_{IC} \ge \omega_{BC} > \omega_{eq}$.
## ❖ Initial conditions

The loss function $\mathcal{L}_{IC}$ quantifies the error between the prescribed initial values and those predicted by the network. Evaluated over a set of $N_{IC}$ points $(t_0, \vec{x}_i)$, it is written as:

$$ \mathcal{L}_{IC} = \frac{1}{N_{IC}} \sum_{i=1}^{N_{IC}} \{U_{\theta}(t_0, \vec{x}_i) - U_0(\vec{x}_i)\}^2, \quad (5.3) $$

where $U_0$ denotes the set of initial values of the function being sought.

## ❖ Boundary conditions

The formulation of the term $\mathcal{L}_{BC}$ depends on the type of boundary condition considered. Suppose that $N_{BC}$ points are used to evaluate the error associated with these conditions. For example, in the case of a Neumann condition imposing a zero normal derivative on the boundary, one may use:

$$ \mathcal{L}_{BC} = \frac{1}{N_{BC}} \sum_{i=1}^{N_{BC}} \partial_{\vec{x}_i} U_{\theta}(t_i, \vec{x}_i)^2, \quad (5.4) $$

where $\partial_{\vec{x}_i}$ is the derivative along the outward normal direction at the boundary. In the case of a Dirichlet condition imposing a specific value on the boundary, one may use the loss function:

$$ \mathcal{L}_{BC} = \frac{1}{N_{BC}} \sum_{i=1}^{N_{BC}} \{U_{\theta}(t_i, \vec{x}_i) - U_{\partial\Omega}(\vec{x}_i)\}^2, \quad (5.5) $$

where $U_{\partial\Omega}(\vec{x}_i)$ is the imposed value at the boundary point $\vec{x}_i$.

## ❖ Partial Differential Equations (PDE)

The term $\mathcal{L}_{eq}$ penalizes the failure of the neural network to satisfy the PDEs. Considering $N_{eq}$ points in the space-time domain, the standard formulation is:

$$ \mathcal{L}_{eq} = \frac{1}{N_{eq}} \sum_{i=1}^{N_{eq}} \{\partial_t U_{\theta,i} + \vec{\nabla} \cdot \mathbf{F}[U_{\theta,i}] - \mathbf{S}[U_{\theta,i}]\}^2, \quad (5.6) $$

where $U_{\theta,i} = U_{\theta}(t_i, \vec{x}_i)$ is the neural-network prediction at the space-time coordinates $(t_i, \vec{x}_i)$. This formulation enables a pointwise evaluation of the differential-equation residual.
<figure><img src="image_63.png" /><figcaption>Figure 5.2: Overview of different adaptations of the differential equation loss $\mathcal{L}_{eq}$ for addressing the shock problem.</figcaption></figure>

However, it shows its limitations in cases where the solution contains discontinuities, such as shocks in the Euler equations. Since standard activation functions are continuous, the solution learned by the network will also be continuous, making it difficult to approximate a discontinuous profile. Moreover, at discontinuities, the local PDEs are no longer strictly valid, as mentioned in Section 2.1. Hydrodynamics. This amounts to looking for an entropic solution, whereas a shock corresponds to a weak non-entropic solution. To overcome these difficulties, several strategies have been proposed, summarized in figure 5.2.

1. Weak formulation

A first approach consists of using an integral formulation of the PDEs [164, 165], which mitigates the effects of discontinuities:

$$ \mathcal{L}_{eq} = \frac{1}{N_{eq}} \sum_{i=1}^{N_{eq}} \left\{ \int_{\Omega} \int_{t_0}^{t_f} \left( \partial_t U_{\theta}(t, \vec{x}) + \vec{\nabla} \cdot \mathbf{F}[U_{\theta}(t, \vec{x})] - \mathbf{S}[U_{\theta}(t, \vec{x})] \right) \phi_i(t, \vec{x}) \, dt \, d\vec{x} \right\}^2, \quad (5.7) $$

where $\phi_i$ are test functions that vanish outside the domain. By taking $\phi_i$ as Dirac distributions, one recovers the classical formulation 5.6. This approach improves robustness to discontinuities but significantly increases computational cost.

2. Adaptive weighting of residuals

The impact of discontinuities can also be reduced by weighting the importance of the differential equations according to the gradient, using a factor $\lambda_i$:

$$ \mathcal{L}_{eq} = \frac{1}{N_{eq}} \sum_{i=1}^{N_{eq}} \lambda_i \left\{ \partial_t U_{\theta}(t_i, \vec{x})_i + \vec{\nabla} \cdot \mathbf{F}[U_{\theta,i}] - \mathbf{S}[U_{\theta,i}] \right\}^2, \quad (5.8) $$

Two possible choices for $\lambda_i$:
a) **Liu (2023) [166]:** reduction of the importance of compressible regions

$$ \lambda_i = \frac{1}{1 + \epsilon (|\vec{\nabla} \cdot \vec{v}| - |\vec{\nabla} \cdot \vec{v}|)}, \qquad (5.9) $$

where $\epsilon$ is a tunable hyperparameter and $\vec{v}$ the fluid's velocity,

b) **Ferrer-Sánchez (2024) [167]:** more general expression for the one-dimensional case:

$$ \lambda_i = \frac{1}{1 + \alpha_1 |\partial_x \rho|^{\eta_1} + \alpha_2 |\partial_x v|^{\eta_2} + \alpha_3 |\partial_x p|^{\eta_3}}, \qquad (5.10) $$

where $\alpha_1, \alpha_2, \alpha_3, \eta_1, \eta_2, \eta_3$ are hyperparameters to be tuned.

This approach is fast and significantly improves shock prediction, but requires careful tuning of the hyperparameters, and it may result to poor predictions of the shock dynamics.

la

### 3. Artificial viscosity

Another option consists in adding a diffusive term to the PDEs to smooth the discontinuity:

$$ \mathcal{L}_{eq} = \frac{1}{N_{eq}} \sum_{i=1}^{N_{eq}} \left\{ \partial_t U_{\theta,i} + \vec{\nabla} \cdot F[U_{\theta,i}] - S[U_{\theta,i}] - \nu(t_i, x_i) \Delta U_{\theta,i} \right\}^2, \qquad (5.11) $$

with $\nu$ a viscosity coefficient, which may be constant, parameterized, or learned [168]. This method requires computing second derivatives, which increases the training cost.

### 4. Relaxed networks

Finally, some approaches propose relaxing the strict definition of the fluxes [169], by predicting the fluxes separately using a second neural network:

$$ \mathcal{L}_{eq} = \mathcal{L}_{eq}^* + \epsilon_{dissip} \mathcal{L}_{dissip}, \qquad (5.12) $$

where $\epsilon_{dissip}$ is a tunable hyperparameter and:

$$ \mathcal{L}_{eq}^* = \frac{1}{N_{eq}} \sum_{i=1}^{N_{eq}} \left\{ \partial_t U_{\theta,i} + \vec{\nabla} \cdot F_{\theta,i} - S[U_{\theta,i}] \right\}^2, \qquad (5.13) $$

$$ \mathcal{L}_{dissip} = \frac{1}{N_{eq}} \sum_{i=1}^{N_{eq}} \{F_{\theta,i} - F[U_{\theta,i}]\}^2, \qquad (5.14) $$

where $F_{\theta,i}$ is predicted directly by a network. In some cases, this formulation improves the quality of the solution, at the price of increased computational cost and additional hyperparameters to tune.
### 5.1.2 Selecting data

Although PINNs provide a continuous solution in the space-time domain, it is nevertheless necessary to discretize this domain in order to evaluate the different components of the loss function. Three distinct datasets must be generated to train the neural network:

* one set for evaluating the initial conditions,

* one set for evaluating the boundary conditions,

* one set for evaluating the PDEs inside the domain.

The first two sets are generally created by random sampling, which is often sufficient to capture the essential information. By contrast, the selection of the points used for evaluating the PDEs requires particular attention.

As in classical numerical methods, where mesh quality strongly influences the accuracy of the solution, the choice and distribution of training points (in both space and time) have a crucial impact on the convergence and performance of PINNs. However, unlike traditional methods, this constraint is not as explicit: simply increasing the number of points does not necessarily improve the solution. Their distribution must be carefully designed according to the expected physical dynamics.

This stems from the formulation of the PINN loss function, which relies primarily on minimizing the residuals of the differential equations rather than explicitly comparing with a reference solution. This approach may lead to difficulties in correctly satisfying initial or boundary conditions and in ensuring proper propagation of information throughout the domain. In some cases, it may even violate the principle of causality, thereby disturbing the temporal propagation of the initial conditions [170].

To address these limitations, several sampling strategies have been developed in the literature. They can be grouped into three broad categories:

#### 1. Static non-adaptive strategies

These methods use a fixed set of training points, chosen before the start of training.

* **Grid sampling (figure 5.3a):** consists in regularly selecting spatio-temporal points $(t, x)$, similarly to a classical Cartesian mesh. This method is simple but generally yields the poorest results.
<figure><img src="image_64.png" /><figcaption>Figure 5.3: Examples of sampling strategies used in PINNS</figcaption></figure>

* **Uniform random sampling (figure 5.3b):** points are selected independently and randomly according to a uniform distribution. It slightly improves results compared to grid sampling, but carries a significant risk of undersampling in certain regions.

* **Latin hypercube sampling (figure 5.3c):** ensures better statistical coverage of the domain than uniform sampling, while retaining randomness. Although widely used in the literature [68, 171, 172], its performance remains comparable to that of uniform sampling.

* **Sobol sampling (figure 5.3d):** a quasi-random, low-discrepancy method that aims to minimize the variance of the integral estimator in the loss function (5.2). It provides a more homogeneous coverage of the space, accelerates convergence, improves accuracy, and ensures better representation of the global domain.

* **Hammersley sampling (figure 5.3e):** another quasi-random, low-discrepancy method, constructed from a Hammersley sequence. It offers a particularly uniform distribution of points in the training space, while remaining simple to generate. Its
coverage properties are comparable to those of the Sobol sequence and it is often used as an efficient alternative for reducing variance and improving the convergence of PINNs.

## 2. Non-adaptive resampling strategies

In this approach, the training points are periodically renewed (for example, every *n* epochs) using one of the previous schemes. This update helps avoid overfitting to a fixed configuration and significantly improves solution quality [172];

## 3. Adaptive strategies:

These methods dynamically modify the sampling according to the evolution of training or detected errors.

* **Residual-based sampling:** points are selected according to the local magnitude of the PDE residuals, emphasizing poorly resolved regions [173]. This method is effective but may under-represent certain regions of the domain.

* **Residual- and gradient-based sampling:** extends the previous method by also accounting for the local amplitude of solution gradients [174]. It is particularly useful for capturing strong-gradient phenomena such as shocks.

* **Probabilistic sampling:** instead of selecting only high-residual points, this approach samples points according to a probability distribution proportional to the local error [175]. It introduces beneficial diversity during training, maintaining better domain coverage while still favoring regions with large errors.

* **Neural-network-based sampling:** uses a generative model to learn an optimal sampling distribution, thereby improving the resolution of complex regions of the domain [176].

To assess the relevance of this approach and identify the methodological choices best suited to radiative-hydrodynamics simulations, I will, in the following section, present in detail the strategies adopted and the tests I performed.

## 5.2 Extrapolation of radiative shock simulations with Physics-Informed Neural Networks

In this section, I present the preliminary experiments carried out to evaluate the relevance of PINN-based neural networks for simulating radiative-hydrodynamics phenomena. PINNs
<figure><img src="image_65.png" /><figcaption>Figure 5.4: Distribution of the data used to train the PINN neural network. The blue points correspond to data originating from the reference simulation, the orange points to the locations where boundary conditions are imposed, and the green points to the positions where the PDE residuals are evaluated. Time and space are rescaled according to the transformations defined in the sections devoted to the equations in Sections 5.2.1. *Hydrodynamic tests* and 5.2.2. *Radiative hydrodynamic tests*.</figcaption></figure>

remain an actively developed approach, particularly regarding the treatment of discontinuities such as those that arise in shock waves. This difficulty makes their application to radiative hydrodynamics especially challenging.

The work presented here is therefore exploratory and will require further investigation before it can be fully integrated into radiative-hydrodynamics simulation codes. I would like to thank Corentin Pecontal, an intern at the Observatoire de la Côte d'Azur, for his invaluable assistance in carrying out these experiments while I was in the process of writing this manuscript.

The strategy adopted here differs slightly from that presented in Section 5.2. *Extrapolation of radiative shock simulations with Physics-Informed Neural Networks*. Rather than being limited to a set of initial conditions specified at time $t_0$, we assume that a dataset from a reference simulation is available, covering a time interval from $t_0$ to $t_i$, sampled at $N_t$ regularly spaced output times. The objective is then to extrapolate the evolution of the system beyond $t_i$ using a PINN. To train the network, three datasets are used, corresponding to the sets shown in figure 5.4:

1. a dataset extracted from a reference simulation, acting as the starting point (replacing the classical initial conditions mentioned in Section 5.1. *Introduction to PINNs*). These
data cover the time interval $[t_0; t_i]$;

2. a dataset imposing the boundary conditions, selected over the time interval $[t_i; t_f]$ for the purely hydrodynamic test (figure 5.4a), and over $[t_0; t_f]$ for the radiative-hydrodynamics test (figure 5.4b);

3. a set of points on which the PDE residuals are evaluated, distributed over the time interval $[t_i; t_f]$ for the purely hydrodynamic test (figure 5.4a), and over $[t_0; t_f]$ for the radiative-hydrodynamics test (figure 5.4b).

We consider here a physical configuration similar to that described in Section 4.2.1. *Simulation description*, with the only difference that we study a shock characterized by a higher Mach number. This choice is motivated by our ultimate objective, which is to simulate stationary radiative shocks at very high Mach number, a study that could not be performed in Section 4.2.1. *Simulation description* due to computational cost constraints.

The system considered corresponds to a fluid initially characterized by a density $\rho$, a temperature $T$, and a velocity $v$ directed toward a wall, leading to the formation of a shock that propagates in the opposite direction over time. The fluid, composed of argon, is modeled as a monatomic ideal gas, with adiabatic index $\gamma = 5/3$ and mean molecular mass $\mu m_H = 39.948$ u. The initial conditions of the fluid are summarized in Table 5.1. In this purely hydrodynamic framework, the resulting shock corresponds to a Mach number of about $M \approx 16.1$. The only explicitly imposed boundary condition concerns the wall located at $x = 0$ m, where the fluid velocity is kept equal to zero.

<table><caption>Table 5.1: Initial conditions of the test case used to evaluate the PINNs.</caption><thead><tr><th>Density [kg/m³]</th><th>Velocity [m/s]</th><th>Temperature [K]</th><th>Temperature [eV]</th></tr></thead><tbody><tr><td>1</td><td>-50 000</td><td>50 000</td><td>4.3</td></tr></tbody></table>

### 5.2.1 Hydrodynamic tests

First, we performed a series of tests in a purely hydrodynamic regime, with the aim of determining the learning strategy of the neural network and identifying the factors that influence the quality of the extrapolation. This section successively presents the formulation of the equations used, the selected loss function, the data distribution, the network architecture, and finally an analysis of the influence of the learning parameters on the results obtained.
The study domain is defined over the spatial interval $x \in [x_{min}; x_{max}] = [0; 10]$ m and the temporal window $t \in [t_0; t_f] = [0; 2 \times 10^{-4}]$ s.

## ❖ The equations

First of all, we reformulate the Euler equations (2.8) in a way that is easier to handle with a PINN, by explicitly expanding the derivatives of density, temperature, and velocity:

$$ \left\{ \begin{aligned} \partial_t \rho + v \partial_x \rho + \rho \partial_x v &= 0 \\ \partial_t v + v \partial_x v + \mathcal{R} \left( T \frac{\partial_x \rho}{\rho} + \partial_x T \right) &= 0 , \\ \partial_t T + v \partial_x T + (\gamma - 1) T \partial_x v &= 0 \end{aligned} \right. \qquad (5.15) $$

where $\gamma$ is the adiabatic index and $\mathcal{R} = k_B/(\mu m_H)$ is the specific gas constant.

As is often the case in hydrodynamics simulations, a rescaling of the physical quantities is applied in order to reduce differences in magnitude. This normalization facilitates training of the neural network by allowing it to predict values close to unity. The residuals of equations (5.15) then read:

$$ \left\{ \begin{aligned} \mathrm{res}_{\rho}(\tilde{t}, \tilde{x}) &= \partial_{\tilde{t}} \tilde{\rho} + \tilde{v} \partial_{\tilde{x}} \tilde{\rho} + \tilde{\rho} \partial_{\tilde{x}} \tilde{v} \\ \mathrm{res}_{\nu}(\tilde{t}, \tilde{x}) &= \partial_{\tilde{t}} \tilde{v} + \tilde{v} \partial_{\tilde{x}} \tilde{v} + \tilde{\mathcal{R}} \left( \tilde{T} \frac{\partial_{\tilde{x}} \tilde{\rho}}{\tilde{\rho}} + \partial_{\tilde{x}} \tilde{T} \right) , \\ \mathrm{res}_{\mathcal{T}}(\tilde{t}, \tilde{x}) &= \partial_{\tilde{t}} \tilde{T} + \tilde{v} \partial_{\tilde{x}} \tilde{T} + (\gamma - 1) \tilde{T} \partial_{\tilde{x}} \tilde{v} \end{aligned} \right. \qquad (5.16) $$

where the quantities with a tilde are redimensionalized according to:

$$ \tilde{t} = t/\dot{t} , \quad \tilde{x} = x/(\dot{v}\tilde{t}) , \quad \tilde{\rho} = \rho/\dot{\rho} , \quad \tilde{v} = v/\dot{v} , \quad \tilde{T} = T/\dot{T} , \quad \tilde{\mathcal{R}} = \mathcal{R}\dot{T}/\dot{v}^2 , $$

The constants $\dot{t}$, $\dot{\rho}$, $\dot{v}$, and $\dot{T}$ denote the rescaling factors associated with time, density, velocity, and temperature, respectively. In the experiments presented here, we used the following values: $\dot{t} = 10^{-4}$, $\dot{\rho} = 1$, $\dot{v} = 10^4$, and $\dot{T} = 10^5$. For simplicity of notation, we will henceforth omit the tildes, and all quantities will be considered rescaled by default, unless stated otherwise.

## ❖ The cost function

The cost function used to train the neural network comprises three distinct components: (i) the PDE residuals, denoted $\mathcal{L}_{PDE}$, (ii) the error between the network predictions and the
simulation data, denoted $\mathcal{L}_{sim}$, and (iii) the error on the imposed boundary condition, denoted $\mathcal{L}_{BC}$.

The first component, $\mathcal{L}_{PDE}$, assesses the consistency of the predictions with the underlying physical equations. It is computed using the previously defined residuals, and incorporates the adaptive weighting strategy proposed by Liu (2023) [166], designed to mitigate the impact of the shock discontinuity during training. It is written as:

$$ \mathcal{L}_{PDE} = \frac{1}{N_{PDE}} \sum_{i=1}^{N_{PDE}} \left\{ \frac{\mathrm{res}_{\rho}(t_i, x_i)^2 + \mathrm{res}_{\nu}(t_i, x_i)^2 + \mathrm{res}_{\mathrm{T}}(t_i, x_i)^2}{1 + \lambda(|\partial_{\bar{x}}v_i| - \partial_{\bar{x}}v_i)} \right\}, $$

where $N_{PDE}$ is the number of residual evaluation points located at $(t_i, x_i)$, $\mathrm{res}_{\rho}$, $\mathrm{res}_{\nu}$, and $\mathrm{res}_{\mathrm{T}}$ denote the PDE residuals evaluated from the network outputs, and $\lambda$ is a hyperparameter that will be discussed later.

The second component, $\mathcal{L}_{sim}$, measures the discrepancy between the network predictions and the simulation data, using logarithmic or hyperbolic transformations adapted to each variable. It is written as:

$$ \mathcal{L}_{sim} = \frac{1}{N_{sim}} \sum_{i=1}^{N_{sim}} \left\{ \left[ \log(\rho_{pred}(t_i, x_i)) - \log(\rho_{sim}(t_i, x_i)) \right]^2 + \left[ \log(T_{pred}(t_i, x_i)) - \log(T_{sim}(t_i, x_i)) \right]^2 + \left[ \mathrm{asinh}(v_{pred}(t_i, x_i)) - \mathrm{asinh}(v_{sim}(t_i, x_i)) \right]^2 \right\}, $$

where $N_{sim}$ denotes the number of simulation points $(t_i, x_i)$, and where the functions log and asinh correspond respectively to the base-10 logarithm and the inverse hyperbolic sine. These transformations homogenize the relative influence of physical quantities, which often differ by several orders of magnitude, in the error computation.

The third component, $\mathcal{L}_{BC}$, enforces the boundary condition that the velocity must be zero at $x = 0$ m. It is defined as:

$$ \mathcal{L}_{BC} = \frac{1}{N_{BC}} \sum_{i=1}^{N_{BC}} \mathrm{asinh}(v_{pred}(t_i, 0))^2, $$

where $N_{BC}$ is the number of time instants $t_i$ at which the condition is evaluated. The use of the asinh function here ensures consistency with the formulation of $\mathcal{L}_{sim}$.

The total loss function combining these three contributions is:

$$ \mathcal{L} = \omega_{PDE} \mathcal{L}_{PDE} + \omega_{sim} \mathcal{L}_{sim} + \omega_{BC} \mathcal{L}_{BC}, \quad (5.17) $$
where the coefficients $\omega_{PDE}$, $\omega_{sim}$, and $\omega_{BC}$ control the relative importance of each term. In this study, we chose the values $\omega_{PDE} = 1$, $\omega_{sim} = 100$, and $\omega_{BC} = 10$, so as to initially orient training toward learning from the simulation data, then impose the boundary conditions, and finally progressively encourage the physical extrapolation of the model through the PDEs.

### ❖ The data

For each component of the loss function, a specific dataset is used (see figure 5.4a):

1. $\mathcal{L}_{sim}$, **simulation data (set of blue points in figure 5.4a):** it is evaluated over the entire spatial domain, but only over the time interval $[t_0; t_i]$, corresponding to the supervised learning phase. Sampling is regular in space, with 100 points uniformly distributed in space and $N_T$ snapshots uniformly distributed in time over $[t_0; t_i]$. The corresponding data come from the analytical solution of the problem, presented in Appendix D. *Hydrodynamic shock on a wall*;

2. $\mathcal{L}_{PDE}$, **residual points (set of green points in figure 5.4a):** this component measures the PDE residuals over the extrapolation region, i.e., for $t \in [t_i; t_f]$ and for all $x \in [x_{min}; x_{max}]$. A total of 10 000 $(t, x)$ points are randomly sampled in this domain using a Hammersley quasi-random sequence;

3. $\mathcal{L}_{BC}$, **boundary conditions (set of orange points in figure 5.4a):** this component enforces the boundary condition at $x = 0$ during the extrapolation phase. We used 100 time points, equally spaced over the interval $t \in [t_i; t_f]$, all located at $x = 0$.

This procedure relies on two free parameters defining the training sets: the number of time snapshots $N_T$ used for supervision on $[t_0; t_i]$, and the ratio $\alpha = (t_i - t_0)/(t_f - t_0)$, which determines the fraction of the temporal interval allocated to the supervised phase relative to the extrapolation phase. An example configuration of these datasets is illustrated in figure 5.4a, with the choices $N_T = 30$ and $\alpha = 0.3$. The influence of these two parameters on the performance of the PINN will be analyzed in the following.

### ❖ The neural network architecture

The neural network architecture used in these experiments consists of three distinct parts (see figure 5.5): an input-normalization layer, an MLP network, and an output-transformation
<figure><img src="image_66.png" /><figcaption>Figure 5.5: Architecture of the neural network used to represent the hydrodynamic shock with the PINN strategy. It is composed of three parts: an input-normalization layer, an MLP neural network, and an output-transformation layer. tanh and lin correspond respectively to the hyperbolic-tangent activation function and the identity function.</figcaption></figure>

layer.

The purpose of input normalization is to map the temporal and spatial variables into the interval [0; 1], so as to avoid biases in the training of the MLP arising from scale differences. The transformations applied to the raw inputs (t, x) are:

$$x_1 = \frac{t - t_0}{t_f - t_0} \qquad (5.18)$$

$$x_2 = \frac{x - x_{\min}}{x_{\max} - x_{\min}} \qquad (5.19)$$

The core of the network, based on an MLP architecture, consists of an input layer with two neurons, each using the hyperbolic-tangent activation function. This layer is followed by a variable number of intermediate hidden layers, all with the same number of neurons and also using the hyperbolic-tangent activation. Finally, the output layer contains three neurons and uses the identity activation function.

The transformation applied to the outputs of the MLP serves two purposes: (i) to guarantee the positivity of physical quantities such as density $\rho$ and temperature T; (ii) to introduce adjustable multiplicative and additive factors for each output, enabling adaptation to the specific orders of magnitude of the predicted variables. The raw outputs of the MLP, denoted $(y_1, y_2, y_3)$, are therefore transformed as:

$$\rho = |w_{\rho}| \text{SP}(y_1) + |b_{\rho}|, \qquad (5.20)$$
<figure><img src="image_67.png" /><figcaption>Figure 5.6: Evaluation of the $L_{2,global}$ error for different numbers of hidden layers and different numbers of neurons per hidden layer.</figcaption></figure>

$$T = |w_T| \text{SP}(y_2) + |b_T|, \qquad (5.21)$$

$$v = |w_v| y_3 + |b_v|, \qquad (5.22)$$

where SP denotes the SoftPlus function, which guarantees the positivity of its output, and where the coefficients $w_\rho$, $w_T$, $w_v$, $b_\rho$, $b_T$, and $b_v$ are trainable parameters optimized during the training of the neural network.

In order to definitively fix the architecture of the MLP, we evaluated its ability to represent a solution exhibiting a sharp discontinuity, characteristic of the shock under consideration, by varying both the number of hidden layers and the number of neurons per hidden layer. To do so, we restricted the training to a purely supervised phase (parameter $\alpha = 1$, cost function $\mathcal{L} = \mathcal{L}_{IC}$), using only data taken from an independent test set composed of regularly spaced points over the space-time domain $(t, x) \in [t_0; t_f] \times [x_{\min}; x_{\max}]$. More precisely, we considered a regular Cartesian grid of 256 x 256 points. The density, velocity, and temperature profiles used as references come from the analytical solution of the problem, presented in Appendix D. Hydrodynamic shock on a wall. To assess the performance of the neural networks, we evaluated the following relative quadratic errors, computed over this test set:

$$L_{2,\rho} = \frac{\sum_{i=1}^{N_{test}} \{\rho_{pred}(t_i, x_i) - \rho_{sim}(t_i, x_i)\}^2}{\sum_{i=1}^{N_{test}} \rho_{sim}(t_i, x_i)^2}, \qquad (5.23)$$
$$L_{2,T} = \frac{\sum_{i=1}^{N_{test}} \{T_{pred}(t_i, x_i) - T_{sim}(t_i, x_i)\}^2}{\sum_{i=1}^{N_{test}} T_{sim}(t_i, x_i)^2}, \quad (5.24)$$

$$L_{2,v} = \frac{\sum_{i=1}^{N_{test}} \{v_{pred}(t_i, x_i) - v_{sim}(t_i, x_i)\}^2}{\sum_{i=1}^{N_{test}} v_{sim}(t_i, x_i)^2}, \quad (5.25)$$

where $X_{pred}$ denotes the quantities predicted by the neural network, $X_{sim}$ the reference quantities obtained from the analytical solution, and $N_{test} = 65 536$ the total number of test points. A global error, denoted $L_{2,global}$, can also be defined as the arithmetic mean of the three previous errors:

$$L_{2,global} = \frac{L_{2,\rho} + L_{2,T} + L_{2,v}}{3}. \quad (5.26)$$

The neural networks were trained successively using the Adam optimizer for 10 000 epochs (with a learning rate $\eta_{Adam} = 5 \times 10^{-3}$), followed by the L-BFGS optimizer for an additional 2 000 iterations (with a learning rate $\eta_{LBFGS} = 0.1$). The global errors $L_{2,global}$ obtained for the different architectures are shown in figure 5.6. The results show that all tested configurations are able to represent the discontinuity with relatively low error, except for certain configurations with only a single hidden layer. Consequently, and as a compromise between expressiveness and training cost, we retained an architecture with three hidden layers containing one neuron each for the remainder of the experiments.

## ❖ The impact of the different training parameters

Let us summarize here the three study parameters introduced in the previous sections, whose influence on the performance of the neural networks is analyzed: (i) the hyper-parameter $\lambda$, which appears in the adaptive weighting of the PDE residuals, (ii) the number of temporal sampling points $N_T$ used in the simulation data, and (iii) the ratio $\alpha$, which quantifies the relative duration of the supervised phase (training on simulation data) within the total temporal training window.

To study the impact of these parameters, we carried out a systematic exploration by training several neural networks for different combinations of values. Each training run followed the same protocol: an initial optimization phase using the Adam algorithm for 10 000 epochs (learning rate $\eta_{Adam} = 10^{-3}$), followed by a refinement stage using the L-BFGS algorithm for an additional 2 000 epochs (learning rate $\eta_{LBFGS} = 10^{-2}$). For each configuration tested, the retained model is the one that reached the lowest value of the global loss function $\mathcal{L}$.

In accordance with the methodology described above, we quantified the accuracy of the
<figure><img src="image_68.png" /><figcaption>Figure 5.7: Global $L_2$ error for different values of the hyper-parameter $\lambda$, while varying the number of temporal sampling points $N_T$ and the ratio $\alpha = (t_i - t_0)/(t_f - t_0)$. The simulation with $N_T = 1$ and $\alpha = 0$ corresponds to the classical PINN strategy, with the adaptive PDE weighting proposed by Liu (2023).</figcaption></figure>

predictions using the global error $L_{2,global}$, computed on the same test set. The results, shown in figure 5.7, indicate that increasing either $N_T$ or the ratio $\alpha$ tends to significantly improve the quality of the predictions. This highlights a key point: a small value of $\alpha$, corresponding to a short supervised phase, can be compensated for by a higher temporal density of simulation data. In other words, it is not necessary to extend the duration of the supervised phase to ensure reliable extrapolation, provided that it contains a sufficient number of temporal samples. Moreover, the analysis reveals the crucial role of the hyper-parameter $\lambda$: the best performances are consistently obtained for large values, typically $\lambda = 10^3$, suggesting that a stronger weighting of the differential residuals is essential to properly capture the extrapolated system dynamics. As illustrated in figure 5.8, the neural network accurately reproduces the shock structure when using $\lambda = 10^3$, $\alpha = 0.1$, and $N_T = 40$. The two curves are indistinguishable, demonstrating the quality of the extrapolation.
<figure><img src="image_69.png" /><figcaption>Figure 5.8: Comparison of the density, temperature, velocity, radiation temperature, and reduced flux profiles predicted by the neural network and obtained from the simulation, at the final time $t_f = 2 \times 10^{-4}$ s. The network was trained with the parameters $\alpha = 0.1$, $N_T = 40$, and $\lambda = 10^3$. The variables are shown in their rescaled form.</figcaption></figure>

Finally, it is worth noting that the so-called “classical” PINN configurations, namely $\alpha = 0$ and $N_T = 1$, systematically yield the poorest results, regardless of the value of $\lambda$. This behavior can be explained by the adaptive weighting strategy proposed by Liu (2023) [166]: although it effectively mitigates the impact of discontinuities during training, it tends, in turn, to weaken the representation of shock dynamics when the PDE-based constraints are not sufficiently informed by supervised data.

### 5.2.2 Radiative hydrodynamic tests

We then applied the same approach in order to extrapolate a radiative shock simulation, while keeping the same initial conditions as before. However, this part of the work remains
incomplete due to technical difficulties, that I describe below together with the methodological choices adopted.

The computational domain considered here extends spatially over $x \in [x_{min}; x_{max}] = [0; 30]$ m and temporally over $t \in [t_0; t_f] = [0; 2 \times 10^{-4}]$ s. The enlargement of the spatial domain was motivated by the need to correctly capture the entire radiative precursor. In this configuration, photon diffusion in the material is neglected, and we assume that their mean free path is independent of both radiation frequency and the hydrodynamic conditions of the fluid. This assumption makes it possible to assign the same value to the Planck and Rosseland mean opacities, namely $\kappa_R = \kappa_P = 1/\ell$, where $\ell$ denotes the photon mean free path. For the case studied, we set $\ell = 1$ m.

In what follows, I successively present the formulation of the radiative hydrodynamics equations employed, the details of the loss function used, the training data adopted, the neural network architecture, and finally the difficulties encountered with the solution currently obtained.

## ❖ The equations

First, let us rewrite the radiative hydrodynamics equations (2.99), adopting the M1-gray model for the description of radiation, in order to make explicit the time derivatives of the following physical variables: the density $\rho$, the fluid velocity $v$, the gas temperature $T$, the radiation temperature $T_R$, and the reduced flux $f_R$:

$$ \left\{ \begin{aligned} \partial_t \rho + v \partial_x \rho + \rho \partial_x v &= 0 \\ \partial_t v + v \partial_x v + \mathcal{R} \left( T \frac{\partial_x \rho}{\rho} + \partial_x T \right) &= S_v \\ \partial_t T + v \partial_x T + (\gamma - 1) T \partial_x v &= c S_T, \quad (5.27) \\ \partial_t T_R + c \left\{ f_R \partial_x T_R + \frac{T_R}{4} \partial_x f_R \right\} &= c S_{T_R} \\ \partial_t f_R + c \left\{ (\chi'_R - f_R) \partial_x f_R + \frac{4(\chi_R - f_R^2)}{T_R} \partial_x T_R \right\} &= c S_{f_R} \end{aligned} \right. $$

where $c$ is the speed of light, $\chi_R$ is the Eddington factor of the M1-gray model, and $\chi'_R$ its first derivative. These two quantities are given by:

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

$$ \chi'_R = \frac{2f_R}{\sqrt{4 - 3f_R^2}}, \quad (5.29) $$
Moreover, $S_v$, $S_T$, $S_{T_R}$ and $S_{f_R}$ denote the source terms that describe the interaction between radiation and matter. Under the assumptions made concerning the opacities, these source terms can be written, to order 1 in $v/c$, as:

$$S_v = \frac{a_R T_R^4}{\rho \ell} \left[ f_R - \frac{v}{c} \left\{ \chi_R + \left( \frac{T}{T_R} \right)^4 \right\} \right] , \quad (5.30)$$

$$S_T = \frac{(\gamma - 1) a_R T_R^4}{\rho \mathcal{R} \ell} \left[ 1 - \left( \frac{T}{T_R} \right)^4 - \frac{2\nu f_R}{c} \right] , \quad (5.31)$$

$$S_{T_R} = -\frac{T_R}{4\ell} \left[ 1 - \left(\frac{T}{T_R}\right)^4 - \frac{\nu f_R}{c} \right] , \quad (5.32)$$

$$S_{f_R} = -\frac{1}{\ell} \left[ f_R \left( \frac{T}{T_R} \right)^4 - \frac{v}{c} \left\{ \chi_R - f_R^2 + \left( \frac{T}{T_R} \right)^4 \right\} \right] , \quad (5.33)$$

where $a_R$ is the radiation constant. As before, the neural networks are trained using the rescaled version of these radiative hydrodynamics equations. The residuals of the differential equations then take the form:

$$\left\{
\begin{aligned}
\mathrm{res}_{\rho}(\tilde{t}, \tilde{x}) &= \partial_{\tilde{t}} \tilde{\rho} + \tilde{v} \partial_{\tilde{x}} \tilde{\rho} + \tilde{\rho} \partial_{\tilde{x}} \tilde{v} \\
\mathrm{res}_{v}(\tilde{t}, \tilde{x}) &= \partial_{\tilde{t}} \tilde{v} + \tilde{v} \partial_{\tilde{x}} \tilde{v} + \tilde{R} \left( \tilde{T} \frac{\partial_{\tilde{x}} \tilde{\rho}}{\tilde{\rho}} + \partial_{\tilde{x}} \tilde{T} \right) - \tilde{S}_v \\
\mathrm{res}_{T}(\tilde{t}, \tilde{x}) &= \tilde{c}^{-1} \partial_{\tilde{t}} \tilde{T} + \frac{\tilde{v}}{\tilde{c}} \partial_{\tilde{x}} \tilde{T} + \frac{(\gamma-1) \tilde{T}}{\tilde{c}} \partial_{\tilde{x}} \tilde{v} - \tilde{S}_T \\
\mathrm{res}_{T_R}(\tilde{t}, \tilde{x}) &= \tilde{c}^{-1} \partial_{\tilde{t}} \tilde{T}_R + f_R \partial_{\tilde{x}} \tilde{T}_R + \frac{\tilde{T}_R}{4} \partial_{\tilde{x}} f_R - \tilde{S}_{T_R} \\
\mathrm{res}_{f_R}(\tilde{t}, \tilde{x}) &= \tilde{c}^{-1} \partial_{\tilde{t}} f_R + \left( \partial_{f_R} \chi_R - f_R \right) \partial_{\tilde{x}} f_R + \frac{4(\chi_R - f_R^2)}{\tilde{T}_R} \partial_{\tilde{x}} \tilde{T}_R - \tilde{S}_{f_R}
\end{aligned}
\right. , \quad (5.34)$$

where the rescaled source terms are:

$$\tilde{S}_v = \frac{\tilde{a}_R \tilde{T}_R^4}{\tilde{\rho} \tilde{\ell}} \left[ f_R - \frac{\tilde{v}}{\tilde{c}} \left\{ \chi_R + \left( \frac{\tilde{T}}{\tilde{T}_R} \right)^4 \right\} \right] ,$$

$$\tilde{S}_T = \frac{(\gamma - 1) \tilde{a}_R \tilde{c} \tilde{T}_R^4}{\tilde{\rho} \tilde{R} \tilde{\ell}} \left[ 1 - \left( \frac{\tilde{T}}{\tilde{T}_R} \right)^4 - \frac{2\tilde{v}f_R}{\tilde{c}} \right] ,$$

$$\tilde{S}_{T_R} = \frac{\tilde{T}_R}{4\tilde{\ell}} \left[ 1 - \left( \frac{\tilde{T}}{\tilde{T}_R} \right)^4 - \frac{\tilde{v}f_R}{\tilde{c}} \right] ,$$

$$\tilde{S}_{f_R} = \frac{1}{\tilde{\ell}} \left[ f_R \left( \frac{\tilde{T}}{\tilde{T}_R} \right)^4 - \frac{\tilde{v}}{\tilde{c}} \left\{ \chi_R - f_R^2 + \left( \frac{\tilde{T}}{\tilde{T}_R} \right)^4 \right\} \right] .$$

The quantities with a tilde are rescaled variables defined as:

$$\tilde{t} = t/\dot{t} ,$$

$$\tilde{x} = x/(\dot{v}\tilde{t}) , \quad \tilde{\rho} = \rho/\dot{\rho} ,$$

$$\tilde{v} = v/\dot{v} , \quad \tilde{T} = T/\dot{T} , \quad \tilde{R} = \mathcal{R} \dot{T}/\dot{v}^2 ,$$
$$ \tilde{T}_R = T_R/\tilde{T} \ , \quad \tilde{c} = c/\tilde{v} \ , \quad \tilde{a}_R = a_R \tilde{T}^4/(\tilde{\rho} \tilde{v}^2) \ , $$

and $\tilde{i}$, $\tilde{\rho}$, $\tilde{v}$ and $\tilde{T}$ are scaling constants. In what follows, we adopt the values $\tilde{i} = 10^{-4}$, $\tilde{\rho} = 10$, $\tilde{v} = 10^4$ and $\tilde{T} = 10^5$. To simplify the notation, we will again omit the tildes hereafter; all quantities will be implicitly considered rescaled unless stated otherwise.

## ❖ The cost function

As in the hydrodynamic case, the loss function includes the PDE residual component $\mathcal{L}_{PDE}$, the prediction error relative to the simulation data $\mathcal{L}_{sim}$, and the boundary-condition error $\mathcal{L}_{BC}$. These three components can be written as:

$$ \begin{aligned} \mathcal{L}_{PDE} &= \frac{1}{N_{PDE}} \sum_{i=1}^{N_{PDE}} \left\{ \frac{\mathrm{res}_{\rho}(t_i, x_i)^2 + \mathrm{res}_{v}(t_i, x_i)^2 + \mathrm{res}_{T}(t_i, x_i)^2 + \mathrm{res}_{R}(t_i, x_i)^2 + \mathrm{res}_{f_R}(t_i, x_i)^2}{1 + \lambda(|\partial_{\tilde{x}}v_i| - \partial_{\tilde{x}}v_i)} \right\}, \\ \mathcal{L}_{sim} &= \frac{1}{N_{sim}} \sum_{i=1}^{N_{sim}} \left\{ \begin{aligned} & [\log(\rho_{pred}(t_i, x_i)) - \log(\rho_{sim}(t_i, x_i))]^2 + \\ & [\log(T_{pred}(t_i, x_i)) - \log(T_{sim}(t_i, x_i))]^2 + \\ & [\mathrm{asinh}(v_{pred}(t_i, x_i)) - \mathrm{asinh}(v_{sim}(t_i, x_i))]^2 + \\ & [\log(T_{R,pred}(t_i, x_i)) - \log(T_{R,sim}(t_i, x_i))]^2 + \\ & [\mathrm{asinh}(\alpha_{f_R} f_{R,pred}(t_i, x_i)) - \mathrm{asinh}(\alpha_{f_R} f_{R,sim}(t_i, x_i))]^2 \end{aligned} \right\}, \\ \mathcal{L}_{BC} &= \frac{1}{N_{BC}} \sum_{i=1}^{N_{BC}} \left\{ \mathrm{asinh}(v_{pred}(t_i, 0))^2 + \mathrm{asinh}(\alpha_{f_R} f_{R,pred}(t_i, 0))^2 \right\}, \end{aligned} $$

where the quantities $X_{pred}$ correspond to the values predicted by the neural network, and the quantities $X_{sim}$ denote the values taken from the simulation data. The parameter $\lambda$ is a training hyper-parameter to be tuned, whereas $\alpha_{f_R}$ is a factor introduced to increase the weight of the reduced flux in the loss function and, as we will see in the section devoted to the architecture, it also helps the neural network predict values with the correct order of magnitude. Indeed, the values of $f_R$ are very small, which tends to make their learning more difficult; introducing $\alpha_{f_R}$ therefore compensates for this imbalance. Empirically, we observed that choosing $\alpha_{f_R} = 10^5$ yields the best results.

The total loss function retains the same form as that used in the purely hydrodynamic test (see equation (5.17)), and we adopted the same weights $\omega_{sim} = 10$, $\omega_{PDE} = 1$, and $\omega_{BC} = 1$, for the same reasons as in the purely hydrodynamic case.
❖ The data

We performed this simulation using the HADES code, in a plane-parallel configuration. The spatial domain is discretized into 3 000 cells along the x-axis and 6 cells along the y-axis. One hundred output times were stored, uniformly distributed over the full duration of the simulation. The box dimensions are [$x_{min} = 0$; $x_{max} = 30$] m in the x direction and [$y_{min} = 0$; $y_{max} = 0.06$] m in the y direction, and the evolution is computed from $t_0 = 0$ s up to $t_f = 2 \times 10^{-4}$ s.

As in the purely hydrodynamic case, three distinct datasets are used to train the neural network (see figure 5.4b):

1. HADES simulation data (blue points).

These are used to evaluate the cost term $\mathcal{L}_{sim}$. This dataset covers the entire spatial domain and all times within [$t_0$; $t_i$]. The network is trained on $N_T$ timesteps selected in this interval, and all spatial grid points are retained ($N_x = 3000$).

2. Points for the radiation hydrodynamics equations (green points).

They are used to evaluate the PDE residuals in the cost term $\mathcal{L}_{PDE}$. These points are distributed over the entire spatial domain and over all times in [$t_0$; $t_f$]. We sampled 10 000 points randomly in this domain using a quasi-random Hammersley sequence.

3. Points for boundary conditions (orange points).

These are used in the cost term $\mathcal{L}_{BC}$. The points are located at $x = 0$ m and span all times in [$t_0$; $t_f$]. We use 100 time points in this dataset.

In this case, we chose to extend the last two datasets to the entire temporal domain ($t \in [t_0; t_f]$), rather than restricting them only to the extrapolation interval ($t \in [t_i; t_f]$), as was done in the purely hydrodynamic test. This modification proved essential: it allows the neural network to correctly reproduce the simulation data and to satisfy the physical equations for times $t \in [t_0; t_i]$. Without this constraint, the extrapolation obtained beyond $t_i$ becomes completely unphysical. Two parameters therefore need to be explored, as in the purely hydrodynamic case:

• $\alpha = (t_i - t_0)/(t_f - t_0)$, which sets the fraction of the time window devoted to supervised training;

• $N_T$, the number of simulation timesteps used in this supervised phase.
<figure><img src="image_70.png" /><figcaption>Figure 5.9: Neural-network architecture used to represent the radiative shock with the PINN strategy. It is composed of two MLP neural networks: the first predicts the density and the velocity, while the second predicts the gas temperature, the radiation temperature, and the reduced flux. Each of these networks consists of three parts: an input-normalization layer, an MLP core, and an output-transformation layer. The symbols tanh and lin denote the hyperbolic-tangent activation function and the identity function, respectively.</figcaption></figure>

❖ The neural network architecture

For these experiments, we used an architecture composed of two independent MLP neural networks (see figure 5.9):

• the first predicts the fluid density $\rho$ and velocity $v$,

• the second predicts the gas temperature $T$, the radiation temperature $T_R$, and the reduced flux $f_R$.

We opted for this separation of the neural networks because of the similarity observed
between the profiles of $\rho$ and $v$ on the one hand, and of $T$ and $T_R$ on the other. Moreover, we chose to assign the prediction of the reduced flux $f_R$ to the network responsible for predicting $T$ and $T_R$. Indeed, our tests showed that, for a sufficiently expressive architecture, it is possible to obtain an accurate prediction of this quantity despite its different profile.

Each network consists of three modules: an input-normalization layer, the MLP core, and finally an output-transformation layer. The input normalization layer rescales the spatio-temporal variables into a common interval [0; 10], which facilitates training by limiting scale effects. This strategy proved more effective than the one used in the purely hydrodynamic case, and is written:

$$x_1 = 10 \frac{t - t_0}{t_f - t_0} \qquad (5.35)$$

$$x_2 = 10 \frac{x - x_{\min}}{x_{\max} - x_{\min}} \qquad (5.36)$$

The output transformation layer of the MLPs serves the same purposes as in the purely hydrodynamic case: ensuring the positivity of physical quantities and adapting their scales. The transformations are defined as:

$$\rho = |w_{\rho}| \text{SP}(y_1) + |b_{\rho}|, \qquad (5.37)$$

$$v = |w_v| y_2 + |b_v|, \qquad (5.38)$$

$$T = |w_T| \text{SP}(y_3) + |b_T|, \qquad (5.39)$$

$$T_R = |w_T| \text{SP}(y_4) + |b_T|, \qquad (5.40)$$

$$f_R = \frac{|w_{f_R}| y_5}{\alpha_{f_R}}. \qquad (5.41)$$

Here, SP denotes the SoftPlus activation, chosen to guarantee positivity of the relevant output variables. The coefficients $w_\rho$, $w_T$, $w_v$, $w_{f_R}$ and their associated biases $b_\rho$, $b_T$, $b_v$ are free parameters optimized during training. In particular, the scaling factor $\alpha_{f_R}$, already introduced in the loss term $\mathcal{L}_{sim}$, is intended to help the network predict the very small values often taken by the reduced flux $f_R$.

To evaluate prediction quality, we introduce the following relative squared errors, computed on a test set containing all HADES simulation data:

$$L_{2,T_R} = \frac{\sum_{i=1}^{N_{test}} \{T_{R,pred}(t_i, x_i) - T_{R,sim}(t_i, x_i)\}^2}{\sum_{i=1}^{N_{test}} T_{R,sim}(t_i, x_i)^2}, \qquad (5.42)$$

$$L_{2,f_R} = \frac{\sum_{i=1}^{N_{test}} \{f_{R,pred}(t_i, x_i) - f_{R,sim}(t_i, x_i)\}^2}{\sum_{i=1}^{N_{test}} f_{R,sim}(t_i, x_i)^2}, \qquad (5.43)$$
<figure><img src="image_71.png" /><figcaption>(a) Evaluation of the $L_{2,MLP_1}$ error for different numbers of hidden layers and different numbers of neurons per hidden layer.</figcaption></figure>

<caption>(b) Evaluation of the $L_{2,MLP_2}$ error for various numbers of hidden layers and neurons per hidden layer.</caption>

Figure 5.10: Evaluation of the $L_2$ error of the MLP$_1$ and MLP$_2$ neural networks as a function of the number of hidden layers and the number of neurons per hidden layer.

where the quantities indexed “pred” come from the neural network predictions, and those indexed “sim” come from the test set. The total number of samples is $N_{test} = 303 000$. From these individual errors, we define aggregated errors:

$$L_{2,MLP_1} = \frac{L_{2,\rho} + L_{2,v}}{2}, \qquad (5.44)$$

$$L_{2,MLP_2} = \frac{L_{2,T} + L_{2,T_R} + L_{2,f_R}}{3}, \qquad (5.45)$$

$$L_{2,global} = \frac{2L_{2,MLP_1} + 3L_{2,MLP_2}}{5}, \qquad (5.46)$$

where $L_{2,\rho}$, $L_{2,v}$, and $L_{2,T}$ were defined in equations (5.23), (5.25), and (5.24). $L_{2,MLP_1}$ evaluates the performance of MLP$_1$, $L_{2,MLP_2}$ that of MLP$_2$, and $L_{2,global}$ measures the overall performance of the strategy.

In this architecture, two hyper-parameters remain to be chosen for each network: the number of hidden layers and the number of neurons per hidden layer. As in the previous study, we explored several configurations, training each architecture on the full simulation dataset (setting $\alpha = 1$ and $N_T = 101$). Training proceeds in two stages:

* a first phase using the Adam optimizer for 10 000 epochs (learning rate $\eta_{Adam} = 5 \times 10^{-3}$),

* followed by a second phase using the L-BFGS optimizer for 2 000 additional iterations (learning rate $\eta_{LBFGS} = 10^{-2}$).
<figure><img src="image_72.png" /><figcaption>(a) Density</figcaption></figure>

<figure><img src="image_73.png" /><figcaption>(c) Velocity</figcaption></figure>

<figure><img src="image_74.png" /><figcaption>(b) Gas temperature</figcaption></figure>

<figure><img src="image_75.png" /><figcaption>(d) Radiative temperature</figcaption></figure>

<figure><img src="image_76.png" /><figcaption>(e) Reduced flux</figcaption></figure>

Figure 5.11: Comparison of the density, temperature, velocity, radiative temperature, and reduced flux profiles predicted by the neural network and obtained from the simulation, at the final time $t_f = 2 \times 10^{-4}$ s. The network was trained exclusively on simulation data.

Once the training phases were completed, we evaluated the prediction errors in order to analyze the performance of the different neural-network architectures. To do so, we examined separately the errors $L_{2,MLP_1}$ and $L_{2,MLP_2}$, associated respectively with each network. For MLP1, the results shown in figure 5.10a indicate that the error $L_{2,MLP_1}$ decreases as both the number of hidden layers and the number of neurons per layer increase. Based on these observations, we selected an architecture composed of two hidden layers, each containing 20 neurons, since this configuration yielded a low error on the test set. Similar conclusions were reached for the MLP2 network, as illustrated by the results in figure 5.10b. We therefore also adopted, for this network, an architecture with two hidden layers of 20 neurons each.
Figure 5.11 shows that the neural networks trained from the simulation data are, overall, able to correctly reproduce the profiles of the physical quantities considered. The predictions are satisfactory for the fluid velocity, the temperatures T and $T_R$, as well as for the reduced flux $f_R$. However, the density is less accurately captured, particularly in the shocked region of the fluid.

### ❖ Training difficulties observed with the PINN strategy

We performed a test with the parameters $\lambda = 10^3$, $\alpha = 0.5$ and $N_T = 50$, representative of the best results I obtained so far. The analysis of the PINN neural-network predictions (see figure 5.12) shows that it faithfully reproduces the dynamics over the time interval $t \in [0; 1 \times 10^{-4}]$ s (in physical units), i.e. the region on which it was trained using data from the HADES code.

However, beyond this interval, at times such that $t \in [1 \times 10^{-4}; 2 \times 10^{-4}]$ s, although the network is still able to correctly capture the propagation of the shock, two limitations clearly appear: (i) it does not correctly reproduce the dynamics of the radiative precursor front, which it extrapolates almost linearly whereas it should slow down; (ii) it provides inaccurate predictions of the hydrodynamic and radiative values downstream of the shock.

These difficulties can be explained by two major characteristics of the physical situation studied: (i) the presence of a shock, which introduces a discontinuity in the hydrodynamic fields; (ii) the very stiff nature of the equations, particularly at the edge of the radiative precursor, where the source terms reach high values, typically $S_v \approx 10^3$ m/s$^2$, $cS_T \approx 10^{11}$ K/s, $cS_{T_R} \approx 10^{12}$ K/s and $cS_{f_R} \approx 10^5$ s$^{-1}$, to be compared with the values of the physical quantities in this region, which are $v \approx 10^4$ m/s, $T \approx T_R \approx 10^4$ K, and $f_R \approx 10^{-3} \in [-1; 1]$.

In the experiment presented here, the evolution equations for the temperature T, the radiative temperature $T_R$, and the reduced flux $f_R$ were divided by the speed of light c (see equation (5.34)). This normalization artificially reduces the stiffness and makes the network training considerably easier. However, we believe that it decreases the relative importance of the corresponding equations in the loss function and may be one of the causes of the poor representation of the radiative precursor dynamics, since these three quantities are precisely those most strongly influenced by the precursor. Without this weighting by c, training becomes extremely difficult, to the point that the neural network is no longer even able to reproduce even the simulation data.

Moreover, in this work, we applied the adaptive-residual weighting method proposed by
<figure><img src="image_77.png" /><figcaption>Figure 5.12: Comparison between the prediction of the PINN neural network and the data from the HADES simulation. The neural network was trained on the HADES simulation data for times $t \in [0; 1 \times 10^{-4}]$ s.</figcaption></figure>
Liu (2023) [166], which improves the representation of discontinuities in the hydrodynamic quantities. Nevertheless, it remains necessary to develop more suitable strategies capable of simultaneously handling equation stiffness and discontinuities, in order to better capture the structure of the radiative precursor. Designing such an approach represents an essential direction for future work in this field.

Finally, in future developments, it will be crucial to ensure that the reduced flux $f_R$ and the ratio $T/T_R$ are predicted with very high accuracy by the network, since these quantities directly enter the source terms $S_v$, $cS_T$, $cS_{TR}$ and $cS_{fR}$. As can be seen from their expressions (see equations (5.30)–(5.33)), even a small error in the estimation of $f_R$ or $T/T_R$ can lead to significant deviations in the estimation of these source terms.

## 5.3 Synthesis

In order to overcome the limitation imposed by the CFL condition of the HADES code, I explored the use of PINN neural networks with the objective of extrapolating, to later times, data produced by pre-existing numerical simulations, in particular those obtained with HADES.

In the case of pure hydrodynamics, this approach proved satisfactory: it makes it possible to efficiently extrapolate the simulated profiles and to faithfully recover the structure of the hydrodynamic shock. In contrast, the results obtained for radiative shocks remain, at this stage, less satisfactory. This difficulty appears to be linked to the combination of two factors: on the one hand, the presence of a pronounced discontinuity at the shock front, and on the other hand, the strong stiffness of the source terms describing the interactions between radiation and fluid within the radiative precursor. These terms are extremely sensitive to the ratio between the gas temperature and the radiation temperature, $T/T_R$, as well as to the reduced flux $f_R$, which complicates the training of the neural network.

To improve the quality of extrapolation in the context of radiative shocks, it therefore appears necessary to develop a method capable of simultaneously handling both the shock discontinuity and the stiffness of these source terms, in order to ensure a coherent and stable representation of the fluid–radiation coupling.
<figure><img src="image_78.png" /></figure>
# CONCLUSIONS AND PERSPECTIVES

## Conclusions

The simulation code HADES is a powerful and accurate tool for studying the interactions between a fluid and radiation, particularly thanks to its ability to accurately model radiative transfer using two different models: the M1-gray model and the M1-multigroup model. However, this accuracy comes at a high computational cost, especially in the case of the M1-multigroup model. One of the main bottlenecks lies in the computation of the Eddington factor in the M1-multigroup model, an essential ingredient of the model's closure relation but one for which no analytical expression exists. Until now, three approaches were possible: search algorithms, which are very accurate but extremely costly in computational time; interpolation on precomputed databases, which is fast and accurate but requires prior knowledge of the radiative quantities and their orders of magnitude, a condition rarely met, especially in a multigroup context; and finally, the use of the analytical Eddington factor from the M1-gray model in multigroup simulations, a fast but very inaccurate approach.

After a detailed study of the dependence of the Eddington factor on the radiative quantities and on the frequency bounds of a considered group, I developed a method based on the use of artificial intelligence, in which neural networks were trained using accurate data produced by search algorithms. This approach turns out to be about 3 000 times faster than these same search algorithms, while being significantly more accurate than using the factor from the M1-gray model. Moreover, it is fully general, since it requires no prior knowledge of the
characteristics of the radiation, unlike interpolation techniques.

This improvement of the HADES code enabled me to conduct a detailed study of the influence of the spectral character of radiation on the dynamics and structure of radiative shocks. I highlighted three major effects induced by fine spectral modelling with the M1-multigroup model, compared with the M1-gray model: a slower shock, a lower temperature and a higher density in the shocked region, as well as an enlargement of the radiative precursor, especially of its non-equilibrium zone. Furthermore, the comparison between the jump relations obtained from the M1-multigroup simulations and those predicted by a model based on the diffusion approximation reveals non-negligible discrepancies, on the order of 3–15% for the density and 1–8% for the temperature. These results emphasise the need for caution when using simplified analytical models based on the diffusion approximation.

Beyond numerical considerations, these results have major astrophysical implications, since radiative shocks are ubiquitous throughout the stellar life cycle. During star formation, for instance, the accretion shocks associated with the appearance of the first hydrostatic core play a crucial role in determining the thermal structure of protostars. Our simulations indicate that when spectral effects are treated more accurately, the energy dissipation is greater than predicted by the M1-gray model, implying cooler and more compact protostars.

In supernova remnants, hydrodynamic instabilities such as the Vishniac instability [41] or the detonation-wave instability observed experimentally by Grun [177] are strongly influenced by radiative losses, which modify the dynamics of the shocked layers. A faithful modelling of radiative transport, made possible by the M1-multigroup formalism, could help reproduce the filamentary structures that are observed but not yet obtained in current simulations.

Finally, in binary systems containing a magnetised accreting white dwarf, simulations based on parametrised cooling functions predict quasi-periodic oscillations in the light curve, resulting from a radiative-shock cooling instability [75]. However, these models struggle to reproduce observations across the full electromagnetic spectrum [178], which may stem from an overly simplified treatment of radiative processes. The M1-multigroup model developed in this work represents a methodological advance that could improve the realism of such simulations. Nevertheless, although such results are now accessible with the HADES code, a systematic study of radiative shocks across a wide variety of configurations, particularly for varying Mach numbers, different opacity laws, or in 2D geometry, remains difficult, if not impossible, due to the still very high computational cost. This limitation arises mainly from the Courant–Friedrichs–Lewy condition imposed by the explicit schemes used in HADES, which requires extremely small time steps to ensure the stability of the radiative transfer resolution.
To overcome this obstacle, I explored the potential of using PINN neural networks, which differ from classical architectures by explicitly embedding the physical equations within their loss function. The goal is to enable temporal extrapolation of simulations at a cost far lower than that required to continue the computation directly with HADES. To assess the potential of this approach, I first tested it on the case of a strong hydrodynamic shock. The results showed that it yielded a faithful extrapolation of the simulation data, including at times significantly beyond those initially provided.

However, its application to radiative shocks proved more challenging: the extrapolation performance is, at this stage, noticeably less satisfactory. This difficulty arises from the intrinsic complexity of the problem, which combines the presence of a discontinuity with the pronounced stiffness of the differential equations. This situation highlights the need for further methodological development to adapt the PINN approach effectively to the specificities of radiative hydrodynamics.

## Perspectives

The application of artificial intelligence to radiative hydrodynamics simulation is a field rich with promising perspectives. Much work remains to be done to fully exploit these approaches, both from a methodological standpoint and to improve the physical understanding of the modelled phenomena. As an illustration, several research directions deserve to be explored in the coming years:

* Improving PINN-based strategies to allow for more reliable extrapolation of radiative hydrodynamics simulations, particularly in complex contexts such as radiative shocks. A complementary line of research would be to develop uncertainty or prediction-error estimators, using probabilistic architectures such as Bayesian Neural Networks, for instance, in order to quantify the confidence in extrapolated predictions;

* Exploring the use of neural operators (Neural Operators), which enable learning solution mappings of parametrised problems in operator form, for example by directly predicting the next time step from the current state as a replacement for classical numerical schemes. Such an approach can significantly reduce the computational cost of simulation codes and is already being used in other simulation frameworks, as suggested in recent work by Azizzadenesheli et al. (2024) [179];
* Investigating the possibility of using machine learning methods to improve the computation of the mean opacities $κ_E$ and $κ_F^{(i)}$ appearing in the M1-gray and M1-multigroup models. The goal would be to avoid the common approximation of replacing them with Planck and Rosseland means, and instead compute more realistic weighted opacities from opacity tables (e.g. TOPS), depending on the local shape of the specific intensity in the M1 model;

* Considering the use of neural networks trained on experimental data or on data from reference codes that solve the specific intensity (for example via characteristics or Monte Carlo methods) in order to provide a more accurate and realistic estimate of the closure relation (Eddington tensor). Such an approach could potentially overcome some intrinsic limitations of the M1 model, such as the poor representation of the interaction of two light beam, by enriching the formalism with learned closure relations better suited to complex regimes in 2D or 3D configurations.
BASIC CONCEPTS ON THE PLANCK FUNCTION

In this appendix, we detail several calculation elements related to the Planck function, which will be particularly useful in Appendix B. *Numerical calculation of the radiative quantities*, devoted to the numerical aspects of radiative transfer. Before delving into these calculations, let us recall that the Planck function describes the spectral distribution of the energy intensity emitted by a blackbody in thermal equilibrium at a given temperature. In dimensionless form, it can be written as follows:

$$b(x) = \frac{15}{\pi^4} \frac{x^3}{e^x - 1}, \qquad (A.1)$$

We will first present the analytical expressions of the successive derivatives of the function $b$, then detail the computation of its antiderivative, as well as the way to approximate certain of its integrated forms, which are required for the numerical treatments of radiative transfer.

## A.1 Derivatives of the Planck function

Let us begin by presenting the first three derivatives of the Planck function, which appear in particular in asymptotic expansions and certain numerical approximations. The three first derivatives of $b$ are written as follows:

$$b'(x) = \frac{15}{\pi^4} x^2 \left( \frac{3-x}{(e^x-1)} - \frac{x}{(e^x-1)^2} \right),$$

$$b''(x) = \frac{15}{\pi^4} x \left( \frac{2x^2}{(e^x-1)^3} + 3x \frac{x-2}{(e^x-1)^2} + \frac{x^2-6x+6}{e^x-1} \right),$$
$$b^{(3)}(x) = -\frac{15}{\pi^4} \left( \frac{6x^3}{(e^x - 1)^4} + 6x^2 \frac{2x - 3}{(e^x - 1)^3} + x \frac{7x^2 - 27x + 18}{(e^x - 1)^2} + \frac{x^3 - 9x^2 + 18x - 6}{e^x - 1} \right).$$

## A.2 Antiderivatives of the Planck function

Let us now introduce the antiderivative $\mathcal{B}$ of the Planck function $b$, which is expressed in the following manner:

$$\mathcal{B}(x) = \int_{0}^{x} b(z) \, dz \quad . \qquad (A.2)$$

This antiderivative admits an exact expression involving polylogarithm functions:

$$\mathcal{B}(x) = \frac{15e^{-x}}{\pi} [-6\mathrm{Li}_4(e^{-x}) - 6x\mathrm{Li}_3(e^{-x}) - 3x^2\mathrm{Li}_2(e^{-x}) - x^3\mathrm{Li}_1(e^{-x})] + C^P, \quad (A.3)$$

where $C^P$ is an integration constant. By convention, it is set to $C^P = 1$, so that $\mathcal{B}(+\infty) = 1$. The polylogarithm function of order $n$, denoted $\mathrm{Li}_n$, is defined for $|x| < 1$ by the series:

$$\mathrm{Li}_n(x) = \sum_{k=1}^{\infty} \frac{x^k}{k^n} \quad . \qquad (A.4)$$

However, this exact expression cannot be used directly in numerical contexts due to the infinite nature of the series. It is therefore necessary to truncate it to a finite number of terms, denoted $N_P$, which leads to using the following approximation:

$$\mathrm{Li}_n^{N_P}(x) = \sum_{k=1}^{N_P} \frac{x^k}{k^n} \quad . \qquad (A.5)$$

This approximation is, however, problematic in certain cases, particularly for derived quantities such as the integral of $z \mapsto \mathcal{B}(z)/z^4$, used in the computation of the radiative energy. Indeed, the truncation may introduce a non-physical divergence, as can be seen in figure A.1. To remedy this difficulty, Hung Chinh Nguyen (2011) [37] proposes an alternative method based on the Taylor expansion of the Planck function close to $x = 0$. This approach makes it possible to obtain an analytical approximation of $\mathcal{B}$ by term-by-term integration of the Taylor expansion of $b$:

$$b_T^{N_T}(x) = \frac{15}{\pi^4} \sum_{k=2}^{N_T} T_k x^k, \qquad (A.6)$$

where $N_T$ is the order of the retained series, and the coefficients $T_k$ are those of the Taylor expansion of the function $b$ and are such that:

$$T_2 = 1, \qquad T_3 = -\frac{1}{2}, \qquad T_4 = \frac{1}{12}, \qquad T_5 = 0,$$
<figure><img src="image_79.png" /><figcaption>Figure A.1: Logarithm of the relative error on the estimation of the functions $\mathcal{B}$, $\mathcal{H}_2$, $\mathcal{H}_3$ and $\mathcal{H}_4$, obtained using two methods: the approximation by polylogarithms truncated at $N_P$ terms (solid lines), and the approximation by a Taylor expansion of order $N_T$ for $\mathcal{B}$ (dashed lines). The error is computed with respect to the truncated polylogarithm approximation, using $N_P = 100 000$.</figcaption></figure>

$$ 
\begin{array}{l@{\quad}l@{\quad}l@{\quad}l}
T_6 = -\frac{1}{720} & , & T_7 = 0 & , & T_8 = \frac{1}{30240} & , & T_9 = 0 & , \\
T_{10} = -\frac{1}{1209600} & , & T_{11} = 0 & , & T_{12} = \frac{1}{47900160} & . & & 
\end{array}
$$

The integration of this Taylor expansion then provides an approximation of $\mathcal{B}$, valid for small values of $x$:

$$ \mathcal{B}_T^{N_T}(x) = \frac{15}{\pi^4} \sum_{k=3}^{N_T+1} \frac{T_{k-1}x^k}{k} + C^T . \qquad (A.7) $$

where $C^T$ is an integration constant that we set to zero in order to ensure that $\mathcal{B}_T^{N_T}(0) = 0$. By combining the two approaches, Taylor expansion of order $N_T$ for $x < x^*$ and polylogarithm functions truncated at $N_P$ for $x > x^*$, it is possible to construct an accurate global approximation
<figure><img src="image_80.png" /><figcaption>Figure A.2: Explanatory diagram of the computation of the integrals $H_i$, using Taylor expansions, truncated polylogarithm functions, or a combination of both.</figcaption></figure>

of the function $\mathcal{B}$. For example, by taking $N_P = 20$, $N_T = 12$, and a transition point $x^* = 1$, the maximum approximation error remains below $10^{-11}$, as illustrated in figure A.1a.

## A.3 Functions derived from the Planck function

Let us introduce a family of antiderivatives associated to the Planck function $\mathcal{B}$, defined by:

$$ H_i(x) = \int \frac{\mathcal{B}(z)}{z^i} dz . \qquad (A.8) $$

In this section, I revisit and expand upon the results presented by Hung Chinh Nguyen (2011) [37], which I have further developed in order to enable a more efficient computation of these integrals in the case where the integration bounds $x_1$ and $x_2$ are very close. I will also present the expressions of the successive derivatives of these integrals.
## ❖ Numerical computation of the functions $H_4$, $H_3$ and $H_2$

Using the approximations of the Planck function $\mathcal{B}$ through truncated polylogarithm functions (equation (A.5)) or through a Taylor series expansion (equation (A.6)), the functions $H_i$ can be estimated using the following expressions:

$$ H_i(x) = \begin{cases} \mathcal{H}_{i,P}^{N_P}(x) & \text{if } x \le x^* \\ \mathcal{H}_{i,T}^{N_T}(x) & \text{otherwise} \end{cases}. \qquad (A.9) $$

where the subscripts $P$ and $T$ respectively denote the truncated polylogarithm functions and the Taylor expansion, and $x^*$ is the matching point between their two domains of validity. This computation is illustrated in figure A.2. For the specific cases of the functions $H_2$, $H_3$, and $H_4$, one obtains the following explicit expressions for the truncated polylogarithm functions:

$$ \mathcal{H}_{2,P}^{N_P}(x) = \frac{15}{\pi^4} \left[ 6x^{-1} \mathrm{Li}_4^{N_P}(e^{-x}) + 4 \mathrm{Li}_3^{N_P}(e^{-x}) + x \mathrm{Li}_2^{N_P}(e^{-x}) \right] - \frac{1}{x} + C_2^P, $$

$$ \mathcal{H}_{3,P}^{N_P}(x) = \frac{15}{\pi^4} \left[ 3x^{-2} \mathrm{Li}_4^{N_P}(e^{-x}) + 3x^{-1} \mathrm{Li}_3^{N_P}(e^{-x}) + \mathrm{Li}_2^{N_P}(e^{-x}) \right] - \frac{1}{2x^2} + C_3^P, $$

$$ \mathcal{H}_{4,P}^{N_P}(x) = \frac{15}{\pi^4} \left[ 2x^{-3} \mathrm{Li}_4^{N_P}(e^{-x}) + 2x^{-2} \mathrm{Li}_3^{N_P}(e^{-x}) + x^{-1} \mathrm{Li}_2^{N_P}(e^{-x}) \right] - \frac{1}{3x^3} + C_4^P, $$

and, for the Taylor expansions:

$$ \mathcal{H}_{2,T}^{N_T}(x) = \frac{15}{\pi^4} \sum_{k=2}^{N_T} \frac{T_k x^k}{k(k+1)} + C_2^T, $$

$$ \mathcal{H}_{3,T}^{N_T}(x) = \frac{15}{\pi^4} \sum_{k=1}^{N_T-1} \frac{T_{k+1} x^k}{k(k+2)} + C_3^T, $$

$$ \mathcal{H}_{4,T}^{N_T}(x) = \frac{15}{\pi^4} \left[ \frac{T_2 \ln(x)}{3} + \sum_{k=1}^{N_T-2} \frac{T_{k+2} x^k}{k(k+3)} \right] + C_4^T. $$

The terms $C_i^P$ and $C_i^T$ are integration constants. As for the function $\mathcal{B}$, one can verify that choosing $N_P = 20$, $N_T = 12$, and $x^* = 1$ ensures excellent accuracy (maximum error below $10^{-12}$; see figure A.1). For the computation of an integral of the form $\int_{x_1}^{x_2} \frac{\mathcal{B}(z)}{z^i} dz$ with $x_2 \ge x_1$, we introduce:

$$ x_0 = \frac{x_2 + x_1}{2}, $$

$$ x_T = \frac{x^* + x_1}{2}, $$

$$ x_P = \frac{x_2 + x^*}{2}, $$
$$ \delta = \frac{x_2 - x_1}{2x_0} , $$

$$ \delta_T = \frac{x^* - x_1}{2x_T} , $$

$$ \delta_P = \frac{x_2 - x^*}{2x_P} , $$

and we write:

$$ \mathcal{H}_i(x_1; x_2) = \int_{x_1}^{x_2} \frac{\mathcal{B}(z)}{z^i} dz = \delta e^{-x_1} x_0^{1-i} \mathcal{H}_i^*(x_1; x_2) , $$

where the function $\mathcal{H}_i^*$ depends on the domain in which the bounds $x_1$ and $x_2$ lie:

$$ \mathcal{H}_i^*(x_1; x_2) = \begin{cases} \mathcal{H}_{i,P}^*(x_1; x_2) & \text{if } x_1 > x^* \text{ and } x_2 > x^* \\ \mathcal{H}_{i,T}^*(x_1; x_2) & \text{if } x_1 \le x^* \text{ and } x_2 \le x^* \\ \frac{\delta_T}{\delta} \left(\frac{x_0}{x_T}\right)^{i-1} \mathcal{H}_{i,T}^*(x_1; x^*) + e^{x_1 - x^*} \frac{\delta_T}{\delta} \left(\frac{x_0}{x_P}\right)^{i-1} \mathcal{H}_{i,P}^*(x^*; x_2) & \text{otherwise} \end{cases} . $$

Letting $x_l$ and $x_h$ denote the lower and upper bounds (equal respectively to $x_1$ or $x^*$, and to $x_2$ or $x^*$), and letting $x_m$ and $\delta_m$ be variables such that $x_m = (x_h + x_l)/2$ and $\delta_m = (x_h - x_l)/(2x_m)$, the functions $\mathcal{H}_{i,P}^*$ and $\mathcal{H}_{i,T}^*$ are given by:

$$ \mathcal{H}_{i,P}^*(x_l; x_h) = e^{x_l} \delta_m^{-1} x_m^{i-1} \left( \mathcal{H}_{i,P}^{N_P}(x_h) - \mathcal{H}_{i,P}^{N_P}(x_l) \right) , $$

$$ \mathcal{H}_{i,T}^*(x_l; x_h) = e^{x_l} \delta_m^{-1} x_m^{i-1} \left( \mathcal{H}_{i,T}^{N_T}(x_h) - \mathcal{H}_{i,T}^{N_T}(x_l) \right) . $$

For the functions $\mathcal{H}_2$, $\mathcal{H}_3$, and $\mathcal{H}_4$, the differences $\mathcal{H}_{i,P}^*$ can be efficiently evaluated by expanding the truncated polylogarithm functions:

$$ \mathcal{H}_{2,P}^*(x_l; x_h) = - \frac{15}{\pi^4} \sum_{k=1}^{N_P} e^{-(k-1)x_l} \left[ 6 \frac{h_k^1(x_m, \delta_m)}{k^4} + 4x_m \frac{\Delta_k(x_m, \delta_m)}{k^3} + x_m^2 \frac{h_k^2(x_m, \delta_m)}{k^2} \right] + \frac{2e^{x_l}}{1 - \delta_m^2} , $$

$$ \mathcal{H}_{3,P}^*(x_l; x_h) = - \frac{15}{\pi^4} \sum_{k=1}^{N_P} e^{-(k-1)x_l} \left[ 3 \frac{h_k^3(x_m, \delta_m)}{k^4} + 3x_m \frac{h_k^1(x_m, \delta_m)}{k^3} + x_m^2 \frac{\Delta_k(x_m, \delta_m)}{k^2} \right] + \frac{2e^{x_l}}{(1 - \delta_m^2)^2} , $$

$$ \mathcal{H}_{4,P}^*(x_l; x_h) = - \frac{15}{\pi^4} \sum_{k=1}^{N_P} e^{-(k-1)x_l} \left[ 2 \frac{h_k^4(x_m, \delta_m)}{k^4} + 2x_m \frac{h_k^3(x_m, \delta_m)}{k^3} + x_m^2 \frac{h_k^1(x_m, \delta_m)}{k^2} \right] + \frac{2e^{x_l}(3 + \delta_m^2)}{3(1 - \delta_m^2)^3} , $$

where the auxiliary functions $h_k^j$ are defined by:

$$ h_k^1(x_m, \delta_m) = \frac{\Delta_k(x_m, \delta_m) + \Sigma_k(x_m, \delta_m)}{1 - \delta_m^2} , $$

$$ h_k^2(x_m, \delta_m) = \Delta_k(x_m, \delta_m) - \Sigma_k(x_m, \delta_m) , $$

$$ h_k^3(x_m, \delta_m) = \frac{(1 + \delta_m^2)\Delta_k(x_m, \delta_m) + 2\Sigma_k(x_m, \delta_m)}{(1 - \delta_m^2)^2} , $$
$$h_k^4(x_m, \delta_m) = \frac{(1 + 3\delta_m^2)\Delta_k(x_m, \delta_m) + (3 + \delta_m^2)\Sigma_k(x_m, \delta_m)}{(1 - \delta_m^2)^3}.$$

The functions $\Delta_k$ and $\Sigma_k$ are given by:

$$\Delta_k(x_m, \delta_m) = \frac{1 - e^{-2kx_m\delta_m}}{\delta_m}, \qquad (A.10)$$

$$\Sigma_k(x_m, \delta_m) = 1 + e^{-2kx_m\delta_m}, \qquad (A.11)$$

For $kx_m\delta_m \ll 1$, one may use the approximation:

$$\Delta_k(x_m, \delta_m) \approx 2kx_m \left( 1 - kx_m\delta_m + \frac{2(kx_m\delta_m)^2}{3} \right).$$

Finally, the differences $H_{i,T}^*$ can also be efficiently evaluated using the following expressions:

$$H_{2,T}^*(x_l; x_h) = \frac{15e^{x_l}}{\pi^4} \sum_{k=3}^{N_T+1} \sum_{l=0}^{\lfloor \frac{k-2}{2} \rfloor} \binom{k-1}{2l+1} \frac{2T_{k-1}x_m^k \delta_m^{2l}}{(k-1)k},$$

$$H_{3,T}^*(x_l; x_h) = \frac{15e^{x_l}}{\pi^4} \sum_{k=3}^{N_T+1} \sum_{l=0}^{\lfloor \frac{k-3}{2} \rfloor} \binom{k-2}{2l+1} \frac{2T_{k-1}x_m^k \delta_m^{2l}}{(k-2)k},$$

$$H_{4,T}^*(x_l; x_h) = \frac{15e^{x_l}}{\pi^4} \left\{ \sum_{k=4}^{N_T+1} \sum_{l=0}^{\lfloor \frac{k-4}{2} \rfloor} \binom{k-3}{2l+1} \frac{2T_{k-1}x_m^k \delta_m^{2l}}{(k-3)k} + \frac{T_2 x_m^3}{3\delta_m} \ln \left( \frac{1+\delta_m}{1-\delta_m} \right) \right\},$$

where $\binom{N}{k}$ denotes the binomial coefficient "N choose k", and the function $\lfloor \cdot \rfloor$ is the floor function. For small values of $\delta_m$, the logarithm in $H_{4,T}$ can be approximated by:

$$\frac{1}{\delta_m} \ln \left( \frac{1 + \delta_m}{1 - \delta_m} \right) \approx 2 \left( 1 + \frac{\delta_m^2}{3} \right).$$

### ❖ Computation of combined expressions of the functions $H_4$, $H_3$, and $H_2$

Let $x_1$ and $x_2$ be two real numbers such that $x_2 \ge x_1$, and let $x_0 = (x_1 + x_2)/2$ and $\delta = (x_2 - x_1)/(2x_0)$. In the following appendix, we will need to evaluate the following combinations:

$$\mathcal{E}(x_1; x_2) = x_0^3 \mathcal{H}_4(x_1; x_2), \qquad (A.12)$$
$$ \mathcal{F}(x_1; x_2) = x_0^2 \mathcal{H}_3(x_1; x_2) - x_0^3 \mathcal{H}_4(x_1; x_2) , \qquad (A.13) $$

$$ \mathcal{P}(x_1; x_2) = x_0 \mathcal{H}_2(x_1; x_2) - 2x_0^2 \mathcal{H}_3(x_1; x_2) + x_0^3 \mathcal{H}_4(x_1; x_2) . \qquad (A.14) $$

In practice, these quantities are evaluated using the following expressions, valid for any configuration of the bounds $x_1$ and $x_2$ relative to a threshold value $x^*$:

$$ \mathcal{E}(x_1; x_2) = \delta e^{-x_1} \mathcal{E}^*(x_1; x_2) , $$

$$ \mathcal{F}(x_1; x_2) = \delta^3 e^{-x_1} \mathcal{F}^*(x_1; x_2) , $$

$$ \mathcal{P}(x_1; x_2) = \delta^3 e^{-x_1} \mathcal{P}^*(x_1; x_2) . $$

with:

$$ \mathcal{E}^*(x_1; x_2) = \begin{cases} \mathcal{H}_{4,P}^*(x_1; x_2) & \text{if } x_1 > x^* \text{ and } x_2 > x^* \\ \mathcal{H}_{4,T}^*(x_1; x_2) & \text{if } x_1 \le x^* \text{ and } x_2 \le x^* \\ \mathcal{H}_4^*(x_1; x_2) & \text{otherwise} \end{cases} , $$

$$ \mathcal{F}^*(x_1; x_2) = \begin{cases} \mathcal{F}_P^*(x_1; x_2) & \text{if } x_1 > x^* \text{ and } x_2 > x^* \\ \mathcal{F}_T^*(x_1; x_2) & \text{if } x_1 \le x^* \text{ and } x_2 \le x^* \\ \delta^{-2} (\mathcal{H}_3^*(x_1; x_2) - \mathcal{H}_4^*(x_1; x_2)) & \text{otherwise} \end{cases} , $$

$$ \mathcal{P}^*(x_1; x_2) = \begin{cases} \mathcal{P}_P^*(x_1; x_2) & \text{if } x_1 > x^* \text{ and } x_2 > x^* \\ \mathcal{P}_T^*(x_1; x_2) & \text{if } x_1 \le x^* \text{ and } x_2 \le x^* \\ \delta^{-2} (\mathcal{H}_2^*(x_1; x_2) - 2\mathcal{H}_3^*(x_1; x_2) + \mathcal{H}_4^*(x_1; x_2)) & \text{otherwise} \end{cases} . $$

Thus, the integral $\mathcal{E}^*$ is already known, since it is computed in the same way as $\mathcal{H}_4^*$. We still need to derive explicit the expressions of $\mathcal{F}^*$ and $\mathcal{P}^*$. The quantities $\mathcal{F}_P^*, \mathcal{F}_T^*, \mathcal{P}_P^*, \text{ and } \mathcal{P}_T^*$ are defined from the functions $\mathcal{H}_{i,P}^*$ and $\mathcal{H}_{i,T}^*$ as follows:

$$ \mathcal{F}_P^*(x_1; x_2) = \mathcal{H}_{2,P}^*(x_1; x_2) - \mathcal{H}_{4,P}^*(x_1; x_2) , $$

$$ \mathcal{F}_T^*(x_1; x_2) = \mathcal{H}_{3,T}^*(x_1; x_2) - \mathcal{H}_{4,T}^*(x_1; x_2) , $$

$$ \mathcal{P}_P^*(x_1; x_2) = \mathcal{H}_{2,P}^*(x_1; x_2) - 2\mathcal{H}_{2,P}^*(x_1; x_2) + \mathcal{H}_{4,P}^*(x_1; x_2) , $$

$$ \mathcal{P}_T^*(x_1; x_2) = \mathcal{H}_{2,T}^*(x_1; x_2) - 2\mathcal{H}_{3,T}^*(x_1; x_2) + \mathcal{H}_{4,T}^*(x_1; x_2) . $$

From the explicit expressions of the functions $\mathcal{H}_{i,P}^*$ given earlier, one obtains the following efficient formulas for the computation of $\mathcal{F}_P^*$ and $\mathcal{P}_P^*$:

$$ \mathcal{F}_P^*(x_1; x_2) = \frac{15}{\pi^4} \sum_{k=1}^{N_P} e^{-(k-1)x_1} \left[ \frac{z_k^1(x_0, \delta)}{k^4} + x_0 \frac{z_k^2(x_0, \delta)}{k^3} + x_0^2 \frac{z_k^3(x_0, \delta)}{k^2} \right] - \frac{8e^{x_1}}{3(1 - \delta^2)^3} , $$
$$ \mathcal{P}_P^*(x_1; x_2) = -\frac{15}{\pi^4} \sum_{k=1}^{N_P} e^{-(k-1)x_1} \left[ \frac{q_k^1(x_0, \delta)}{k^4} + x_0 \frac{q_k^2(x_0, \delta)}{k^3} + x_0^2 \frac{h_k^1(x_0, \delta)}{k^2} \right] + \frac{2e^{x_1}(1+3\delta^2)}{3(1-\delta^2)^3}, $$

where the auxiliary function $h_k^1$ was introduced in the previous section, and where the functions $z_k^i$ and $q_k^i$ are defined as follows:

$$ z_k^1(x_0, \delta) = \frac{8(\Sigma_k(x_0, \delta) + \Delta_k(x_0, \delta)) + (1-\delta^2)(1+3\delta^2)\mathcal{T}_k(x_0, \delta)}{(1-\delta^2)^3}, $$

$$ z_k^2(x_0, \delta) = \frac{(1 - \delta^2)\mathcal{T}_k(x_0, \delta) + 4\Delta_k(x_0, \delta)}{(1 - \delta^2)^2}, $$

$$ z_k^3(x_0, \delta) = \frac{\Delta_k(x_0, \delta)}{1 - \delta^2}, $$

$$ q_k^1(x_0, \delta) = 2 \frac{(1 + 3\delta^2)(\Delta_k(x_0, \delta) + \Sigma_k(x_0, \delta)) - (1 - \delta^2)(1 - 3\delta^2)\mathcal{T}_k(x_0, \delta)}{(1 - \delta^2)^3}, $$

$$ q_k^2(x_0, \delta) = \frac{4\delta^2\Delta_k(x_0, \delta)}{(1 - \delta^2)^2}. $$

The functions $\Sigma_k$ and $\Delta_k$ are given respectively by equations (A.11) and (A.10), and the function $\mathcal{T}_k$ is defined by:

$$ \mathcal{T}_k(x_0, \delta) = \frac{kx_0\Sigma_k(x_0, \delta) - \Delta_k(x_0, \delta)}{\delta^2}. \qquad (A.15) $$

For small values of $kx_0\delta$, this function admits the following approximation:

$$ \mathcal{T}_k(x_0, \delta) \approx \frac{2(kx_0)^3}{3} \left( 1 - kx_0\delta + \frac{3(kx_0\delta)^2}{5} \right). $$

Similarly, the explicit expressions of $\mathcal{H}_{i,T}^*$ make it possible to efficiently evaluate $\mathcal{F}_T^*$ and $\mathcal{P}_T^*$:

$$ \mathcal{F}_T^*(x_1; x_2) = \frac{15e^{x_1}}{\pi^4} \left\{ \sum_{k=2}^{\left\lfloor \frac{N_T}{2} \right\rfloor} \frac{2x_0^{2k+1} T_{2k} \delta^{2(k-2)}}{(2k-1)(2k+1)} + \sum_{k=6}^{N_T+1} \sum_{l=0}^{\left\lfloor \frac{k-6}{2} \right\rfloor} \binom{k-2}{2l+3} \frac{4x_0^k T_{k-1}(l+1)\delta^{2l}}{(k-3)(k-2)k} - \frac{T_2 x_0^3}{3\delta^3} \left[ \ln \left( \frac{1+\delta}{1-\delta} \right) - 2\delta \right] \right\}, $$

$$ \mathcal{P}_T^*(x_1; x_2) = \frac{15e^{x_1}}{\pi^4} \left\{ \sum_{k=2}^{\left\lfloor \frac{N_T+1}{2} \right\rfloor} \frac{x_0^{2k} T_{2k-1} \delta^{2(k-2)}}{(2k-1)k} + \sum_{k=2}^{\left\lfloor \frac{N_T}{2} \right\rfloor} \frac{2(2k-3)x_0^{2k+1} T_{2k} \delta^{2(k-2)}}{(2k-1)(2k+1)} + \sum_{k=6}^{N_T+1} \sum_{l=0}^{\left\lfloor \frac{k-6}{2} \right\rfloor} \binom{k-1}{2l+3} \frac{4(l+1)(2l+1)x_0^k T_{k-1} \delta^{2l}}{(k-3)(k-2)(k-1)k} + \frac{T_2 x_0^3}{3\delta^3} \left[ \ln \left( \frac{1+\delta}{1-\delta} \right) - 2\delta \right] \right\}. $$
For small values of $\delta$, the logarithmic term is approximated by:

$$ \frac{1}{\delta^3} \left\{ \ln \left( \frac{1+\delta}{1-\delta} \right) - 2\delta \right\} \approx \frac{2}{3} \left( 1 + \frac{3\delta^2}{5} \right) . $$

The main advantage of the formulations developed in this section for the quantities $\mathcal{E}$, $\mathcal{F}$, and $\mathcal{P}$ lies in the fact that they explicitly highlight the proportionality factors in $\delta$ and $e^{-x_1}$. This explicit factorization facilitates the computation of these quantities, and in particular that of the ratios $\mathcal{F}/\mathcal{E}$ and $\mathcal{P}/\mathcal{E}$, which will be used in Appendix B. *Numerical calculation of the radiative quantities.*

However, in the case where $x_1 < x^*$ and $x_2 > x^*$, no such explicit dependence on $\delta$ can be extracted for $\mathcal{F}$, $\mathcal{P}$, and $\mathcal{E}$. This limitation then degrades the accuracy of the computations of the ratios $\mathcal{F}/\mathcal{E}$ and $\mathcal{P}/\mathcal{E}$ when $\delta$ becomes very small.

## ❖ Derivatives of the functions $H_4$, $H_3$, and $H_2$

Let us now present the odd-order derivatives of the functions $H_2$, $H_3$, and $H_4$. The first, third, and fifth derivatives of the function $H_2$ are written as follows:

$$ H'_2(x) = \frac{\mathcal{B}(x)}{x^2}, $$

$$ H_2^{(3)}(x) = \frac{h_2^{(3)}(x)}{x^4} = \frac{x^2 b'(x) - 4 x b(x) + 6 \mathcal{B}(x)}{x^4}, $$

$$ H_2^{(5)}(x) = \frac{h_2^{(5)}}{x^6} = \frac{x^4 b^{(3)}(x) - 8 x^3 b''(x) + 36 x^2 b'(x) - 96 x b(x) + 120 \mathcal{B}(x)}{x^6}, $$

those of the function $H_3$ are given by:

$$ H'_3(x) = \frac{\mathcal{B}(x)}{x^3}, $$

$$ H_3^{(3)}(x) = \frac{h_3^{(3)}}{x^5} = \frac{x^2 b'(x) - 6 x b(x) + 12 \mathcal{B}(x)}{x^5}, $$

$$ H_3^{(5)}(x) = \frac{h_3^{(5)}}{x^7} = \frac{x^4 b^{(3)}(x) - 12 x^3 b''(x) + 72 x^2 b'(x) - 240 x b(x) + 360 \mathcal{B}(x)}{x^7}, $$

while those of the function $H_4$ are written as follows:

$$ H'_4(x) = \frac{\mathcal{B}(x)}{x^4}, $$

$$ H_4^{(3)}(x) = \frac{h_4^{(3)}(x)}{x^6} = \frac{x^2 b'(x) - 8 x b(x) + 20 \mathcal{B}(x)}{x^6}, $$
$$ \mathcal{H}_4^{(5)}(x) = \frac{h_4^{(5)}}{x^8} = \frac{x^4 b^{(3)}(x) - 16x^3 b''(x) + 120x^2 b'(x) - 480 x b(x) + 840 \mathcal{B}(x)}{x^8} . $$

In the expressions above, the functions $h_i^{(j)}$ denote the numerators of the corresponding derivatives, defined by the general relation $h_i^{(j)}(x) = x^{i+j-1}\mathcal{H}_i^{(j)}(x)$.
<figure><img src="image_81.png" /></figure>
<figure><img src="image_82.png" /></figure>

## NUMERICAL CALCULATION OF THE RADIATIVE QUANTITIES

As presented in Chapter 2. *Radiative hydrodynamics*, the computation of the integrals associated with the radiation energy, radiation flux, and radiation pressure is essential in the framework of the M1-multigroup model, as it enables a precise evaluation of its closure relation. These integrals make it possible to compute the reduced flux $f_g$ and, more importantly, the Eddington factor $\chi_g$, from the Lagrange multipliers $(\alpha_{0,g}, \beta_g)$. An accurate determination of $\chi_g$ relies on evaluating these integrals in search algorithms, described in Appendix C. *Computing the Eddington factor using search algorithms*. This approach is largely based on the work of Hung Chinh Nguyen (2011) [37], which I have adapted and, in some cases, improved in order to increase accuracy when the initial estimates proved insufficient.

The objective of this section is to present in detail the numerical computation of the integrals of the radiative quantities associated to group g: the energy $E_g$, the flux $\vec{F}_g$, and the pressure $P_g$. To simplify notation, we denote by $v_1$ and $v_2$ the frequency bounds of the considered group. By convention, the x-axis is chosen as the propagation direction of the radiative flux (see figure B.1). Let us first recall the expression of the specific intensity (2.62), which, in the framework of the M1-multigroup model, can be written as:

$$I_v(\theta; \alpha_{0,g}, \beta_g) = \sum_{g=1}^{G} \mathbb{1}_g(v) I_{v,g}(\theta; \alpha_{0,g}, \beta_g),$$
<figure><img src="image_83.png" /><figcaption>Figure B.1: Form of the specific intensity $I_{v,g}$ for different values of $\beta_g$, assuming $h\nu\alpha_{0,g}/k_B = 1$. The reference frame is chosen such that the x-axis corresponds to the direction of the radiative flux $\vec{F}_g$.</figcaption></figure>

where $I_{v,g}$ denotes the specific intensity within frequency group g, given by:

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

The parameters $\alpha_{0,g}$ and $\beta_g$ are the Lagrange multipliers associated with group g, with $\alpha_{0,g} \in \mathbb{R}^+$ and $\beta_g \in ]-1; 1[$. The angle $\theta$ denotes the angle between a considered direction and the propagation axis of the radiative flux (see figure B.1). In the adopted reference frame, the radiative flux and radiative pressure are expressed as follows:

$$\left\{ \begin{array}{l} \vec{F}_g = F_g \vec{x} = \operatorname{sgn}(F_g) f_g c E_g \vec{x}, \\ \mathbb{P}_g = \begin{bmatrix} \mathbb{P}_g & 0 \\ 0 & \mathbb{P}_g^* \end{bmatrix} = \begin{bmatrix} \chi_g & 0 \\ 0 & \frac{1-\chi_g}{2} \end{bmatrix} E_g, \end{array} \right.$$

where $f_g$ and $\chi_g$ denote the reduced flux and the Eddington factor, respectively. These two quantities vary within the following intervals: $f_g \in [0; 1[$ and $\chi_g \in [1/3; 1[$. As indicated in equations (2.68), (2.69), and (2.70), the computation of the radiative quantities can then be
written as:

$$
\begin{cases}
E_g = \frac{4h\pi}{c^3} \int_{v_1}^{v_2} \int_{-1}^{1} v^3 \left[ \exp\left(\frac{hv}{k_B}\alpha_{0,g}(1+\beta_g\mu)\right) - 1 \right]^{-1} d\mu \, dv , \\
F_g = \frac{4h\pi}{c^2} \int_{v_1}^{v_2} \int_{-1}^{1} v^3 \mu \left[ \exp\left(\frac{hv}{k_B}\alpha_{0,g}(1+\beta_g\mu)\right) - 1 \right]^{-1} d\mu \, dv , \\
P_g = \frac{4h\pi}{c^3} \int_{v_1}^{v_2} \int_{-1}^{1} v^3 \mu^2 \left[ \exp\left(\frac{hv}{k_B}\alpha_{0,g}(1+\beta_g\mu)\right) - 1 \right]^{-1} d\mu \, dv .
\end{cases}
$$

These expressions can also be rewritten in a dimensionless form by introducing the follow-
ing redimensionalized quantities:

$$
\left\{
\begin{aligned}
\hat{E}_g &= \frac{1}{2} \int_{-1}^{1} \frac{\mathcal{B}(x_1; x_2)}{\hat{\alpha}_0^4 (1 + \beta_g \mu)^4} \mathrm{d}\mu \, , \\
\hat{F}_g &= \frac{1}{2} \int_{-1}^{1} \mu \frac{\mathcal{B}(x_1; x_2)}{\hat{\alpha}_0^4 (1 + \beta_g \mu)^4} \mathrm{d}\mu \, , \\
\hat{P}_g &= \frac{1}{2} \int_{-1}^{1} \mu^2 \frac{\mathcal{B}(x_1; x_2)}{\hat{\alpha}_0^4 (1 + \beta_g \mu)^4} \mathrm{d}\mu \, .
\end{aligned}
\right.
$$

Here, the redimensionalized quantities are defined as follows:

$$
\hat{E}_g = E_g / a_R \hat{T}^4, \quad \hat{F}_g = F_g / a_R c \hat{T}^4, \quad \hat{P}_g = P_g / a_R \hat{T}^4, \quad \hat{\alpha}_0 = \alpha_{0,g} \hat{T}^4,
$$

where $\hat{T}$ denotes a reference temperature and $a_R$ the radiation constant. $\mathcal{B}$ is the integral of the Planck function $b$ evaluated between $x_1$ and $x_2$ (see Appendix A. Basic concepts on the Planck function). The bounds $x_1$ and $x_2$ of the integral, are given by $x_1 = \hat{v}_1\alpha_0(1+\beta_g\mu)$ and $x_2 = \hat{v}_2\alpha_0(1+\beta_g\mu)$, where $\hat{v} = \frac{hv}{(k_B\hat{T})}$. To simplify notations, we will hereafter omit the “$\hat{\cdot}$” when referring to dimensionless quantities.

Finally, introducing $z = v\alpha_0(1 + \beta_g\mu)$, the dimensionless radiative quantities can be expressed as follows:

$$
\left\{
\begin{aligned}
E_g &= E_{v2} - E_{v1}, \\
F_g &= F_{v2} - F_{v1}, \\
P_g &= P_{v2} - P_{v1},
\end{aligned}
\right.
$$

where :

$$
E_{v_j} = \frac{\nu_j^3}{2\alpha_0\beta_g} \int_{z_1(v_j)}^{z_2(v_j)} \frac{\mathcal{B}(z)}{z^4} \, dz\, ,
$$

$$
F_{v_j} = \frac{\nu_j^2}{2\alpha_0^2\beta_g^2} \int_{z_1(v_j)}^{z_2(v_j)} \frac{\mathcal{B}(z)}{z^3} \, dz - \frac{\nu_j^3}{2\alpha_0\beta_g^2} \int_{z_1(v_j)}^{z_2(v_j)} \frac{\mathcal{B}(z)}{z^4} \, dz\, ,
$$

$$
P_{v_j} = \frac{\nu_j}{2\alpha_0^3\beta_g^3} \int_{z_1(v_j)}^{z_2(v_j)} \frac{\mathcal{B}(z)}{z^2} \, dz - \frac{2\nu_j^2}{2\alpha_0^2\beta_g^3} \int_{z_1(v_j)}^{z_2(v_j)} \frac{\mathcal{B}(z)}{z^3} \, dz + \frac{\nu_j^3}{2\alpha_0\beta_g^3} \int_{z_1(v_j)}^{z_2(v_j)} \frac{\mathcal{B}(z)}{z^4} \, dz\, .
$$
and where $z_1(v) = \alpha_0 v(1 - \beta_g)$ and $z_2(v) = \alpha_0 v(1 + \beta_g)$. Here, one recognizes the functions $H_2$, $H_3$, and $H_4$, introduced in Appendix A.3. Functions derived from the Planck function. These terms can also be rewritten in the following form:

$$ 
\begin{aligned}
E_{v_j} &= \frac{\mathcal{E}(z_1(v_j); z_2(v_j))}{2\alpha_0^4\beta_g}, \\
F_{v_j} &= \frac{\mathcal{F}(z_1(v_j); z_2(v_j))}{2\alpha_0^4\beta_g^2}, \\
P_{v_j} &= \frac{\mathcal{P}(z_1(v_j); z_2(v_j))}{2\alpha_0^4\beta_g^3}.
\end{aligned}
$$

Thus, the reduced flux and the Eddington factor can be expressed as:

$$ E_g = \frac{\mathcal{E}(z_1(v_2); z_2(v_2)) - \mathcal{E}(z_1(v_1); z_2(v_1))}{2\alpha_0^4\beta_g}, \qquad (B.1) $$

$$ f_g = \left| \frac{1}{\beta_g} \frac{\mathcal{F}(z_1(v_2); z_2(v_2)) - \mathcal{F}(z_1(v_1); z_2(v_1))}{\mathcal{E}(z_1(v_2); z_2(v_2)) - \mathcal{E}(z_1(v_1); z_2(v_1))} \right|, \qquad (B.2) $$

$$ \chi_g = \frac{1}{\beta_g^2} \frac{\mathcal{P}(z_1(v_2); z_2(v_2)) - \mathcal{P}(z_1(v_1); z_2(v_1))}{\mathcal{E}(z_1(v_2); z_2(v_2)) - \mathcal{E}(z_1(v_1); z_2(v_1))}, \qquad (B.3) $$

where the functions $\mathcal{E}$, $\mathcal{F}$, and $\mathcal{P}$ are defined by equations (A.12), (A.13), and (A.14). Four cases can be distinguished for the computation of these quantities (see figure B.2):

1. **"General" case:** This corresponds to the non-isotropic case (i.e., $|\beta_g| > 0$) in which the frequency bounds of the group, $v_1$ and $v_2$, are well separated;

2. **Case close to isotropy ($\beta_g \ll 1$):** In this case, since $z_1(v) \approx z_2(v)$, the terms are generally well computed thanks to the expansions presented in Appendix A. Basic concepts on the Planck function, except when $z_1(v_1) < 1$ and $z_2(v_2) > 1$ and/or $z_1(v_2) < 1$ and $z_2(v_2) > 1^1$. We however assume that the frequency bounds $v_1$ and $v_2$ remain well separated;

3. **Case close to the spectral limit ($v_1 \approx v_2$):** In this situation, differences according to the frequencies bounds, such as $\mathcal{E}(z_1(v_2); z_2(v_2)) - \mathcal{E}(z_1(v_1); z_2(v_1))$, become very small and are likely to be numerically inaccurate;

4. **Case close to both isotropy and the spectral limit:** In this case, the difficulties of the two previous cases occur simultaneously.

<footnote n="1">¹The threshold value $x^*$ used here is 1.</footnote>
<figure><img src="image_84.png" /><figcaption>Figure B.2: Different cases for computing the quantities $E_g$, $f_g$, and $\chi_g$. $\delta_g$ denotes the group narrowness and depends on the frequency bounds of the group according to the expression $\delta_g = \nu_1/\nu_2$.</figcaption></figure>

In this work, I addressed only Cases 1 and 2, while Cases 3 and 4, which correspond to regimes close to the spectral limit, were not considered. This limitation required restricting the group narrowness $\delta_g$ to values smaller than 0.9999. I now present the method used to compute the quantities $E_g$, $f_g$, and $\chi_g$, first in a general framework and then under an approximation close to isotropy.

## B.1 “General” case

To accurately compute the quantities $E_g$, $f_g$, and $\chi_g$, we replace the functions $\mathcal{E}$, $\mathcal{F}$, and $\mathcal{P}$ by the functions $\mathcal{E}^*$, $\mathcal{F}^*$, and $\mathcal{P}^*$, defined by the equations (see Appendix A.3. Functions derived from the Planck function):

$$ \mathcal{E}^*(x_1; x_2) = \delta^{-1} e^{x_1} \mathcal{E}(x_1; x_2) , $$

$$ \mathcal{F}^*(x_1; x_2) = \delta^{-3} e^{x_1} \mathcal{F}(x_1; x_2) , $$

$$ \mathcal{P}^*(x_1; x_2) = \delta^{-3} e^{x_1} \mathcal{P}(x_1; x_2) . $$
The equations (B.1), (B.2), and (B.3) can then be rewritten as:

$$E_g = \frac{e^{-z_1(v_1)}}{2\alpha_0^4} [e^{-X_{1,2}\mathcal{E}^*}(z_1(v_2); z_2(v_2)) - \mathcal{E}^*(z_1(v_1); z_2(v_1))] , \quad (B.4)$$

$$f_g = \left| \beta_g \frac{e^{-X_{1,2}\mathcal{F}^*(z_1(v_2); z_2(v_2)) - \mathcal{F}^*(z_1(v_1); z_2(v_1))}}{e^{-X_{1,2}\mathcal{E}^*(z_1(v_2); z_2(v_2)) - \mathcal{E}^*(z_1(v_1); z_2(v_1))}} \right| , \quad (B.5)$$

$$\chi_g = \frac{e^{-X_{1,2}\mathcal{P}^*(z_1(v_2); z_2(v_2)) - \mathcal{P}^*(z_1(v_1); z_2(v_1))}}{e^{-X_{1,2}\mathcal{E}^*(z_1(v_2); z_2(v_2)) - \mathcal{E}^*(z_1(v_1); z_2(v_1))}} , \quad (B.6)$$

where $X_{1,2} = \alpha_0(v_2 - v_1)(1 - \beta_g)$. For the computation of the functions $\mathcal{E}^*$, $\mathcal{F}^*$, and $\mathcal{P}^*$, we rely on the expressions developed in Appendix A. Basic concepts on the Planck function. This formulation offers several numerical advantages: it avoids any division by $\beta_g$, thereby eliminating the risk of divergence or error amplification when $\beta_g$ becomes small; it also prevents the multiplication of terms by decaying exponentials, which are numerically rounded to zero when $\alpha_0\nu$ is large, making the calculations infeasible. However, in cases where $z_1(v_1) < 1 < z_2(v_1)$ or $z_1(v_2) < 1 < z_2(v_2)$, the expressions used inevitably involve divisions by $\beta_g$. These configurations lead to amplified numerical errors and a significant loss of accuracy as $\beta_g$ tends to zero, as illustrated in figures B.3a and B.3b. To overcome this limitation, it is necessary to resort to an alternative method that is better suited to this regime.

## B.2 Case close to isotropy

In the near-isotropic case, $\beta_g$ takes very small values, and thus one can compute the Taylor expansion of the radiation energy, the reduced flux, and the Eddington factor with respect to $\beta_g$. Introducing $\xi_j = \alpha_0\nu_j$, the Taylor expansions of the integral functions $\mathcal{H}_2$, $\mathcal{H}_3$, and $\mathcal{H}_4$ are:

$$\mathcal{H}_2(z_1(v_j); z_2(v_j)) = \frac{2 \beta_g}{\xi_j} \mathcal{B}(\xi_j) + \frac{\beta_g^3}{3\xi_j} h_2^{(3)}(\xi_j) + \frac{\beta_g^5}{60\xi_j} h_2^{(5)}(\xi_j) + O(\beta_g^6) ,$$

$$\mathcal{H}_3(z_1(v_j); z_2(v_j)) = \frac{2 \beta_g}{\xi_j^2} \mathcal{B}(\xi_j) + \frac{\beta_g^3}{3\xi_j^2} h_3^{(3)}(\xi_j) + \frac{\beta_g^5}{60\xi_j^2} h_3^{(5)}(\xi_j) + O(\beta_g^6) ,$$

$$\mathcal{H}_4(z_1(v_j); z_2(v_j)) = \frac{2 \beta_g}{\xi_j^3} \mathcal{B}(\xi_j) + \frac{\beta_g^3}{3\xi_j^3} h_4^{(3)}(\xi_j) + \frac{\beta_g^5}{60\xi_j^3} h_4^{(5)}(\xi_j) + O(\beta_g^6) ,$$

where the functions $h_i^{(j)}$ are defined in Appendix A.3. Functions derived from the Planck function. Thus, the Taylor expansions of the functions $\mathcal{E}$, $\mathcal{F}$, and $\mathcal{P}$ can be written as:

$$\mathcal{E}(z_1(v_j); z_2(v_j)) = 2 \beta_g \mathcal{B}(\xi_j) + \frac{\beta_g^3}{3} h_4^{(3)}(\xi_j) + \frac{\beta_g^5}{60} h_4^{(5)}(\xi_j) + O(\beta_g^6) ,$$
<figure><img src="image_85.png" /><figcaption>(a) Evolution of the reduced flux without using the Taylor expansions for small values of $|\beta_g|$.</figcaption></figure>

<figure><img src="image_86.png" /><figcaption>(b) Evolution of the Eddington factor without using the Taylor expansions for small values of $|\beta_g|$.</figcaption></figure>

<figure><img src="image_87.png" /><figcaption>(c) Evolution of the reduced flux using the Taylor expansions for small values of $|\beta_g|$.</figcaption></figure>

<figure><img src="image_88.png" /><figcaption>(d) Evolution of the Eddington factor using the Taylor expansions for small values of $|\beta_g|$.</figcaption></figure>

Figure B.3: Evolution of the reduced flux and the Eddington factor as functions of $|\beta_g|$, with and without the use of the Taylor expansions presented in Section B.2. Case close to isotropy, for different values of $\alpha_0$ and using the frequency bounds $v_1 = 0.1$ and $v_2 = 10$.

$$ \mathcal{F}(z_1(v_j); z_2(v_j)) = \frac{\beta_g^3}{3} f_3(\xi_j) + \frac{\beta_g^5}{60} f_5(\xi_j) + \mathcal{O}(\beta_g^6) , $$

$$ \mathcal{P}(z_1(v_j); z_2(v_j)) = \frac{2\beta_g^3}{3} \mathcal{B}(\xi_j) + \frac{\beta_g^5}{60} p_5(\xi_j) + \mathcal{O}(\beta_g^6) , $$

where the auxiliary functions $f_i$ and $p_i$ are given by:

$$ f_3(x) = h_4^{(3)}(x) - h_3^{(3)}(x) = 2xb(x) - 8\mathcal{B}(x) , $$

$$ f_5(x) = h_4^{(5)}(x) - h_3^{(5)}(x) = 4[x^3b''(x) - 12x^2b'(x) + 60xb(x) - 120\mathcal{B}(x)] , $$

$$ p_5(x) = h_2^{(5)}(x) - 2h_3^{(5)}(x) + h_4^{(5)}(x) = 12[x^2b'(x) - 8xb(x) + 20\mathcal{B}(x)] , $$

and the successive derivatives of the Planck function $b$ are given in Appendix A.1. *Derivatives of the Planck function.* Thus, a second-order Taylor expansion of the quantities of interest $E_g$,
$f_g$, and $\chi_g$ can be derived:

$$E_g = \frac{1}{\alpha_0^4} \left\{ \mathcal{B}(\xi_2) - \mathcal{B}(\xi_1) + \frac{\beta_g^2}{6} (h_4^{(3)}(\xi_2) - h_4^{(3)}(\xi_1)) \right\} \quad (B.7)$$

$$f_g = \left| \frac{\beta_g f_3(\xi_2) - f_3(\xi_1)}{6 \mathcal{B}(\xi_2) - \mathcal{B}(\xi_1)} \right| \quad (B.8)$$

$$\chi_g = \frac{1}{3} \left\{ 1 + \frac{\beta_g^2}{120} \frac{3[p_5(\xi_2) - p_5(\xi_1)] - 20[h_4^{(3)}(\xi_2) - h_4^{(3)}(\xi_1)]}{\mathcal{B}(\xi_2) - \mathcal{B}(\xi_1)} \right\} \quad (B.9)$$

As illustrated in figure B.3, using the Taylor expansions prevents the amplification of numerical errors when $|\beta_g|$ is small. In practice, these expressions are therefore used to evaluate the radiative quantities $E_g$, $f_g$, and $\chi_g$ in this regime.
C

COMPUTING THE EDDINGTON FACTOR USING SEARCH ALGORITHMS

In a simulation performed with the HADES code, based on the M1-multigroup model, each spectral group *g* is characterized at a given time by two fundamental quantities: the reduced flux *f*<sub>*g*</sub> and the radiative energy *E*<sub>*g*</sub>. One of the main objectives of the model is to determine, for each group, the Eddington factor *χ*<sub>*g*</sub>, which closes the radiative system by relating the radiative pressure tensor to the radiative energy. The three quantities, *E*<sub>*g*</sub>, *f*<sub>*g*</sub>, and *χ*<sub>*g*</sub>, are not independent. Indeed, they can all be expressed in terms of two fundamental parameters, denoted *α*<sub>0,*g*</sub> and *β*<sub>*g*</sub>, known as Lagrange multipliers. These multipliers naturally arise within a variational framework that consists in minimizing the entropy of the radiation field under constraints (energy and flux conservation), in accordance with the minimum entropy principle. The associated numerical procedure, which allows one to evaluate the integral expressions involving these multipliers, is detailed in Appendices A. *Basic concepts on the Planck function* and B. *Numerical calculation of the radiative quantities*, where the spectral integration methods and the approximations employed are presented.

According to the work of Turpault (2003) [97], it has been shown that for any physically
admissible pair (E_g, f_g), there exists a unique pair of Lagrange multipliers (α_0,g, β_g) that allows
one to reconstruct the radiation field compatible with these constraints. It is therefore natural to
devise a search algorithm aimed at recovering this unique pair of multipliers, and subsequently
deducing the corresponding Eddington factor χ_g. In this section, we assume that the frequency
bounds of group g are fixed at v_1 and v_2. The integral expressions depend explicitly on
these bounds as well as on the multipliers. For clarity, we introduce the following functional
notations:

• $E_g (\alpha_{0,g}, \beta_g)$: radiation energy associated to group $g$,

• $f_g (\alpha_{0,g}, \beta_g)$: reduced flux associated to group $g$,

• $\chi_g (\alpha_{0,g}, \beta_g)$: corresponding Eddington factor.

These functions are defined from the spectral integrals of the radiation field and are
evaluated numerically using the methods described in Appendix B. *Numerical calculation of
the radiative quantities*. Finally, we denote by $E_g^*$ and $f_g^*$ the target values of the radiation
energy and the reduced flux, obtained from the simulation, for which we seek to determine
the corresponding Eddington factor $\chi_g$. The search for the Lagrange multipliers therefore
constitutes an essential intermediate step in the consistent calculation of this factor.

C.1 Search algorithms

Historically, a first algorithm was proposed by Hung Chinh Nguyen in his PhD thesis (2011) [37], in the context of the development of the HADES code. This algorithm relies on a hybrid approach that combines the robustness of the bisection method with the fast convergence of Newton's method. More precisely, the bisection method is used to determine a suitable value of the Lagrange multiplier $\beta_g$, ensuring that the radiative flux constraint is satisfied, while Newton's method is then employed to efficiently compute the multiplier $\alpha_{0,g}$, associated with energy conservation. The entire procedure is summarized in the algorithm presented below.
Algorithm C.1 : Bissection-Newton algorithm

We seek to determine the pair $(x, \beta)$, with $x = \ln(\alpha_0)$, such that $\ln(E_g(e^x, \beta)) = \ln(E_g^*)$ and $f_g(e^x, \beta) = f_g^*$. We denote by $i$ the index of the main iteration, by $j$ the index of the dissection or Newton method iterations, and by $x^{(i)}$ and $\beta^{(i)}$ the multipliers estimated at iteration $i$. The procedure is as follows:

* **Special case:** If $f_g^* = 0$, then the Eddington factor is directly $\chi_g = 1/3$,

* **Otherwise,** we proceed through successive iterations:

1. **Initialization:** the multipliers are initialized from the values obtained with the M1-gray model (see Dubroca and Feugeas, 1999 [87]).

2. **Iterative loop:**

- **Step 1 - Bisection method for $\beta^{(i+1)}$:**

Fixing $x^{(i)}$, we search for $\beta^{(i+1)}$ such that $f_g(e^{x^{(i)}}, \beta^{(i+1)}) = f_g^*$. The search is carried out using the bisection method over the initial interval $[0; 1[$. The search stops when the condition

$$ \left| \left( \beta_{j+1}^{(i+1)} - \beta_j^{(i+1)} \right) / \beta_j^{(i+1)} \right| \le \epsilon_D $$

- **Step 2 - Newton search for $x^{(i+1)}$:**

Fixing $\beta^{(i+1)}$, we search for $x^{(i+1)}$ such that $\ln(E_g(\alpha_0^{(i+1)}, \beta^{(i+1)})) = \ln(E_g^*)$, where $\alpha_0^{(i+1)} = e^{x^{(i+1)}}$. A Newton algorithm in $\mathbb{R}$ is used. The search stops when is met the condition

$$ \left| \left( \alpha_{0,j+1}^{(i+1)} - \alpha_{0,j}^{(i+1)} \right) / \alpha_{0,j}^{(i+1)} \right| \le \epsilon_N $$

- **Step 3 - Stopping criterion:**

Steps 1 and 2 are repeated as long as the sum of the relative differences

$$ \left| \left( \alpha_0^{(i+1)} - \alpha_0^{(i)} \right) / \alpha_0^{(i)} \right| + \left| \left( \beta^{(i+1)} - \beta^{(i)} \right) / \beta^{(i)} \right| $$

remains greater than $\epsilon_T$.

* Once convergence is reached, the corresponding Eddington factor $\chi_g$ is computed.

The parameters $\epsilon_D$, $\epsilon_N$, and $\epsilon_T$ are small constants setting the tolerances for the steps of the bisection method, Newton's method, and of the global algorithm, respectively. In this work, they have been set to $10^{-8}$, $10^{-8}$, and $10^{-5}$.
The initial search algorithm, although functional, exhibits a notable limitation: it determines the Lagrange multipliers sequentially, estimating them one at a time. This approach slows down the computation and reduces its efficiency. To address this issue, I developed a new algorithm based on a *line search* strategy, capable of estimating both multipliers, $\alpha_{0,g}$ and $\beta_g$, simultaneously. By exploiting the structure of the problem and evaluating the constraints jointly, this method significantly improves the convergence speed while maintaining good accuracy on the associated radiative quantities.
Algorithm C.2 : Line search algorithm

We seek to determine the pair $(x, y)$, with $x = \ln(\alpha_0)$ and $y = y_0 \sigma^{-1}(\beta)$, such that the relations $\ln(E_g(e^x, \sigma(y/y_0))) = \ln(E_g^*)$ and $\ln(f_g(e^x, \sigma(y/y_0))) = \ln(f_g^*)$ are satisfied. We denote by $i$ the index of the main iteration, and by $x^{(i)}$ and $y^{(i)}$ the values of $x$ and $y$ estimated at iteration $i$. The procedure is as follows:

* **Special case:** if $f_g^* = 0$, then the Eddington factor is directly $\chi_g = 1/3$,

* **Otherwise,** we proceed through successive iterations:

1. **Initialization:** the multipliers are initialized from the values obtained with the M1-gray model (see Dubroca and Feugeas, 1999 [87]).

2. **Iterative loop:** we search for $x^{(i+1)}$ and $y^{(i+1)}$ using a line search algorithm, such that the sum of squared logarithmic residuals

$$ \left[ \ln(f_g(\alpha_0^{(i+1)}, \beta^{(i+1)})) - \ln(f_g^*) \right]^2 + \left[ \ln(E_g(\alpha_0^{(i+1)}, \beta^{(i+1)})) - \ln(E_g^*) \right]^2 $$

vanishes, where $\alpha_0^{(i+1)} = e^{x^{(i+1)}}$ and $\beta^{(i+1)} = \sigma(y^{(i+1)}/y_0)$.

3. **Stopping criterion:** the algorithm stops when the sum of the relative variations

$$ \left| \left( \frac{\alpha_{0,g}^{(i+1)} - \alpha_{0,g}^{(i)}}{\alpha_{0,g}^{(i)}} \right) + \left| \left( \frac{\beta_{g}^{(i+1)} - \beta_{g}^{(i)}}{\beta_{g}^{(i)}} \right) \right| $$

becomes smaller than a small constant $\epsilon$.

* Once convergence is achieved, the corresponding Eddington factor $\chi_g$ is computed.

The parameter $\epsilon$ is a small constant that sets the tolerance of the global algorithm, while $y_0$ is a normalization parameter, here set respectively to $10^{-5}$ and 100. The function $\sigma$ denotes the sigmoid function.

## C.2 Comparison of the search algorithms

The figure C.1 illustrates the computation time required to evaluate the Eddington factor as a function of the reduced flux $f_g$ and the dimensionless radiative temperature $T_g^1$. A precise

<footnote n="1">¹The dimensionless radiative temperature $T_g$ is related to the radiation energy $E_g$; its precise definition is given in Section 2.2.3. *The Eddington factor.*</footnote>
<figure><img src="image_89.png" /><figcaption>(a) Computation time of the Eddington factor using Algorithm C.1: Bisection–Newton algorithm.</figcaption></figure>

<caption>(b) Computation time of the Eddington factor using Algorithm C.2: Line-search algorithm.</caption>

Figure C.1: Computation time of the Eddington factor as a function of the reduced flux $f_g$ and the dimensionless radiative temperature $T_g$, for the rescaled frequency bounds $v_1 = 0.1$ and $v_2 = 10$. The results are obtained using the search algorithms C.1 and C.2. The gray areas indicate cases in which the algorithm did not converge.

analysis of this figure raises two main observations:

1. The **line search** algorithm yields an overall shorter computation time than the **Bisection–Newton** algorithm, while providing comparable estimates of the Eddington factor;

2. The **line search** algorithm also proves more robust: it converges across a wide range of conditions, unlike the Bisection–Newton method, whose convergence may fail in certain regions of the parameter space. The convergence failures of the line search are mainly located in the region where $f_g \approx 1$ and $T_g \gtrsim 1$. However, in this specific region, the radiative quantities can be evaluated using the analytical approximations given by equations (2.76), (2.77), and (2.78), detailed in Section 2.2.3. *The Eddington factor*, which makes it possible to circumvent these numerical difficulties.

Owing to its superior performance, both in numerical efficiency and reliability, the line search algorithm was therefore adopted for all computations presented in this work.
HYDRODYNAMIC SHOCK ON A WALL

In this appendix, I develop the analytical solution of the hydrodynamic case used as a reference in chapter 4. *Study of radiative shocks*. The physical configuration of the problem is recalled in figure D.1.

We first determine the hydrodynamic quantities in the downstream region from the upstream conditions. To this end, we use the Rankine–Hugoniot jump relations. In the particular case

<figure><img src="image_90.png" /><figcaption>Figure D.1: Diagram representing the structure of a hydrodynamic shock propagating from a wall, located in x = 0 at a given moment.</figcaption></figure>
where the fluid velocity in the downstream region is zero ($v_{down} = 0$), these relations are written:

$$ (\rho_{up} - \rho_{down})v_{shock} = \rho_{up}v_{up}, \qquad (D.1) $$

$$ \rho_{up}v_{shock}v_{up} = \rho_{up}v_{up}^2 + p_{up} - p_{down}, \qquad (D.2) $$

$$ (E_{up} - E_{down})v_{shock} = (E_{up} + p_{up})v_{up}, \qquad (D.3) $$

where E denotes the total energy of the fluid, given by the expression:

$$ E = \frac{1}{2}\rho v^2 + \frac{p}{\gamma - 1}, \qquad (D.4) $$

The resolution of this system makes it possible to express the velocity of the shock $v_{shock}$, as well as the density and pressure in the downstream medium, as a function of the upstream quantities. We obtain:

$$ v_{shock} = -\frac{(\gamma - 3)M_{lab} + \sqrt{\{(\gamma + 1)M_{lab}\}^2 + 16}}{4M_{lab}}v_{up}, \qquad (D.5) $$

$$ \rho_{down} = \frac{(\gamma + 1)M_{lab} + \sqrt{\{(\gamma + 1)M_{lab}\}^2 + 16}}{(\gamma - 3)M_{lab} + \sqrt{\{(\gamma + 1)M_{lab}\}^2 + 16}}\rho_{up}, \qquad (D.6) $$

$$ p_{down} = \frac{4 + \gamma M_{lab} \left\{ (\gamma + 1)M_{lab} + \sqrt{\{(\gamma + 1)M_{lab}\}^2 + 16} \right\}}{4} p_{up}, \qquad (D.7) $$

where the parameter $M_{lab}$ (analogous to a Mach number, but expressed in the laboratory reference frame) is defined by:

$$ M_{lab} = \frac{|v_{up}|}{\sqrt{\gamma p_{up} / \rho_{up}}} \qquad (D.8) $$

This number should not, however, be confused with the Mach number of the shock itself, which characterizes the shock strength and is given by:

$$ M = \frac{|v_{up} - v_{shock}|}{\sqrt{\gamma p_{up} / \rho_{up}}}, \qquad (D.9) $$

Assuming that the fluid behaves as an ideal gas, the downstream temperature can also be expressed as a function of the upstream temperature:

$$ T_{down} = \frac{p_{down}}{p_{up}} \frac{\rho_{up}}{\rho_{down}} T_{up} \qquad (D.10) $$
Finally, denoting by $t_0$ the initial time of the simulation, and by $\mathcal{H}$ the Heaviside function (equal to 1 for $x \ge 0$ and 0 otherwise), the full analytical solution of this problem is given by:

$$ \begin{cases} \rho = \rho_{down} + (\rho_{up} - \rho_{down})\mathcal{H}[x - v_{shock}(t - t_0)] , \\ v = v_{down} + (v_{up} - v_{down})\mathcal{H}[x - v_{shock}(t - t_0)] , \\ T = T_{down} + (T_{up} - T_{down})\mathcal{H}[x - v_{shock}(t - t_0)] . \end{cases} \qquad (D.11) $$
<figure><img src="image_91.png" /></figure>
BIBLIOGRAPHY

[1] A. Einstein. "Erklärung der Perihelbewegung des Merkur aus der allgemeinen Relativitätstheorie". In: *Sitzungsberichte der Königlich Preussischen Akademie der Wissenschaften* (Jan. 1915), pp. 831–839.

[2] B. P. Abbott et al. "Observation of Gravitational Waves from a Binary Black Hole Merger". In: *Phys. Rev. Lett.* 116 (6 Feb. 2016), p. 061102.

[3] D. D. Ryutov, R. P. Drake, and B. A. Remington. "Criteria for Scaled Laboratory Simulations of Astrophysical MHD Phenomena". In: *The Astrophysical Journal Supplement Series* 127.2 (Apr. 2000), p. 465.

[4] S. Bouquet et al. "From lasers to the universe: Scaling laws in laboratory astrophysics". In: *High Energy Density Physics* 6.4 (2010), pp. 368–380.

[5] É. Falize, C. Michaut, and S. Bouquet. "Similarity properties and scaling laws of radiation hydrodynamic flows in laboratory astrophysics". In: *The Astrophysical Journal* 730.2 (Mar. 2011), p. 96.

[6] V. Tranchant et al. "Generalizing the similitude approach for laboratory astrophysics through equivalence symmetries for the example of radiative waves". In: *Scientific Reports* 15.1 (Mar. 2025), p. 10806.

[7] Chiavassa, A. et al. "Radiative hydrodynamics simulations of red supergiant stars - IV. Gray versus non-gray opacities". In: *Astronomy & Astrophysics* 535 (2011), A22.

[8] Kravchenko, K. et al. "Tomography of cool giant and supergiant star atmospheres - I. Validation of the method". In: *Astronomy & Astrophysics* 610 (2018), A29.

[9] C.-Y. Wang and R. A. Chevalier. "Instabilities and Clumping in Type Ia Supernova Remnants". In: *The Astrophysical Journal* 549.2 (Mar. 2001), p. 1119.

[10] J. Minière et al. "Numerical study of the Vishniac instability in cooled supernova remnants". In: *Astronomy & Astrophysics* 617 (2018), A133.
[11] H. Abu-Shawareb et al. "Lawson Criterion for Ignition Exceeded in an Inertial Fusion Experiment". In: *Phys. Rev. Lett.* 129 (7 Aug. 2022), p. 075001.

[12] A. B. Zylstra et al. "Experimental achievement and signatures of ignition at the National Ignition Facility". In: *Phys. Rev. E* 106 (2 Aug. 2022), p. 025202.

[13] Remington, Bruce A. "Experiments in ICF, materials science, and astrophysics". In: *EPJ Web of Conferences* 59 (2013), p. 01013.

[14] A. Casner et al. "From ICF to laboratory astrophysics: ablative and classical Rayleigh-Taylor instability experiments in turbulent-like regimes". In: *Nuclear Fusion* 59.3 (Dec. 2018), p. 032002.

[15] V. A. Smalyuk et al. "Review of hydrodynamic instability experiments in inertially confined fusion implosions on National Ignition Facility". In: *Plasma Physics and Controlled Fusion* 62.1 (Oct. 2019), p. 014007.

[16] S. Wolf. "MC3D-3D continuum radiative transfer, Version 2". In: *Computer Physics Communications* 150.2 (2003), pp. 99-115.

[17] Pascucci, I. et al. "The 2D continuum radiative transfer problem - Benchmark results for disk configurations". In: *Astronomy & Astrophysics* 417.3 (2004), pp. 793-805.

[18] A. Maselli, A. Ferrara, and B. Ciardi. "CRASH: a radiative transfer scheme". In: *Monthly Notices of the Royal Astronomical Society* 345.2 (Oct. 2003), pp. 379-394.

[19] A. Smith et al. "AREPO-MCRT: Monte Carlo Radiation Hydrodynamics on a Moving Mesh". In: *The Astrophysical Journal* 905.1 (Dec. 2020), p. 27.

[20] P. Coelho. "The role of ray effects and false scattering on the accuracy of the standard and modified discrete ordinates methods". In: *Journal of Quantitative Spectroscopy and Radiative Transfer* 73.2 (2002). Third International Symposium on Radiative Transfer, pp. 231-238.

[21] Y.-F. Jiang, J. M. Stone, and S. W. Davis. "An algorithm for radiation magnetohydrodynamics based on solving the time-dependent transfer equation". In: *The Astrophysical Journal Supplement Series* 213.1 (June 2014), p. 7.

[22] Y.-F. Jiang. "An Implicit Finite Volume Scheme to Solve the Time-dependent Radiation Transport Equation Based on Discrete Ordinates". In: *The Astrophysical Journal Supplement Series* 253.2, 49 (Apr. 2021), p. 49.

[23] Y.-F. Jiang. "Multigroup Radiation Magnetohydrodynamics Based on Discrete Ordinates including Compton Scattering". In: *The Astrophysical Journal Supplement Series* 263.1 (Oct. 2022), p. 4.
[24] Rijkhorst, E.-J. et al. “Hybrid characteristics: 3D radiative transfer for parallel adaptive mesh refinement hydrodynamics”. In: *Astronomy & Astrophysics* 452.3 (2006), pp. 907–920.

[25] N. J. Turner and J. M. Stone. “A Module for Radiation Hydrodynamic Calculations with ZEUS-2D Using Flux-limited Diffusion”. In: *The Astrophysical Journal Supplement Series* 135.1 (July 2001), p. 95.

[26] P. Anninos, P. C. Fragile, and S. D. Murray. “Cosmos: A Radiation-Chemo-Hydrodynamics Code for Astrophysical Problems”. In: *The Astrophysical Journal Supplement Series* 147.1 (July 2003), p. 177.

[27] B. Freytag et al. “Simulations of stellar convection with CO5BOLD”. In: *Journal of Computational Physics* 231.3 (Feb. 2012), pp. 919–959.

[28] B. Freytag. “Advances in the hydrodynamics solver of CO5BOLD”. In: *Memorie della Societa Astronomica Italiana Supplementi* 24 (Jan. 2013), p. 26.

[29] B. Freytag. “Boundary conditions in CO5BOLD”. In: *Memorie della Societa Astronomica Italiana* 88 (Jan. 2017), p. 12.

[30] E. Chatzopoulos and K. Weide. “Gray Radiation Hydrodynamics with the FLASH Code for Astrophysical Applications”. In: *The Astrophysical Journal* 876.2 (May 2019), p. 148.

[31] J. M. Stone and M. L. Norman. “ZEUS-2D: A Radiation Magnetohydrodynamics Code for Astrophysical Flows in Two Space Dimensions. I. The Hydrodynamic Algorithms and Tests”. In: *The Astrophysical Journal Supplement Series* 80 (June 1992), p. 753.

[32] J. C. Hayes et al. “Simulating Radiating and Magnetized Flows in Multiple Dimensions with ZEUS-MP”. In: *The Astrophysical Journal Supplement Series* 165.1 (July 2006), p. 188.

[33] M. Gehmeyr and D. Mihalas. “Adaptive grid radiation hydrodynamics with TITAN”. In: *Physica D: Nonlinear Phenomena* 77.1 (1994). Special Issue Originating from the 13th Annual International Conference of the Center for Nonlinear Studies Los Alamos, NM, USA, 17/21 May 1993, pp. 320–341.

[34] González, M., Audit, E., and Huynh, P. “HERACLES: a three-dimensional radiation hydrodynamics code”. In: *Astronomy & Astrophysics* 464.2 (2007), pp. 429–435.
[35] N. M. H. Vaytet et al. "A numerical model for multigroup radiation hydrodynamics". In: *Journal of Quantitative Spectroscopy and Radiative Transfer* 112.8 (2011), pp. 1323-1335.

[36] R. Kannan et al. "arepo-rt: radiation hydrodynamics on a moving mesh". In: *Monthly Notices of the Royal Astronomical Society* 485.1 (Jan. 2019), pp. 117-149.

[37] H. C. Nguyen. "Simulation de modèles hydrodynamiques et de transfert radiatif intervenant dans la description d'écoulements astrophysiques". PhD thesis. Paris 11, 2011.

[38] C. Michaut, H. Nguyen, and L. Di Menza. "Computational radiation hydrodynamics". In: *Astrophysics and Space Science* 336.1 (2011), p. 175181.

[39] C. Michaut et al. "HADES code for numerical simulations of high-mach number astrophysical radiative flows". In: *High Energy Density Physics* 22 (2017), pp. 77-89.

[40] C. Cavet. "Étude numérique de l'instabilité de Vishniac dans les restes de supernovae". Theses. Université Paris-Diderot - Paris VII, Nov. 2010.

[41] J. Minière. "Etude de l'instabilité de Vishniac et régime radiatif des restes de supernova". Theses. Observatoire de Paris, Jan. 2014.

[42] O. Saincir. "Approximation numérique de modèles d'hydrodynamique radiative et applications à la physique stellaire". PhD thesis. Reims, 2019.

[43] A. Gintrand. "Modélisation analytique et numérique de l'évolution des restes de supernova en phase radiative". PhD thesis. Université Paris sciences et lettres, 2019.

[44] C. Busschaert. "Etude théorique, numérique et expérimentale de la dynamique des chocs d'accrétion dans les variables cataclysmiques magnétiques". Theses. Observatoire de Paris, Oct. 2013.

[45] É. Falize. "Similitudes et autosimilarité en physique des hautes densités d'énergie". Theses. Observatoire de Paris, Oct. 2008.

[46] L. Boireau. "Astrophysique de laboratoire: modélisation analytique et numérique du choc radiatif: expériences au moyen de lasers de puissance". PhD thesis. Paris 6, 2005.

[47] W. S. McCulloch and W. Pitts. "A logical calculus of the ideas immanent in nervous activity". In: *The bulletin of mathematical biophysics* 5.4 (Dec. 1943), pp. 115-133.

[48] A. M. Turing. "Mind". In: *Psychology and philosophy* 59 (Oct. 1950).

[49] C. M. Bishop. *Neural Networks for pattern recognition*. Oxford University Press, 1995.
[50] B. M. Wilamowski. "Neural network architectures and learning algorithms". In: *IEEE Industrial Electronics Magazine* 3.4 (2009), pp. 56–63.

[51] M. Alizadeh and Z. H. Zhu. "A comprehensive survey of space robotic manipulators for on-orbit servicing". In: *Frontiers in Robotics and AI* 11 (2024).

[52] K. Thangavel et al. "Artificial Intelligence for Trusted Autonomous Satellite Operations". In: *Progress in Aerospace Sciences* 144 (2024), p. 100960.

[53] K. A. Haj Youssef, Y. Zayegh, and M. Alkhedher. "Autonomous AI-Controlled Mars Rover Robot". In: *2024 Advances in Science and Engineering Technology International Conferences (ASET)*. 2024, pp. 1–7.

[54] A. Ildirimzade and P. Lio. "Advancement in Autonomous Navigation in Space through Artificial Intelligence: A Systematic Review". In: *London Journal* 449 (2025), 449U.

[55] K. Huang et al. "Artificial Intelligence in Astronomical Optical Telescopes: Present Status and Future Perspectives". In: *Universe* 10.5 (2024).

[56] H. U. Nørgaard-Nielsen and H. E. Jørgensen. "Foreground removal from CMB temperature maps using an MLP neural network". In: *Astrophysics and Space Science* 318.3 (Dec. 2008), pp. 195–206.

[57] C. Baccigalupi et al. "Neural networks and the separation of cosmic microwave background and astrophysical signals in sky maps". In: *Monthly Notices of the Royal Astronomical Society* 318.3 (Nov. 2000), pp. 769–780.

[58] F. Farsian, N. Krachmalnicoff, and C. Baccigalupi. "Foreground model recognition through Neural Networks for CMB B-mode observations". In: *Journal of Cosmology and Astroparticle Physics* 2020.07 (July 2020), p. 017.

[59] K. Aylor et al. "Cleaning our own dust: simulating and separating galactic dust foregrounds with neural networks". In: *Monthly Notices of the Royal Astronomical Society* 500.3 (Oct. 2020), pp. 3889–3897.

[60] Casas, J. M. et al. "CENN: A fully convolutional neural network for CMB recovery in realistic microwave sky simulations". In: *Astronomy & Astrophysics* 666 (2022), A89.

[61] K. Prajapati et al. "Contribution of AI and deep learning in revolutionizing gravitational wave detection". In: *Astronomy and Computing* 48 (2024), p. 100856.

[62] A. Mazumder et al. "Synthetic observations with the Square Kilometre Array: Development towards an end-to-end pipeline". In: *Journal of Astrophysics and Astronomy* 44.1 (Mar. 2023), p. 19.
[63] R. Sortino et al. "Radio astronomical images object detection and segmentation: a benchmark on deep learning methods". In: *Experimental Astronomy* 56.1 (Aug. 2023), pp. 293–331.

[64] Y. Li et al. "AI-assisted superresolution cosmological simulations". In: *Proceedings of the National Academy of Sciences* 118.19 (2021), e2022038118.

[65] C. Foucart, A. Charous, and P. F. Lermusiaux. "Deep reinforcement learning for adaptive mesh refinement". In: *Journal of Computational Physics* 491 (2023), p. 112381.

[66] M. J. Smith and J. E. Geach. "Astronomia ex machina: a history, primer and outlook on neural networks in astronomy". In: *Royal Society Open Science* 10.5 (2023), p. 221454.

[67] H. Karimabadi, J. Wilkes, and D. A. Roberts. "The need for adoption of neural HPC (NeuHPC) in space sciences". In: *Frontiers in Astronomy and Space Sciences* 10 (2023).

[68] M. Raissi, P. Perdikaris, and G. Karniadakis. "Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations". In: *Journal of Computational Physics* 378 (2019), pp. 686–707.

[69] I. Lagaris, A. Likas, and D. Fotiadis. "Artificial neural networks for solving ordinary and partial differential equations". In: *IEEE Transactions on Neural Networks* 9.5 (1998), pp. 987–1000.

[70] A. G. Baydin et al. "Automatic differentiation in machine learning: a survey". In: *J. Mach. Learn. Res.* 18.1 (Jan. 2017), pp. 5595–5637.

[71] S. A. E. G. Falle. "Catastrophic cooling in supernova remnants". In: *Monthly Notices of the Royal Astronomical Society* 195.4 (Aug. 1981), pp. 1011–1028.

[72] M. D. Smith. "The stability of radiative shocks". In: *Monthly Notices of the Royal Astronomical Society* 238.1 (May 1989), pp. 235–252.

[73] J. M. Blondin and D. F. Cioffi. "The Growth of Density Perturbations in Radiative Shocks". In: *The Astrophysical Journal* 345 (Oct. 1989), p. 853.

[74] E. Falize et al. "Analytical structure of steady radiative shocks in magnetic cataclysmic variables". In: *Astrophysics and Space Science* 322.1 (Aug. 2009), pp. 71–75.
[75] C. Busschaert et al. "POLAR project: Numerical modeling of the accretion column". In: *High Energy Density Physics* 9.1 (2013), pp. 42–46.

[76] Commerçon, B. et al. "Physical and radiative properties of the first-core accretion shock". In: *Astronomy & Astrophysics* 530 (2011), A13.

[77] N. Baker and R. Kippenhahn. "The Pulsations of Models of δ Cephei Stars." In: *Zeitschrift für Astrophysik* 54 (Jan. 1962), p. 114.

[78] R. F. Christy. "Pulsation Theory". In: *Annual Review of Astronomy and Astrophysics* 4 (Jan. 1966), p. 353.

[79] H. Shibahashi. "Theory of asteroseismology". In: *Journal of Astrophysics and Astronomy* 26.2 (June 2005), pp. 139–160.

[80] Flock, M. et al. "Radiation magnetohydrodynamics in global simulations of protoplanetary discs". In: *Astronomy & Astrophysics* 560 (2013), A43.

[81] J. D. Melon Fuksman et al. "A Two-moment Radiation Hydrodynamics Scheme Applicable to Simulations of Planet Formation in Circumstellar Disks". In: *The Astrophysical Journal* 906.2 (Jan. 2021), p. 78.

[82] V. V. Sobolev. *Course in theoretical astrophysics*. Hayka (Москва), 1985.

[83] G. B. Rybicki and A. P. Lightman. *Radiative processes in astrophysics*. John Wiley & Sons Inc., 1979.

[84] C. D. Levermore and G. C. Pomraning. "A flux-limited diffusion theory". In: *The Astrophysical Journal* 248 (Aug. 1981), pp. 321–334.

[85] G. L. Olson, L. H. Auer, and M. L. Hall. "Diffusion, P1, and other approximate forms of radiation transport". In: *Journal of Quantitative Spectroscopy and Radiative Transfer* 64.6 (2000), pp. 619–634.

[86] C. D. Levermore. "Relating Eddington factors to flux limiters". In: *Journal of Quantitative Spectroscopy and Radiative Transfer* 31.2 (1984), pp. 149–160.

[87] B. Dubroca and J.-L. Feugeas. "Etude théorique et numérique d'une hiérarchie de modèles aux moments pour le transfert radiatif". In: *Comptes Rendus de l'Académie des Sciences - Series I - Mathematics* 329.10 (1999), pp. 915–920.

[88] D. S. Kershaw. *Flux limiting nature's own way – A new method for numerical solution of the transport equation*. Tech. rep. Lawrence Livermore National Lab. (LLNL), Livermore, CA (United States), July 1976.
[89] G. N. Minerbo. "Maximum entropy Eddington factors". In: *Journal of Quantitative Spectroscopy and Radiative Transfer* 20.6 (1978), pp. 541–545.

[90] C. D. Levermore. *Chapman–Enskog approach to flux-limited diffusion theory*. Tech. rep. Lawrence Livermore National Lab (LLNL), Livermore, CA (United States), 1979.

[91] G. C. Pomraning and C. D. Levermore. "A flux limited diffusion theory". In: *J Appl Phys* 248 (1981), p. 321.

[92] G. C. Pomraning. *The equations of radiation hydrodynamics*. Pergamon Press, Oxford, 1973.

[93] R. Turpault. "Construction d'un modèle M1-multigroupe pour les équations du transfert radiatif". fr. In: *Comptes Rendus. Mathématique* 334.4 (2002), pp. 331–336.

[94] R. Turpault. "A consistent multigroup model for radiative transfer and its underlying mean opacities". In: *Journal of Quantitative Spectroscopy and Radiative Transfer* 94.3 (2005), pp. 357–371.

[95] D. Mihalas and B. Weibel-Mihalas. *Foundations of radiation hydrodynamics*. Courier Corporation, 1999.

[96] J. Fort. "Information-theoretical approach to radiative transfer". In: *Physica A: Statistical Mechanics and its Applications* 243.3 (1997), pp. 275–303.

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

[98] G. Radureau, C. Michaut, and A. I. Comport. "AI-based computation method for the Eddington factor in the M1-multigroup model". In: *Phys. Rev. E* 111 (3 Mar. 2025), p. 035301.

[99] J. R. Buchler. "Radiation hydrodynamics in the fluid frame". In: *Journal of Quantitative Spectroscopy and Radiative Transfer* 22.3 (1979), pp. 293–300.

[100] R. B. Lowrie, D. Mihalas, and J. E. Morel. "Comoving-frame radiation transport for nonrelativistic fluid velocities". In: *Journal of Quantitative Spectroscopy and Radiative Transfer* 69.3 (2001), pp. 291–304.

[101] D. Mihalas and L. Auer. "On laboratory-frame radiation hydrodynamics". In: *Journal of Quantitative Spectroscopy and Radiative Transfer* 71.1 (2001), pp. 61–97.