# Higher order methods for Radiative Transfer in Astrophysical simulations: Pn vs M1

Mei Palanque, Pierre Ocvirk, Emmanuel Franck, Pierre Gerhard, Dominique Aubert, Olivier Marchal

## ► To cite this version:

Mei Palanque, Pierre Ocvirk, Emmanuel Franck, Pierre Gerhard, Dominique Aubert, et al. Higher order methods for Radiative Transfer in Astrophysical simulations: Pn vs M1. 2025. hal-05184932

HAL Id: hal-05184932

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

Preprint submitted on 24 Jul 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Higher order methods for Radiative Transfer in Astrophysical simulations: $P_n$ vs $M_1$

M. Palanque¹², P. Ocvirk<fnref n="1" />, E. Franck<fnref n="2" />, P. Gerhard<fnref n="3" />, D. Aubert<fnref n="1" />, and O. Marchal<fnref n="1" />,

Received ?; accepted ?

## ABSTRACT

*Context.* In current cosmological simulations, the radiative transfer modules generally rely on the $M_1$ approximation, which has some glaring flaws related to its fluid-like behaviour, such as spurious pseudo-sources and loss of directionality when radiation fronts from different directions collide. $P_n$, another moment-based model used in other fields of physics, may correct these issues.

*Aims.* We aim at testing out $P_n$ in an astrophysical setting and compare it to $M_1$, in order to see if it can indeed correct $M_1$'s flaws. Also, we want to use $P_n$'s solutions to better pinpoint $M_1$ errors.

*Methods.* We implement a $P_n$ radiation transport method and couple it to a photo-thermo-chemistry module to account for the interaction of ionising radiation with the Hydrogen gas, and benchmark it using tests for radiative transfer models comparison in astrophysics as defined in Iliev et al. (2006).

*Results.* We find that high order $P_n$ (e.g. $P_9$) indeed correct $M_1$'s flaws, while faring as well or even better in some aspects in the tests, in particular when directionality is important or colliding radiation fronts occur. By comparing $P_9$ and $M_1$ radiation fields in an idealised and cosmological test case, we highlight a new, thus far unreported artefact of $M_1$, the 'dark sombrero'. A dark sombrero appears as a spherical photon-deficit shell around the source, at typically 1/3-1/4 of the distance to the next source. The photon density in dark sombreros can be underestimated by a factor up to 2-3. They occur in regions where a source's radiation field connects with that of another source or group of sources. These basic properties (position and amplitude) of the dark sombreros may depend on the sources' relative intensities, positions, spatial resolution, although we have not been able to test this in detail in this study. Moreover, the $M_1$ larger scale photon density also exhibits spurious features, enhancing or reducing photon density in various regions. We use a small reionisation-like test simulation to characterize the relative error in hydrogen neutral fractions between $M_1$ and $P_9$. The relative error is well represented by a gaussian with a dispersion of 0.27 dex in $log_{10}(X_{HI})$. Both aspects are likely related to the photons' collisional behaviour in $M_1$.

**Key words.** Reionisation – Radiative Transfer – Method: numerical – Intergalactic medium

## 1. Introduction

The epoch of reionisation takes place starting 150 My after the big bang, and ends approximately around 1 Gyr after the big bang (redshift 20 to 6) (Barkana & Loeb 2001). During this early period of our universe's history, its entire gas content gets ionised by the recently formed first stars and galaxies. This is an important process of the history of our universe. However, details about how this reionisation took place are still debated, mainly on the matter of the mass range of the galaxies driving it and the nature of ionising sources (stars, compact objects...). This tension originates in the difficult task of deriving the properties of galaxies in such a remote epoch, and the lack of observational constraints for this period. Pending the arrival of such constraints such as the ones that will be brought by the Square Kilometre Array telescope SKA (Dewdney et al. 2009), numerical simulations offer insight into the epoch of reionisation and help us prepare future observational campaigns. Some data hints toward a major role played by the more numerous, small and medium sized galaxies (Lewis et al. 2020; Katz et al. 2019) while some others favour higher mass galaxies (Naidu et al. 2020). There are even proposed scenarios where both the lower and higher ends of

the galaxy mass spectrum drive reionisation at different epochs (Ocvirk et al. 2021).

In all of these simulations, radiative transfer is a central component. With the ionisation of the intergalactic medium (IGM) being the main subject of study of this epoch, the need for a reliable radiative transfer model is central to better describe this process. One of the most commonly used models nowadays is the moment derived model $M_1$ (Levermore 1984; Aubert & Teyssier 2008). However, $M_1$ isn't without flaws, as, for example, it tends to create pseudo-sources at the colliding point of two wave fronts. This stems from the fact that $M_1$'s approximates photons as a fluid, making colliding photon flux add up. It is even more glaringly obvious when observing two perpendicular photon beams cross each other. In a physical setting, the two light beams should not interact with each other, but with the $M_1$ approximation, they merge into one single beam whose direction is the sum of the direction vectors of the two initial beams (Rosdahl et al. 2013). $M_1$'s collisionality could be one of the reasons why the model tends to under estimate the photo-ionisation rate at small scales (Wu et al. 2021) and potentially impact the results of previous simulations.

<footnote n="1">¹ Observatoire Astronomique de Strasbourg, Universite de Strasbourg, CNRS UMR 7550, 11 rue de l'Universite, F-67000 Strasbourg, France</footnote>

<footnote n="2">² Institut de Recherche Mathématique Avancée (IRMA), University of Strasbourg</footnote>

<footnote n="3">³ Direction du Numérique - Pôle CESAR, University of Strasbourg</footnote>
There are many other radiative transfer models in the literature despite the omnipresence of $M_1$ in the field of cosmology simulations (Brunner 2002; Garrett & Hauck 2013). We aim at testing out mainly one of them specifically, the $P_n$ model as described in BRUNNER (2002), Meltz (2015), and observe how it fares compared to $M_1$. The objective is determining if $P_n$ can correct $M_1$'s flaws. We also aim at characterising $M_1$'s artefacts in the observable space (Lyman-$\alpha$ forest transmissions and neutral hydrogen fraction).

## 2. Methodology

In this section, we will describe in more depth how the $P_n$ model works, how it was implemented in our test cases, what kind of chemistry we devised and the code we used to run all of our tests on GPUs.

### 2.1. The $P_n$ model

$P_n$ is a moment derived model of order $n$ of the equation of radiative transfer. Just like $M_1$ or any other moment derived model, they approximate the equation of radiative transfer in the void (Eq. 1) by deriving its moments until reaching the $n$-th order. At this point, a need arises for a closure equation to close the system. This is where the approximations for $M_1$ and $P_n$ diverge.

$$ \frac{1}{c} \partial_t I(r, \Omega, t) + \Omega \cdot \nabla I(r, \Omega, t) = 0, \quad (1) $$

with $\Omega(\theta, \phi)$ the angular unit vector.

$M_1$ is a moment derived model of order 1, which means the closure equation consists in writing the moment of order 2 as a combination of lower order moments. A detailed explanation of the $M_1$ model's closure equation can be found in Aubert & Teyssier (2008). The equation can be projected on a simple $(1, \Omega)$ basis. However, since $P_n$ can technically be derived up to any moment $n$, we need to project our equation on an infinite basis. We'll be using the spherical harmonics basis (Eq. 2, 3) as is done in Meltz (2015), $P$ being the Legendre polynomial.

$$ Y_l^m(\theta, \phi) = \begin{cases} \sqrt{\frac{1(l-m)!}{\pi(l+m)!}} \frac{2l+1}{2} P_l^m(\cos(\phi)) \cos(m\theta), & m > 0 \\ \sqrt{\frac{1}{2\pi}} \frac{2l+1}{2} P_l(\cos(\phi)), & m = 0 \\ -\sqrt{\frac{1(l+m)!}{\pi(l-m)!}} \frac{2l+1}{2} P_l^{-m}(\cos(\phi)) \sin(m\theta), & m > 0 \end{cases} \quad (2) $$

$$ P_l^m(\mu) = (-1)^m \sqrt{((1 - \mu^2)^m)} \frac{d^m P_l}{d\mu^m} \quad (3) $$

We then apply the $P_n$ closure equation, which consists in truncating our projected equations at the order $n$. Let us call $w = \sum_{l=0}^{n} \sum_{m=-l}^{l} w_{l,m}$ our vector solution. In other terms, we approximate our intensity $I$ as:

$$ I(\mathbf{r}, \mathbf{\Omega}) = \sum_{l=0}^{n} \sum_{m=-l}^{l} w_{l,m} Y_{l,m} \quad (4) $$

Article number, page 2

or

$$ \forall l > n, \forall m \in [-l; l], w_{l,m}(\mathbf{r}, \mathbf{r}') = 0 \quad (5) $$

By projecting the whole radiative transfer equation on our basis, we obtain the following equation:

$$ \oint_{4\pi} \Omega \cdot \nabla IY = \oint_{4\pi} \sin(\phi) \cos(\theta) \partial_x IY d\Omega + \oint_{4\pi} \sin(\phi) \sin(\theta) \partial_y IY d\Omega + \oint_{4\pi} \cos(\phi) \partial_z IY d\Omega \quad (6) $$

Recurrence formulas developed in Meltz (2015) show that Eq. 6 can be simplified as:

$$ \frac{1}{c} \partial_t w + J^x \partial_x w + J^y \partial_y w + J^z \partial_z w = 0 \quad (7) $$

Where $I$ is the intensity projected on the spherical harmonics basis, and $J^x, J^y, J^z$ are constant matrices of size $(n+1)^2 \times (n+1)^2$. They are defined in appendix A. This is a very convenient notation since these matrices are constant and only need to be computed once at the beginning of a simulation run. However, they are very sparse despite being quite large, which can take a lot of memory space. The simplicity of the closure equation being a simple truncation means this method is very sensitive to sharp discontinuities, which make it oscillate. On top of that, it isn't based on a physical entropy like $M_1$ is, and thus doesn't guarantee the positivity of the system, which we will delve deeper into these issues in subsection 2.3 and 3.2.

### 2.2. Transport Kernel

We use a classical finite volumes solver for our model, which is described on a three-dimensional grid of fixed size. Each cell can be described by its position in the grid $(i, j, k)$ and its volume $\Delta x * \Delta y * \Delta z$. For all our tests, we'll use cubic cells, hence $\Delta x = \Delta y = \Delta z$. Our time step $\Delta t$ is controlled by a Courant condition (CFL) defined as $CFL = \frac{\Delta t}{h_{min}}$ with $c$ velocity of the light and $h_{min} = \frac{\text{Volume}_{cell}}{\text{Surface}_{cell}} = \frac{\Delta x^3}{6\Delta x^2}$. The CFL must be chosen such as CFL < 1 to ensure stability of our scheme in an explicit case such as ours, determining the value of $\Delta t$. In this paper, the Courant condition will be consistently chosen to be CFL = 0.8, except for a few exceptions detailed at the end of subsection 2.3.

Our code will make use of two separate kernels running successively to perform a complete time step $p$. As such, we will call $p, p+1/2$ and $p+1$ the state of our variables respectively at the beginning of the transport kernel, at the interface between the transport and chemistry kernels, and at the beginning of the next time step $p+1$. This notation will be used throughout this subsection and the next one dedicated to the chemistry kernel.

To compute our explicit transport time step using the $P_n$ model, we use a Rusanov scheme (Rusanov 1962) as advised in Meltz (2015) Sahmim (2005) for its great stability. For comparison purposes, and as opposed to Aubert & Teyssier (2008) we're referencing, our own implementation of $M_1$ will also use that scheme.
Using Eq. 7, and defining S the source function, we can write:

Let us call $w_l$ the factor corresponding to moment of order 1 in our solution vector $w$. We can derive the moments of the radiative transfer equations until the n-th order as follows, with S the source function:

$$ \frac{1}{c} \partial_t w + J^x \partial_x w + J^y \partial_y w + J^z \partial_z w = S \quad (8) $$

This system can be discretised as follows, for a time step p and positions $(i, j, k)$:

$$ \forall l \in [1, n] $$

$$ \frac{w_l^{p+1/2, i, j, k} - w_l^{p, i, j, k}}{\Delta t} + c \left[ \frac{w_{l+1}^{p, i+1/2, j, k} - w_{l+1}^{p, i-1/2, j, k}}{\Delta x} + \frac{w_{l+1}^{p, i, j+1/2, k} - w_{l+1}^{p, i, j-1/2, k}}{\Delta y} + \frac{w_{l+1}^{p, i, j, k+1/2} - w_{l+1}^{p, i, j, k-1/2}}{\Delta z} \right] = S \quad (9) $$

With $i + 1/2$ the flux at the interface between the cells $i$ and $i + 1$ of our grid, in the $x$ direction. This flux can be computed with the Rusanov scheme as follows:

$$ w^{p,i+1/2,j,k} = \frac{1}{2} J^x(w^{p,i+1,j,k} + w^{p,i,j,k}) - \frac{B}{2}(w^{p,i+1,j,k} - w^{p,i,j,k}) \quad (10) $$

In the same way, we can compute the two other directions:

$$ w^{p,i,j+1/2,k} = \frac{1}{2} J^y(w^{p,i,j+1,k} + w^{p,i,j,k}) - \frac{B}{2}(w^{p,i,j+1,k} - w^{p,i,j,k}) \quad (11) $$

$$ w^{p,i,j,k+1/2} = \frac{1}{2} J^z(w^{p,i,j,k+1} + w^{p,i,j,k}) - \frac{B}{2}(w^{p,i,j,k+1} - w^{p,i,j,k}) \quad (12) $$

With B a majoration of the spectral radius of the matrices $J$. Said spectral radius being $]-1; 1[$ (Meltz 2015), we chose $B = 1$ for our implementation. With this scheme in place, we can compute multiple adimensional tests of section 3. However, to accurately test the abilities of $P_n$ in an astrophysical setting, there's a need for a dimensional chemistry kernel. This kernel runs after the transport kernel, returns the temperature and ionised fraction of the simulation, as well as some feedback photon density value $w$ due the absorption of photons by neutral hydrogen atoms and recombination.

## 2.3. Chemistry Kernel

Our chemistry kernel works separately from the transport kernel devised in 2.2. During a time step $p$, it takes in the photon density $w_0^{p+1/2}$ from our transport kernel, defined as our coefficient of order 0, and outputs the updated temperature $T^p(r)$, ionised fraction $x^p(r)$, and photon density variation $dw_0^p(r)$. This last parameter is then used to update all moments in the transport kernel to take into account the photon recombination in the energy budget. As a reminder, the only chemical compound we take into account is hydrogen, and its density $n_H$ stays constant in time as a result of the absence of a hydrodynamical solver in our tests.

With $n_{HI}$ and $n_{HII}$ respectively density of neutral and ionised hydrogen, $n_e$ the density of free electrons equal to $n_{HII}$, and $x$

the fraction of ionised hydrogen, the hydrogen density $n_H$ can be written as follows:

$$ n_H = n_{HI} + n_{HII} = (1-x)n_H + xn_H \quad (13) $$

With that in mind, we aim at solving the following equations to obtain the updated $w_0, w_k, x$ and $T$:

$$ \frac{\partial w_0}{\partial t} = -(1-x)n_H c \sigma w_0 + x^2 n_H^2 (\alpha_A - \alpha_B) \quad (14) $$

$$ \frac{\partial w_k}{\partial t} = -(1-x)n_H c \sigma w_k, \forall k \in [1, n] \quad (15) $$

$$ \frac{\partial x}{\partial t} = (1-x)c\sigma w_0 - x^2 n_H^2 \alpha_A - x(1-x)n_H\beta \quad (16) $$

$$ \frac{\partial E}{\partial t} = \mathcal{H} - \mathcal{L} \quad (17) $$

Where $\alpha_A(T)$, $\alpha_B(T)$ and $\beta(T)$ correspond respectively to the case A and case B recombination rates, and the HI collisional ionisation coefficient of HII as defined in Hui & Gnedin (1997), $\sigma$ is the effective HI cross section at $e_{HI}$, $E(T, t)$ the thermal energy of the gas and $\mathcal{H}(x, w_0, n_H)$ and $\mathcal{L}(T, n_H)$ are the photoionisation heating rate and cooling rate. In all of our test cases of section 4, we consider $10^5$ K black body sources to comply with literature test cases (Iliev et al. 2006) implying a cross-section $\sigma = 1.63 \times 10^{-22} \text{ m}^2$ at $e_{HI} = 29.61$ (Aubert & Teyssier 2008), except for test 1 where we consider $3.10^4$ K black body sources and $\sigma = 6.3 \times 10^{-22} \text{ m}^2$ at $e_{HI} = 13.6$ eV (Osterbrock 1974).

We use an implicit scheme to solve the photon density variation and ionised fraction, since these processes happen at a small time scale as described in Aubert & Teyssier (2008). However, for sake of simplicity, our temperature will use an explicit scheme.

As we said previously in subsection 2.1, and as we'll delve deeper into in subsection 3.2, $P_n$, in essence, does not guarantee that the photon density $w_0^{p+1/2}$ at the end of a transport step is positive. This does not concern us while in the transport kernel, but it cannot stay this way when entering the chemistry kernel. Our way around this possible non positivity of the density is to approximate it as zero in negative regions, as shown in Eq. 18. This might seem like a strong choice, but as shown in Appendix B, the negative photon density output is mostly negligible as long as $P_n$'s order is sufficiently high in spherical cases.

$$ \bar{w}_0^{p+1/2} = \max(w_0^{p+1/2}, 0) \quad (18) $$

In the chemistry kernel, we'll only ever use this truncated version $\bar{w}_0$ of the photon density.

To follow the rate of reionisation, the neutral fraction of gas is a key observable that can be derived from the 21cm emission of neutral hydrogen. This parameter can be defined as $x_{HI} = 1-x$ where $x$ is the fraction of ionised hydrogen in our simulation.

We can write the coupled equation of the first moment in 1D and the ionising equation as follows:

$$ \frac{dw_0}{dt} + \frac{dw_1}{dr} = \dot{w}_0 + \dot{w}_0^{\text{rec}} - n_H \sigma c w_0 (1-x) \quad (19) $$

Where $\dot{w}_0$ and $\dot{w}_0^{\text{rec}}$ correspond respectively to ionising sources and ionising photons from recombination. Assuming
that $\alpha_A n_H^2 x^2 = \dot{w}_0^{rec} + \alpha_B n_H^2 x^2$, and using Eq. 16, we can rewrite Eq. 19 as:

$$ \frac{dw_0}{dt} + \frac{dw_1}{dr} = \dot{w}_0 - \alpha_B n_H^2 x^2 + \beta n_H^2 x(1-x) - n_H \frac{dx}{dt} \quad (20) $$

Now, we can approximate $\frac{dw_0}{dt} = \frac{w_0^{p+1} - w_0}{\Delta t}$. We'll also write $x^p = x^{p+1/2} = x$ and $x^{p+1} = X$ for sake of comprehension. Doing this, and knowing that $w_0^{p+1/2} = w_0^p + \dot{w}_0 - \frac{dw_1}{dr}$ is the solution to the transport equation in the void already solved in section 2.2, we can rewrite Eq. 20 as:

$$ w_0^{p+1} = w_0^{p+1/2} + \beta n_H^2 (1-X) X \Delta t - \alpha_B n_H^2 X^2 \Delta t - n_H(X-x) \quad (21) $$

Replacing $w_0^{p+1}$ in Eq. 19 by this new definition, and using the truncated version of the photon density $\bar{w}_0^{p+1/2}$ instead of the real value $w_0^{p+1/2}$ to ensure positivity, we obtain a third degree polynomial in $X$ defined as follows:

$$ mX^3 + nX^2 + pX + q = 0 \quad (22a) $$

$$ m = (\alpha_B + \beta)n_H^2 \Delta t \quad (22b) $$

$$ n = n_H - \frac{(\alpha_B + \beta)n_H}{\sigma c} - \alpha_B n_H^2 \Delta t - 2\beta n_H^2 \Delta t \quad (22c) $$

$$ p = \beta n_H^2 \Delta t - n_H(1+x) - \bar{w}_0 - \frac{1}{\sigma c \Delta t} + \frac{\beta n_H}{\sigma c} \quad (22d) $$

$$ q = \bar{w}_0 + n_H x + \frac{x}{\sigma c \Delta t} \quad (22e) $$

Solving this system by finding the only real root between 0 and 1 will give the updated value of the ionised fraction $x^{p+1}$ in the cell. Given our knowledge of this updated ionised fraction and Eq. 21, we define the updated photon density as follows:

$$ w_0^{p+1} = w_0^{p+1/2} + dw_0^{p+1/2} \quad (23) $$

$$ \text{with } dw_0^{p+1/2} = \beta n_H^2 (1-X) X \Delta t - \alpha_B n_H^2 X^2 \Delta t - n_H(X-x) $$

This way, we add the photon density variation $dw_0^{p+1/2}$ derived from the truncated photon density $\bar{w}_0^{p+1/2}$ to the real photon density $w_0^{p+1/2}$. As such, even if $dw_0^{p+1/2}$ is strictly positive, $w_0^{p+1}$ can be negative if $w_0^{p+1/2}$ already was.

In the same way as done in Eq. 19, we can write all coefficients of order superior to 0 as follows:

$$ \frac{dw_l}{\Delta t} + \frac{dw_{l+1}}{dr} = -n_H \sigma w_l (1-x) \quad (24) $$

As such, using an implicit scheme and the updated $X$ previously computed, the updated value of all coefficients of order superior to 0 is written:

$$ w_l^{p+1} = \frac{w_l^{p+1/2}}{1 + n_H \sigma c \Delta t (1 - X)} \quad (25) $$

With $w_l^{p+1/2} = w_l^p - \frac{dw_{l+1}^p}{dr} \Delta t$ value of the coefficient at the end of a transport step. This way, all of our coefficients are updated with photon variation due to absorption.

The last parameter to be computed is the temperature of the hydrogen gas. Our code follows its variations through heating

and cooling processes. The heating rate $\mathcal{H}$ mainly involves photoionisation and is given by:

$$ \mathcal{H} = n_H(1-X)\bar{w}_0^{p+1/2} \sigma \Delta e \quad (26) $$

Where $\Delta e = e_{HI} - e_{HII}$ where delta e is the energy leftover in the unbound electron and proton system after a photo-ionisation, available as thermal energy. In our case, this difference is equal to $\Delta e = 29.61 - 13.6$ eV. For stability reasons, we chose to use the updated value of the ionised fraction $X$, but the non-updated value of the photon density $w_0^{p+1/2}$. The cooling rate $\mathcal{L}$ is the result of collisional cooling due to case A and B recombination, collisional ionisation and excitation, and bremsstrahlung effect. We use fits from Hui & Gnedin (1997) and Maselli et al. (2003) to compute these processes.

Our internal gas energy $E = \frac{3}{2}n_{tot}k_B T$ is linked to these two rates by Eq. 17. Knowing that $n_{tot} = n_e + n_H = xn_H + n_H = n_H(1+x)$, it can rewritten as:

$$ \frac{3}{2} k_B n_H \left[ T \frac{\partial(x+1)}{\partial t} + (x+1) \frac{\partial T}{\partial t} \right] = \mathcal{H} - \mathcal{L} \quad (27) $$

Which, when discretising, gives the following explicit equation for the evolution of the temperature, remembering that $T^p = T^{p+1/2}$, $x^p = x^{p+1/2} = x$ and $x^{p+1} = X$:

$$ \frac{T^{p+1} - T^{p+\frac{1}{2}}}{\Delta t} = \frac{1}{X+1} \times \left[ \frac{2[\mathcal{H}(x, X, \bar{w}_0^{p+1/2}) - \mathcal{L}(X, T^{p+\frac{1}{2}})]}{3k_B n_H} - \frac{(X-x)T^{p+1}}{\Delta t} \right] \quad (28) $$

Since computing the updated temperature $T^{p+1}$ requires prior knowledge of the updated ionised fraction $X$, it is computed in last during a chemistry time step.

We now have updated all chemistry variables for the time step $p+1$. Updated values of $w^{p+1}$ are sent back to the transport kernel as new initial conditions for the next time step. We note that we didn't implement any subcycling for our chemistry kernel in this paper despite the possibility of $\Delta t$ being too large to accurately describe the short time scale behaviour of our chemistry. In cases where the time step would be too great to be handled by the kernel, we opted for a diminished CFL, which increases our computational cost. However, subcycling our chemistry kernel would ensure a more optimised implementation of the code and should be used for further use of the model.

In conclusion, our scheme is made of Eq.9 for the transport, and Eq.22,23,25 and 28.

## 2.4. RKMS

Moment based models can quickly become quite computationally expensive as the order of the model increases. Indeed, each moment requires an increasing number of coefficients. The moment of order 0, representing the isotropic component of our radiation, only requires one, but the first order already needs three more (the three directional flux components). For a given order $n$ of a moment based model, the number of coefficients to be computed at each time step, in each cell, is $(n+1)^2$.
<figure><img src="image_2.png" /><figcaption>Fig. 1: Execution time comparison of RKMS compared to our Julia implementation on CPU and GPU depending on the box size</figcaption></figure>

This means that, while $M_1$ only requires the computation of four coefficients at each time step, $P_9$ will need a hundred. This quickly became a problem for our python-based and then julia-based CPU implementations, which struggled to compute our models in a reasonable time frame.

Reduced Kinetic Model Solver (RKMS) is an openCL based, python wrapped, mono-GPU implementation of $P_n$ developed by Pierre Gerhard during his PHD thesis https://github.com/p-gerhard/rkms, (Gerhard 2020). Its fully optimised GPU core makes it incredibly faster than our previous implementations, at the cost of the flexibility. Extensive work on RKMS was necessary to implement the chemistry kernel and the dimensionalisation required for the test cases in section 4, but allowed us to run a lot of tests far more quickly than we expected, as shown in Fig. 1. However, another limit of RKMS is that it is not parallelised on several GPUs, and thus can still take a long time to compute higher orders of $P_n$. As a result, the biggest box size used in this paper will be a 128³ box with $P_9$ implementation.

For our GPU runs, whose results will be shown throughout this paper, we relied heavily on the services of the High Performance Computing Centre of Strasbourg (HPC).

### 3. Results: Radiation only test cases

We start by highlighting the ways in which $P_n$ solves $M_1$'s issues in a series of simple tests. All test cases in this section are purely adimensional radiative transfer models with no chemistry coupling nor hydrodynamics and gravity, ran on a 64³ grid. Their aim is to showcase how $P_n$ compares to $M_1$ in the specific cases where it tends to have issues.

#### 3.1. Continuous Isotropic Sources

Our first and most straightforward test case serves to illustrate the main improvement $P_n$ brings to the table compared to $M_1$. Indeed, as talked about previously, $M_1$ suffers from its fluid-like

approximation of photons resulting in a collisional behaviour. As such, two wave fronts colliding end up adding up instead of just crossing each other as would be physically expected. This tends to create pseudo-sources and ejecta of energy that have no physical reality. This mainly appears when the photon fluxes coming from several sources interact with each other.

To highlight this flaw in $M_1$ and see if $P_n$ can correct it, we put two continuous isotropic sources next to each other in a 64³ box. Both sources are adimensional of intensity equal to 1 per time step. We want to observe the consequence of their wave front interacting in the continuous regime. Knowing the issues $M_1$ usually have in this kind of setting, we expect the apparition of a pseudo-source and beams of energy being emitted from it, orthogonally to the direction of the flux. As for $P_n$, we expect none of this behaviour, with simply two wave fronts crossing each other without interaction.

Results of this test for $M_1$ and $P_9$ are showcased in Fig. 2. Immediately, the pseudo source created by the coalescence of two wave fronts becomes evident in $M_1$, which emits vertically and changes the directionality of the photon flux. It is also visible in the profile made along the line formed by the two sources, where a bump appears in the photon distribution between the sources, highlighting their interaction. This behaviour is completely absent in $P_9$, where the bump is absent in the profile, and no coalescence is apparent in the slice shown in the same figure, on the right. The directionality of the energy is conserved, and the continuous regime is what we should physically expect from two isotropic sources placed next to each other. As such, we can confidently say that $P_n$ corrects this issue of $M_1$ in the isotropic case.

#### 3.2. Continuous non-isotropic sources

The issues of $M_1$ can be highlighted even further in a more extreme test case that shows how detrimental it can be to the physics of a simulation. Indeed, the directionality change highlighted in the previous case can be shown even more strongly. In the next test case, we will put two adimensional non-isotropic continuous sources in a 64³ box. Those directional
<figure><img src="image_3.png" /><figcaption>Fig. 2: Comparison of colliding fronts of two isotropic continuous sources in M₁ and P₉</figcaption></figure>

sources each emit a beam of light that cross at the centre of the box. Physically, we expect the two light beams should just cross each other without interacting, since light is non collisional. However, considering M₁'s peculiar properties, we expect the two beams to coalesce and their resultant directional vector being the sum of the two initial direction vectors.

Results for this test case are shown in Fig. 3 and highlight even more blatantly the flaw of M₁. Here, the two beams merge and the resulting new beam changes direction, in a completely non-physical way. However, we can see that P₉ manages to maintain the crossing behaviour that should be modelled here, and corrects once again the issues of M₁. However, we can also see a lot of oscillations around the sources in P₉ which do not look very physical either. Indeed, as mentioned previously, Pₙ has a hard time with discontinuities, and the non-isotropy of the sources is one such discontinuity. It creates ripple effects in all directions around the sources, which are very low intensity compared to the beam. We also mention the case of temporal discontinuities in Appendix B.

## 4. Results: Coupled radiation and chemistry - Cosmological Radiative Benchmark Tests

In this section, we show how Pₙ fits the requirements for a Cosmological radiative transfer model by making it undergo four separate test cases described in Iliev et al. (2006). These stan-

dardised tests allow for an easy comparison with other radiative transfer models in the literature, but we will keep making a comparison with the M₁ model for sake of comprehension. More in depth analysis of M₁'s results for the same tests can be found in Aubert & Teyssier (2008). All test cases in this section are dimensional, i.e. using physical box sizes and emissivities, with coupled chemistry as described in 2.3 and no hydrodynamics or gravity.

### 4.1. Isothermal Strömgren Sphere

A single ionising source can only ionise a given radius of gas around it before equilibrium is reached and we enter the continuous regime. The sphere of such ionised gas is referred to as the Strömgren Sphere, which radius is given by the following equation:

$$R_s = \frac{3S}{4\pi\alpha_B(T)n_H^2}^{1/3} \quad (29)$$

With S the intensity of the ionising source in photon.s⁻¹, nH the density of the gas in m⁻³, T the gas temperature in K and αB the case B recombination rate of hydrogen. The characteristic evolution time of the front can be derived as tₚ = 1/αB(T)nH. As such, with a unique non variable source and fixed temperature
<figure><img src="image_4.png" /><figcaption>Fig. 3: Beam crossing comparison between M₁ (Left) and P₉ (Right). We can observe the difference between how the beams cross each other in P₉ and average out in M₁</figcaption></figure>

<table><caption>Table 1: Initial Condition Parameters for the benchmarked tests</caption><thead><tr><td>Test Name</td><td>Duration (Myr)</td><td>Cross section (m<sup>-2</sup>)</td><td>T(t<sub>0</sub>) (K)</td><td>n<sub>h</sub> (m<sup>-3</sup>)</td><td>x<sub>HII</sub>(t<sub>0</sub>)</td><td>Source</td></tr></thead><tbody><tr><td>Isothermal Strömgren Sphere</td><td>500</td><td>6.3 &times; 10<sup>-22</sup></td><td>1 &times; 10<sup>4</sup></td><td>1 &times; 10<sup>3</sup></td><td>1.2 &times; 10<sup>-3</sup></td><td>5 &times; 10<sup>48</sup> photons/s</td></tr><tr><td>Adiabatic Strömgren Sphere</td><td>100</td><td>1.63 &times; 10<sup>-22</sup></td><td>100</td><td>1 &times; 10<sup>3</sup></td><td>1.2 &times; 10<sup>-3</sup></td><td>5 &times; 10<sup>48</sup> photons/s</td></tr><tr><td>Shadowing a dense clump</td><td>3</td><td>1.63 &times; 10<sup>-22</sup></td><td>8 &times; 10<sup>4</sup></td><td>2 &times; 10<sup>2</sup></td><td>1.2 &times; 10<sup>-3</sup></td><td>10<sup>10</sup> photon/m<sup>2</sup></td></tr><tr><td>Cosmological map</td><td>4</td><td>1.63 &times; 10<sup>-22</sup></td><td>100</td><td>Cube</td><td>1.2 &times; 10<sup>-3</sup></td><td>List</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Fig. 4: Front radius evolution of P₅ for several reductions of c</figcaption></figure>

and hydrogen density, the theoretical radius of this sphere can easily be computed, and compared to the result of our model. With that in mind, we devise a test similar to the one given in Iliev et al. (2006) consisting of a single continuous isotropic source in the centre of a simulation box, emitting $S = 5 \times 10^{48}$ photons.s⁻¹. The surrounding hydrogen as a fixed density $n_H = 1 \times 10^3$ m⁻³ and temperature $T = 1 \times 10^4$ K. The initial value of ionised fraction is $1.2 \times 10^{-3}$. The grid uses a 65³ cube representing a 13.2 kpc side with reflexive boundary conditions.

Fig. 4 shows the evolution of the position of the ionising front over the Strömgren radius over the length of the simulation, for several values of reduction for the speed of light c, along with a theoretical evolution of the front derived as $[1 - \exp(-t/t_r)]^{1/3}$. Here, the simulation runs for a little bit over four times the characteristic evolution time, for about 500 Myrs. Fig. 5 shows side by side, profiles of the ionised and neutral fractions at 35 Myr on the left, and 500 Myr on the right. Both P₅ and M₁ are compared here, with their relative difference being showcased in the underplot.

With several values of reduced c, we can observe that P_n converges toward the theoretical value of the strömgren sphere in Fig. 4. However, we also can see that the variation in value of c also changes the celerity of the front line, and could impact results in further tests, despite allowing for a far shorter computation time. All in all though, P_n passes this test easily, and matches the results of M₁ very closely, as shown in Fig. 5, as the difference between the two models never oversteps 1%.

## 4.2. Non-isothermal Strömgren Sphere

A similar experiment can be made, but this time with temperature variation being taken into account. The conditions of the test are similar to the one previously described, but with an initial condition in temperature of 100 K and the source being a 10⁵ K black body. This time, the computation of the theoretical strömgren sphere radius is less straightforward, but we use a comparison with M₁ as done previously, as well as with the previous test, to determine how well P_n is doing.
<figure><img src="image_6.png" /><figcaption>Fig. 5: Ionised and Neutral hydrogen profile comparison of M₁ and P₅ around 35 Myr (Left) and 500 Myr (Right) in the isothermal case</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6: Ionised and Neutral hydrogen profile comparison of M₁ and P₅ around 35 Myr (Left) and 100 Myr (Right) with coupled temperature evolution</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7: Comparison of the temperature profiles in the strömgren sphere at three times, 10, 35 and 100 Myr</figcaption></figure>

In Fig. 6, are shown side by side the profiles of the ionised and neutral fractions at 35 Myr on the left, and 100 Myr on the right. As done previously, we compare both models, with their relative difference being showcased in the underplot. Fig.

7 shows the same comparison but on the temperature profile around our source, and at three different time steps instead of two: 10, 35 and 100 Myr.

Once again, P_n matches M_1 very well, both in ionised fraction and in temperature, with maximum differences of around 2% for the former and 5% in the latter. Surprisingly, we observe that this variation in ionised fraction has been divided by a factor two compared to the same test without temperature variation at the same time step. The temperature profile also fits what is physically expected of this test, with a sudden drop when the ionising front is reached.

### 4.3. Shadowing a dense clump

Previous test cases took place in a context in which P_n is known to fare pretty well, i.e. with continuous, isotropic sources that tend to limit the oscillations of the model. We've confirmed that, in such cases, P_n fares as well as M_1. However, the third test in Iliev et al. (2006) makes use of non-isotropic sources that may cause more issue for the model.
<figure><img src="image_9.png" /><figcaption>Fig. 8: Dense clump of hydrogen in the path of a flux of ionising radiations using $M_1$ (top) and $P_9$ (bottom) at 3 Myr</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9: Example slice of hydrogen density map from Iliev et al. (2006) used in this test</figcaption></figure>

The shadowing by a dense clump test aims at mimicking how resistant to ionisation a clump of dense, cold gas would be using our radiative transfer model. It consists in a cold, high-density sphere of hydrogen placed on the path of a constant directional flux of photons. The sphere - hereby referred to as "clump" - should resist full ionisation, as it is too dense for the radiation to properly penetrate it, instead creating a slim shell of ionised gas on its exposed surface. On top of that, as the clump blocks the path of the ionising flux, it should create a "shadow" of non-ionised, non-heated gas in its wake.

For our experiment, we define the background hydrogen density as $0.2 \text{ h}^{-1} \text{ m}^{-3}$ comoving, which correspond to the mean background density of the universe. In the same way, the mean density of a dark matter clump is two hundred times greater at $40 \text{ h}^{-1} \text{ m}^{-3}$. Considering a redshift of $z = 9$ for our test, this turns

into a high density hydrogen sphere of $4 \times 10^4 \text{ m}^{-3}$ plunged in a background gas density of $2e2 \text{ m}^{-3}$. We use a 6.6 kpc box with a clump of radius 0.8 kpc that runs for 3 Myr. The initial temperature is homogeneous at 8000 K, unless in the clump where it is set at 80 K. The ionising flux is emitted along the z direction by a plan of source cells, at a value of $1 \times 10^{10} \text{ s}^{-1} \text{m}^{-2}$. However, as shown in Fig. 3, non-isotrope sources in $P_n$ have modes and oscillations, which means that some of the energy, even if negligible, isn't aligned with the main direction of our beams. Worse, some of this photon density goes in the opposite direction, and can create further oscillations if met with the discontinuity of the absorbent boundary conditions. To ensure that such oscillations don't impact the result of the experiment, we create buffers between the cells used in our model and the border of our box, which dampen the intensity of our oscillations and their impact on the test itself. As such, we will use nested boxes, with, at the centre, our $64^3$, 6.6 kpc box, inside a $128^3$, 13.2 kpc one that will serve as a buffer. Indeed, the boundary conditions can interact with the ripples created by non-isotropic $P_n$ sources and create oscillations which prevents a fair comparison between the two models. This bigger box will serve to minimise the impact of boundary conditions on the simulation box. To ensure a fair comparison, $M_1$ will use the same box as $P_n$. In all plots of this test case, we will only showcase the main simulation box and omit the buffers. The centre of the clump is placed at coordinates [0.5,0.5,0.25] in the main simulation box.

We can observe the results of $P_9$ compare well to $M_1$ in Fig. 8, where both models show a similar ionisation of the clump. On top of that, we can see that $P_9$ creates a stronger shadow behind said clump than $M_1$, which, in this case, is closer to what should be physically expected. Indeed, on top of the already present numerical diffusion caused by the Rusanov scheme used for both models, $M_1$ tends to be an already quite diffusive model, which leads to the flux of photons encroaching on the
<figure><img src="image_11.png" /><figcaption>Fig. 10: Maps of the photon density at 0.4 Myr in $M_1$ (Left) and $P_9$ (Right) with a limited dynamic</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11: Maps of the neutral fraction of hydrogen at 0.4 Myr in $M_1$ (Left) and $P_9$ (Right)</figcaption></figure>

shadow. This was already shown in Aubert & Teyssier (2008), where comparison with a ray-tracing model highlighted this behaviour of the model. In this case, we can confidently say that $P_9$ performs better compared to $M_1$, which is a good sign.

One can however observe the modes in the photon density of $P_9$. This is the consequence of the non-directional components of $P_9$ and the oscillations that result from it. However, this difference compared to $M_1$ has little to no impact on the temperature and neutral fraction maps, and is specific to this situation adverse to $P_n$ that wouldn't be found in most simulations.

4.4. Cosmological density map

The fourth and final test is also the one who comes closest to the future applications this model could have in cosmological simulations of the epoch of reionisation. It consists in testing the behaviour of our radiative transfer model on a fixed 128³ cells grid of hydrogen density provided by the authors of Iliev et al. (2006), along with sixteen isotropic continuous sources of varying intensities placed in that cosmological field. This situation mimics a real simulation cube (without hydro or gravity) at a smaller scale to get a better idea of how the model will react in situations more practical than the previously mentioned test cases. We expect to observe strong differences between $M_1$ and $P_9$ due to the interaction between several sources and the pseudo-sources we've shown in section 3, which
<figure><img src="image_13.png" /><figcaption>Fig. 12: Position of the 50% ionisation front at 0.4 Myr</figcaption></figure>

may cause a delay in the propagation front in M₁ compared to Pₙ. Indeed, since these unwanted interactions tend to change the directionality of part of the photon flux, the energy brought to the main ionising front might be diminished compared to Pₙ, where this problem doesn't arise. We also aim at observing the differences between several orders of Pₙ, especially P₃, P₅, P₇ and P₉, to observe at which point the solution seems to converge. As for Pₙ potential negativity, we expect it to have a minor impact of the test, as this test doesn't hold any discontinuities that might cause the model to strongly oscillate, apart from potential boundary condition effects. Finally, we'll also push the test further than was done in the initial paper to reach the optically thin regime where M₁ is known to struggle, to see how much better Pₙ fares compared to it.

The simulation box represents a 0.5/h Mpc sized cube comoving, which, with $h = 0.7$ in this specific case, and considering a redshift $z = 9$, translates to a 71.4 kpc sized box. The initial temperature is set to be 100 K, and the initial ionised fraction $1.2 \times 10^{-3}$. The 16 sources range from $0.64 \times 10^{52}$ photons.s$^{-1}$ to $7.97 \times 10^{52}$ photons.s$^{-1}$. The full test is described at https://astronomy.sussex.ac.uk/~iti20/RT Comparison Project/tests.html. The binary file containing the hydrogen density can also be found on the same website. For our analysis, we'll observe specifically one slice of this box, at coordinate $z = 64$, of which we show the hydrogen density map in Fig. 9. We first observe the results at $t = 0.4$ Myr as was done in the original papers, before running the simulation until $t = 4.0$ Myr to reach the optically thin regime. It is to note that, at this time frame, M₁ appears to have some issues ionising some part of the simulation box due to boundary conditions. The hydrogen density map being periodic, but not the boundary conditions of the box, we shift the cells to avoid the appearance of that artefact in the probability distribution functions showed later.

At 0.4 Myr, we can already spot some strong differences between the two models. Qualitatively, it appears quite obvious in Fig. 10 that photons in P₉ seem to propagate faster than M₁, but also seem more evenly distributed in the already ionised portion of the slice we observe. This translates in the neutral fraction too as shown in Fig. 11, where P₉ seems to have more chunks of ionised regions than M₁ at the same time. We can quantify this difference, by observing the position of the ionising front at 0.4 Myr, defined as the point where the gas is 50% ionised, or

$x = 0.5$. This is shown in Fig. 12, where it appears even more strikingly that, not only P₉, but also P₇, P₅, and even P₃ are ionising faster than M₁. This may be the result of the interaction between sources in M₁ and the artifacts they create, resulting in less energy contributing to the ionising of neutral hydrogen. We can also point out that all instances of Pₙ seem to converge toward a single solution as the order increases, the highest tested order being P₉. It looks like, based on this metric, P₇ could be a potentially less expensive alternative to P₉ with similar results. To better compare the homogeneity of P₉ and M₁ in ionised regions, we can look at the probability distribution function (PDF) of the photon density and the neutral fraction, which is done if Fig. 13. We observe here a significant difference between M₁ and P₉, with a higher number of highly neutral cells and a higher number of low photon density cells in M₁. We can also point out that P₃ diverges quite a lot from both P₅, P₇ and P₉ both in photon density distribution and neutral fraction distribution, which tends to show that it is too low of an order to use, as it strays too much from the converged solution. Yet, in Fig. 10, it appears already obvious that, in the ionised regions, a large discrepancy is already present between M₁ and P₉, with huge energy ejectas in M₁, mainly in the top right corner, which are absent in the other model, and incompatible with the physically expected results of such an experiment. This effect, that can already be seen in the PDF, is highlighted further in Fig. 14, where we plot the ratio of neutral fraction of Pₙ over M₁ ($x_{HI}^{P_9}/x_{HI}^{M_1}$) along with a histogram of this ratio in the whole simulation box. Three regimes appear:

- a very piqued regime with $\sigma = 0.01$ corresponding to the regions where the gas is neutral in both M₁ and P₉. Those are the areas not yet reached by radiations in both models

- a very broad regime with $\sigma = 3.47$ corresponding to regions where the gas is neutral in either M₁ or P₉. Those areas match the propagation discrepancy between the two models

- a medium regime with $\sigma = 0.26$ corresponding to regions where the gas is ionised in both models. Those areas have reached the optically thin regime

The third regime is especially relevant as it quantifies the theoretical error on $x_{HI}$ due to the radiative transfer model, as would be observed in the Lyman-$\alpha$ forest. In conclusion, we can say that, at 0.4 Myr, Pₙ ionises the medium faster than M₁, and that the inhomogeneities appearing in the ionised region in M₁ are already significant despite being outmatched by the difference in velocity between the two models.

We then observe the simulation at 4.0 Myr, in the optically thin regime. Here, there are no ionising front to observe any more as the whole box should be ionised already. Yet, when we observe neutral fraction maps in Fig. 15, we can very quickly see that M₁ is less ionised than P₉, with bigger patches of more neutral hydrogen in the higher density regions. These patches aren't neutral per se, as their neutral fraction is quite low ($10^{-3}$), but are far less ionised than P₉ even though the simulation has converged by now, as shown by the mean neutral fraction evolution in the left plot of Fig. 16, since both P₉ and M₁'s neutral fractions have reached a plateau or a near plateau. The apparent difference in ionisation can be found in this plot too, as the final mean neutral fraction in M₁ is slightly higher than all the mean values for Pₙ models, who seem to converge toward a value of $2 \times 10^{-6}$. M₁ is also a bit late compared to most Pₙ models, converging 0.2 Myr later, which could also have an impact on simulations as the timing of reionisation is an important open question. This contrast translates in the PDF of the neutral fraction distribution, right plot of Fig. 17. Here, it is
<figure><img src="image_14.png" /><figcaption>Fig. 13: **Left:** Probability distribution function of photon density at 0.4 Myr **Right:** Probability distribution function of neutral fraction at 0.4 Myr</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14: **Left:** Map of neutral fraction of $P_n$ over $M_1$ **Right:** Histogram of this ratio for the whole box, fitted with a sum of 3 log-gaussian curve, at 0.4 Myr</figcaption></figure>

obvious that $M_1$ under-ionises cells even as it stands above $P_n$ orders, which have all converged toward a lower proportion of neutral cells. This graph shows how blatantly bad $M_1$ ends up being once the optically thin regime is reached.

At this point in time, we also physically expect a $1/r^2$ photon density profile around sources, as the fully ionised medium is now transparent to radiations. However, it is quite obvious from Fig. 18 that, while $P_9$ seems to fit that description, $M_1$ is extremely far from it. Here, all of the issues of $M_1$ are shown in plain sight, with large ejecta of energy probably due to the interaction between sources that do not represent the physically expected distribution of photon in that regime, or the unexpected rings of under-densities surrounding the sources,

a known characteristic of the model sometimes referred to as the "dark donut" or the "dark sombrero", and that can be found in real cosmological simulations. This specific artefact is highlighted in Fig. 20, where we can see how the photon density dips around sources in $M_1$, forming said "dark sombrero" shape that is completely unphysical. This phenomenon is ubiquitous in the simulation, as shown in Fig. 21, where we plot all the profiles between source 6 and the six other sources in the simulation. We observe that, as long as two sources are distant enough, the same kind of dips appear in the photon density of $M_1$, where the density of $P_9$ is smooth, as is physically expected at this time in the simulation. This hints toward the fact that this unphysical phenomenon, as well as the ejectas of matter we have mentioned previously, also takes place in more
<figure><img src="image_16.png" /><figcaption>Fig. 15: Maps of the neutral fraction of hydrogen at 4.0 Myr in M_1 (Left) and P_9 (Right)</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 16: Mean neutral fraction evolution over time</figcaption></figure>

complete simulations, and, as mentioned previously with the stark differences in neutral fraction between the models. This is also highlighted in the PDF on the left plot of Fig. 17, where a tail appears in the M_1 distribution which is absent from all P_n ones, showcasing far more cells at low density, and, as such, a far more diffuse distribution of photons in the box. We can quantify this difference in ionisation by observing the plot and histogram in Fig. 19, showcasing the ratio of neutral fractions x_{HI}(P_9)/x_{HI}(M_1) at 4.0 Myr. The sigma of the log-gaussian is σ = 0.27, which is almost the same as in the transitional regime. This error is bound to have an impact on the global output of the neutral fraction, which is the observable that is used to probe the reionisation epoch. It appears that the results given by the M_1 model for this observable might not represent reality as good as P_9 does.

To get a better insight into the observational equivalent of 710 what is presented in this simplified version of a cosmological simulation, we can observe the transmission spectra of our cells, i.e. the fraction of radiation they allow to pass through depending on the wavelength of light, as a makeshift "mock" Lyman-α for- est spectra. The transmission in a cell is computed by taking the inverse exponential of τ, which is defined as follows (Dijkstra 2019):

$$ \tau = n_H(1-x)dx\sigma_{HI}(T) \quad (30) $$

With:

$$ \sigma_{HI}(T) = 5.9 \times 10^{-14} \times \left(\frac{T}{10^4}\right)^{-\frac{1}{2}} \times 10^{-4} \quad (31) $$

However, to ensure we get a representative sample of how the transmission can evolve depending on the medium it is going 720
<figure><img src="image_18.png" /><figcaption>Fig. 17: Left: Probability distribution function of photon density at 4.0 Myr Right: Probability distribution function of neutral fraction at 4.0 Myr</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 18: Maps of the photon density at 4.0 Myr in M₁ (Left) and P₉ (Right)</figcaption></figure>

through, we take four separate lines of sights (LoS) in our box to extract the spectra from. They are shown in Fig. 22, and were chosen with the following criterions:

- LoS 1 goes through two regions with a large contrast in neutral fraction when comparing the two models, as seen in Fig. 19

- LoS 2 goes through source 7, whose "dark sombrero" artifact is very visible in most of our images and in Fig. 20

- LoS 3 goes through the high photon density region created by the cluster of sources 3 to 6, whose interaction with source 7 is shown in most of our images

- LoS 4 goes through a region with high neutral fraction $x_{HI}$, thus less ionised

The results are shown in Fig. 23 for all of the four LoS. We observe far less transmission in $M_1$ when the LoS goes through

less ionised medium, which shows once again the weaker ionising effect of this model. We also highlight the fact that this difference is strongest in areas with a large difference in neutral fractions between the two models, especially with LoS 1. Comparing it to Fig. 19, we can see regions where P₉ ionises more, where $M_1$'s transmission is lower, and regions where $M_1$ over-ionises (LoS 1, 400-500 h⁻¹kpc), i.e. where $M_1$'s transmission is higher than P₉'s. This latter area has been shown to be the result of non-physical interactions in $M_1$, which in turn has a real, very tangible effect in these spectra. What we also observe is the presence of cells with a transmission value of 0 for both models in three of the four LoS, showing that there are still areas that are opaque to radiation at this epoch. However, they are quite sparse, and more common and larger in $M_1$.
<figure><img src="image_20.png" /><figcaption>Fig. 19: Left: Map of neutral fraction of P<sub>n</sub> over M<sub>1</sub> Right: Histogram of this ratio for the whole box, fitted with a log-gaussian curve, at 4.0 Myr</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 20: Photon density profile around source 7 at 4.0 Myr</figcaption></figure>

**5. Conclusions**

P<sub>n</sub> is an alternative moment based radiative transfer model to
the usual M<sub>1</sub> model, based on a projection of the moments on
the spherical harmonics basis and using a quite simple closure
equation. After showcasing how it can correct some of the most
glaring issues in M<sub>1</sub>, such as the interactions between sources,
but also pointing at its weaknesses, namely its great sensibility
to discontinuities and its ability to output negative photon
densities, we compared the two models using benchmarked
tests created for this purpose. It turns out that even in the more
physical and astro-physical context, P<sub>n</sub> fares better than M<sub>1</sub>. To
top it all, P<sub>n</sub> corrects the glaring issues M<sub>1</sub> showcases in the
most complex test, which is also the closest to real simulations
where the two models could be used, especially in the optically
thin regime. As such, P<sub>n</sub> appears to be a viable replacement for
M<sub>1</sub> from a purely physical standpoint.

By comparing P<sub>9</sub> and M<sub>1</sub> radiation fields in an idealised and
cosmological test case, we highlight a new, thus far unreported
artefact of M<sub>1</sub>, the 'dark sombrero'. A dark sombrero appears
in M<sub>1</sub> solutions as a spherical photon-deficit shell around the
source, at typically 1/3-1/4 of the distance to the next source.
The photon density in dark sombreros can be underestimated
by a factor up to 2-3. They occur in regions where a source's
radiation field connects with that of another source or group
of sources. These basic properties (position and amplitude) of
the dark sombreros may depend on the sources' relative inten-
sities, positions, spatial resolution, although we have not been
able to test this in detail in this study. Moreover, the M<sub>1</sub> larger
scale photon density also exhibits spurious features, enhancing
or reducing photon density in various regions. We use a small
reionisation-like test simulation to characterize the relative error
in hydrogen neutral fractions between M<sub>1</sub> and P<sub>9</sub>. We find that
there is a small difference in the timing of the reionisation of
the test box. Also, in its final state M<sub>1</sub> test box is slightly more
neutral than the P<sub>9</sub> realization. In regions where reionisation is
finished in both models, the relative error is well represented by
a gaussian with a dispersion of 0.27 dex in log<sub>10</sub>(x<sub>HI</sub>). Both as-
pects are likely related to the photons' collisional behaviour in
M<sub>1</sub>.

We also compute the Lyman-α transmission on 4 LoS of this
mini-reionisation box, to produce mock Lyman-α forests. We
find the main transmission differences between M₁ and P₉ are
related mostly to the small difference in neutral fraction between
the two models, with M₁ being slightly more opaque than P₉, ex-
cept where spurious overionisation occurs in M₁. We highlight
that this is mostly an experiment and ideally, a proper, large-
scale, well-resolved, Cosmic Dawn - like simulations should be
performed to quantify these transmission differences in a more
realistic scenario.

Even though P<sub>n</sub>'s may display a form of oscillatory 800
behaviour in presence of either quickly time-varying or angle-
discontinuous sources (such strongly collimated / focused
emission), this can be mitigated by increasing P<sub>n</sub>'s order, at
<figure><img src="image_22.png" /><figcaption>Fig. 21: Photon density profiles between source 6 and six other sources of the simulation at 4.0 Myr</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 22: Positions of the Lines of Sight used for the transmission spectra</figcaption></figure>

the cost of a lot of computational power (see Appendix B) or filtering. We’ve shown that models like P<sub>3</sub> or P<sub>5</sub> might have too low orders to be used in simulations, but that P<sub>7</sub> and P<sub>9</sub> seem already close to convergence in the tests we’ve showcased. We consider that there is probably no need to push the order further for our simulations, which, in general, fit the model’s needs quite well, with their isotropic sources and simple geometry. However, further tests could be done comparing P<sub>9</sub>, P<sub>11</sub>, P<sub>13</sub> and P<sub>15</sub> to ensure this convergence in more complex environments.

Of course, all the tests presented in this paper were done in
a simplified context, with only hydrogen, simplified chemical
reactions, and no gravity or hydrodynamics added. Knowing

that these absent processes, especially gravity, already use up
most of the computing power and time in current cosmological
simulations, while the use of M<sub>1</sub> allowed the budget of the
radiative transfer to be relatively low, the question of the cost
P<sub>n</sub> could add to these simulations has to be asked. Furthermore,
even if we’ve highlighted some of the starkest differences
between M<sub>1</sub> and P<sub>n</sub>, we still need to study these differences and
their impact in larger, more complex simulations. As such, a
comparison should be made between fully fledged cosmological
runs with both models to ensure this increased computational
costs translates to a significant variation in our final results.
<figure><img src="image_24.png" /><figcaption>Fig. 23: Transmission spectra along several lines of sight for $M_1$ and $P_9$</figcaption></figure>

*Acknowledgements.* The authors would like to acknowledge the High Performance Computing Center of the University of Strasbourg for supporting this work by providing scientific support and access to computing resources. Part of the computing resources were funded by the Equipex Equip@Meso project (P850) gramme Investissements d'Avenir) and the CPER Alsac calcul/Big Data.

This work of the Interdisciplinary Thematic Institute IRMIA++, as part of the ITI 2021-2028 program of the University of Strasbourg, CNRS and Inserm, was supported by IdEx Unistra (ANR-10-IDEX-0002), and by SFRI-STRAT'US project (ANR-20-SFRI-0012) under the framework of the French Investments for the Future Program.

References

Aubert, D. & Teyssier, R. 2008, Monthly Notices of the Royal Astronomical Society, 387, 295

Barkana, R. & Loeb, A. 2001, Phys. Rep., 349, 125

Brunner, T. A. 2002, Forms of Approximate Radiation Transport, Tech. rep., Sandia National Lab. (SNL-NM), Albuquerque, NM (United States); Sandia National Lab. (SNL-CA), Livermore, CA (United States)

BRUNNER, T. A. 2002

Dewdney, P. E., Hall, P. J., Schilizzi, R. T., & Lazio, T. J. L. W. 2009, Proceedings of the IEEE, 97, 1482

Dijkstra, M. 2019, Saas-Fee Advanced Course, 46, 1

Garrett, C. K. & Hauck, C. D. 2013, Transport Theory and Statistical Physics, 42, 203

Gerhard, P. 2020, PhD thesis, thèse de doctorat dirigée par Helluy, Philippe Mathématiques Strasbourg 2020

Hui, L. & Gnedin, N. Y. 1997, Monthly Notices of the Royal Astronomical Society, 292, 27

Iliev, I. T., Ciardi, B., Alvarez, M. A., et al. 2006, Monthly Notices of the Royal Astronomical Society, 371, 1057

Katz, H., Kimm, T., Haehnelt, M. G., et al. 2019, MNRAS, 483, 1029

Levermore, C. 1984, Journal of Quantitative Spectroscopy and Radiative Transfer, 31, 149

Lewis, J. S. W., Ocvirk, P., Aubert, D., et al. 2020, MNRAS, 496, 4342

Maselli, A., Ferrara, A., & Ciardi, B. 2003, Monthly Notices of the Royal Astronomical Society, 345, 379

Meltz, B. 2015, PhD thesis, University of Paris-Saclay, France

Naidu, R. P., Tacchella, S., Mason, C. A., et al. 2020, ApJ, 892, 109

Ocvirk, P., Lewis, J. S. W., Gillet, N., et al. 2021, MNRAS, 507, 6108

Osterbrock, D. E. 1974, Astrophysics of gaseous nebulae

Rosdahl, J., Blaizot, J., Aubert, D., Stranex, T., & Teyssier, R. 2013, MNRAS, 436, 2188

860
870

Rusanov, V. 1962, USSR Computational Mathematics and Mathematical Physics, 1, 304

Sahmim, S. 2005, Theses, Université Paris-Nord - Paris XIII, mme Anela Kumbaro, Mme Laure Quivy, M. François Alouges (rapporteur), M. Claude Basdevant (président), M. Fayssal Benkhaldoun (Directeur), M. Hervé Guillard (Rapporteur)

Wu, X., McQuinn, M., & Eisenstein, D. 2021, Journal of Cosmology and Astroparticle Physics, 2021, 042
## Appendix A: P_n matrices

In this appendix, we will write the definitions of the matrices $J_x, J_y, J_z$ as defined in the thesis of Bertrand Meltz (Meltz 2015). We won't however go through the extensive proof and recurrence developed in said paper. Let us define $A_l^m$ and $B_l^m$ in the context of spherical harmonics with $l < n, m \in [-l; l]$, as follows:

$$A_l^m = \sqrt{\frac{(l-m)(l+m)}{(2l+1)(2l-1)}} \quad B_l^m = \sqrt{\frac{(l+m-1)(l+m)}{(2l+1)(2l-1)}} \quad (A.1)$$

Then, to each orbital $(l, m)$ we associate an positional index defined as:

$$i(l_i, m_i) = \sum_{l=0}^{l_i-1} \sum_{m=-l}^{l} 1 + \sum_{m=-l_i}^{m_i} 1 = l_i^2 + l_i + m_i + 1 \quad (A.2)$$

Counting all possible orbitals for a model of order n gives us $(n+1)^2$ coefficients, which is the size of our vector and also, consequently the size of our matrices. With $(i, j) \in [1, (n+1)^2]^2$ indices in the matrices each corresponding to an orbital as defined in Eq. A.2, we can define the matrices as follow:

$$\begin{aligned} J_{i,j}^x = \frac{sgn(m_i)}{2} & (1 + (\sqrt{2} - 1)\delta_{m_i,1}) \\ & \times [-\delta_{l_j,l_i-1}\delta_{m_j,m_i-1}B_{l_i}^{m_i} + \delta_{l_j,l_i+1}\delta_{m_j,m_i-1}B_{l_j}^{-m_i}] \\ & + \frac{sgn(m_i + 1)}{2} (1 + (\sqrt{2} - 1)\delta_{m_i,0}) \\ & \times [\delta_{l_j,l_i-1}\delta_{m_j,m_i+1}B_{l_i}^{-m_i} - \delta_{l_j,l_i+1}\delta_{m_j,m_i+1}B_{l_j}^{m_i}] \end{aligned} \quad (A.3)$$

$$\begin{aligned} J_{i,j}^y = \frac{sgn(m_i)}{2} & (1 - \delta_{m_i,1}) \\ & \times [\delta_{l_j,l_i-1}\delta_{m_j,-(m_i-1)}B_{l_i}^{m_i} - \delta_{l_j,l_i+1}\delta_{m_j,-(m_i-1)}B_{l_j}^{-m_i}] \\ & + \frac{sgn(m_i + 1/2)}{2} (1 + (\sqrt{2} - 1)(\delta_{m_i,0} + \delta_{m_i,-1})) \\ & \times [\delta_{l_j,l_i-1}\delta_{m_j,-(m_i+1)}B_{l_i}^{-m_i} - \delta_{l_j,l_i+1}\delta_{m_j,-(m_i+1)}B_{l_j}^{m_i}] \end{aligned} \quad (A.4)$$

$$J_{i,j}^z = \delta_{l_j,l_i-1}\delta_{m_j,m_i}A_{l_i}^{m_i} + \delta_{l_j,l_i+1}\delta_{m_j,m_i}A_{l_j}^{m_i} \quad (A.5)$$

Where sgn is a sign function defined as follows:

$$\forall x \in \mathbb{R}, \quad (A.6)$$

$$sgn(x) = \begin{cases} 1, & x > 0 \\ 0, & x = 0 \\ -1, & x < 0 \end{cases} \quad (A.7)$$

## Appendix B: Order and negativity

As mentioned previously, one of P_n's limitations is its propensity to oscillate in the presence of strong spatial, angular and time discontinuities. Source angular discontinuities are quite unusual in the context of reionisation simulation, since sources are in general isotropic and continuous, and the boundary conditions are generally periodic. Here, we tried very hard to put $P_{20}$ through the worst case scenario and see how it fares.

<figure><img src="image_25.png" /><figcaption>Fig. A.1: Pulse response of M₁ Left: and P₃ Right:</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. A.2: Radial profiles of an impulse response using M₁ and various orders of P_n</figcaption></figure>

This short test consists of a single discontinuous adimensional source. There is no chemistry involved here, we just want to observe how P_n reacts to a pulse-like source that we will define as a dirac, such that the integral of the cell containing it is equal to 1. As such, our source emits as follows during one time step:

$$S = \frac{1}{dxdydz} \quad (B.1)$$

We expect the differences between M₁ and P_n to be quite sharp, with M₁ being strictly positive and non oscillatory, contrary to P_n. However, we're interested in the difference between the various orders of P_n and at what point the order is high enough for these oscillations to dampen enough to be comparable to M₁.

These results can already be seen quite clearly in fig A.1 where the ring of the wave front for M₁ is very visible, but P₃ oscillates and outputs negative energies, which makes the position of the wave front less obvious. This is even more blatant when looking at the profiles for both models, as P₃'s minimum at this time step is -40, its negative amplitude being far superior to its positive part. It is important to note that P₃ is used here on purpose, as its very low order showcases the limits of P_n the most, and demonstrates why there is a need to go to
<figure><img src="image_27.png" /><figcaption>Fig. A.3: Dense clump of hydrogen in the path of a flux of ionising radiations using $P_3$ at 16 Myr</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Fig. A.4: Fraction of negative photon density cells Left: and minimum density Right: over time in the cosmological map test</figcaption></figure>

higher orders for a sensible result, which we did during this paper.

To convince ourselves further, a simple look at Fig. A.2 shows the various profiles of the impulse response for several orders of $P_n$, with $M_1$ in dashed line as a reference. The oscillations of the models are very strong under order 9, where the output becomes almost indistinguishable from $M_1$. As such, we can assume that orders above $P_9$ are able to handle the worst case of $P_n$ the best, as their oscillations are dampened enough for them to be negligible. This is the reason why we focused mainly on $P_9$ in this paper, instead of less expensive orders like $P_7$ or $P_5$. We also point out that it is possible to use filters to dampen the oscillations for almost no computational cost, but it is beyond the scope of this paper.

The issue of the lower orders of $P_n$ appears in other unfavourable cases, for example the clump test of section 4.3, where all sources are directional, thus non isotropic, thus introducing a discontinuity into $P_n$. Fig. A.3 shows the result of this test using $P_3$, and we can clearly note that $P_n$, at low order and similar computational cost to $M_1$, cannot reproduce correctly what we physically expect from this test. $P_3$, which only requires 16 coefficients, creates large visible modes in the photon density, and tends to over-ionise the sphere. On top of that, its modes output negative energies in a large area close to the sources, at the top of the simulation box, which creates artefacts in the neutral fraction and temperature that should not exist. As such, we can say that, to have a reliable $P_n$ model

able to reproduce or outperform $M_1$'s result, a higher order is required, and as such, a higher computational cost.

However, can we ensure the positivity of our model at all times, even at higher orders and in a favourable test? This can be answered partially by Fig. A.4, where we plot the fraction of negative cells to the right in the cosmological map test. One surprising result is that, in this specific test case, the fraction of negative cells seems to be almost independent of the order of our model, with a peak around 7.5% of cells around 0.2 Myr. It means that no matter which model we use, there will be at least some cells outputting a negative energy, which might pose a problem. However, we can note that in the case of this test, these negative cells are transient and do not last until the end of the simulation, but also that these cells are not as negative as others. Indeed, when looking at the right plot of Fig. A.4 that shows the minimum energy value at each time step, it appears obvious that all models above $P_7$ have an almost negligible negative outputs, which in turn won't impact the simulation as much as $P_3$ or $P_5$ would. Indeed, since we approximate the value of a negative cell as 0 (see Eq.18), the impact of this modification will be smaller if the value of our negative cell is already almost a 0. This is yet another reason for favouring higher order $P_n$ even in scenarios where lower orders would seem accurate.