The IDSA and the homogeneous sphere: Issues and possible improvements

In this paper, we are concerned with the study of the Isotropic Diffusion Source Approximation (IDSA) (Baxter et al., Phys. Rev. E 73, 046118, 2006) of radiative transfer. After having recalled well-known limits of the radiative transfer equation, we present the IDSA and adapt it to the case of the homogeneous sphere. We then show that for this example the IDSA suffers from severe numerical difficulties. We argue that these difficulties originate in the min-max switch coupling mechanism used in the IDSA. To overcome this problem we reformulate the IDSA to avoid the problematic coupling. This allows us to access the modeling error of the IDSA for the homogeneous sphere test case. The IDSA is shown to overestimate the streaming component, hence we propose a new version of the IDSA which is numerically shown to be more accurate than the old one. Analytical results and numerical tests are provided to support the accuracy of the new proposed approximation.

1. Introduction. In many astrophysical problems radiation plays a crucial role.
As an example, core-collapse supernovae (CCSNe) explosion mechanisms crucially depend on the radiative transfer of neutrinos. In such practical calculation, the numerical computation of the full transfer equation is usually too costly and one has to rely on approximations.
In the context of CCSNe, Liebendörfer et al. [6] have designed the Isotropic Diffusion Source Approximation (IDSA). This approximation is based on asymptotic solutions of the radiative transfer equations, see [2] for a rigorous derivation of these asymptotic solutions; the coupling is then performed with the so-called Diffusion Source Σ which has the form of a min-max switch. In [2], it is shown that this particular coupling mechanism might be problematic and mathematical issues are discussed. The aim of this paper is to study the behavior of the IDSA on a simple radiative transfer model and try to address some of the mathematical issues highlighted in [2].
The paper is organized as follows. In Section 2 we present the simple radiative transfer model to which we apply the IDSA. We also present the first three angular moments equation associated to it together with the well-known asymptotic limits of these equations, that is the diffusion limit, the reaction limit and the free streaming limit. Section 3 is devoted to the presentation of the IDSA adapted to our model problem. We discuss the three regimes of this approximation together with mathematical issues illustrated by numerical experiments. In Section 4, we discuss the homogeneous sphere test case and its limit of high opacities. We then reformulate the IDSA in order to remove the diffusion source Σ and propose some modifications of the IDSA equations. In Section 5, we perform numerical experiments on the two new formulations of the IDSA and study their modeling error in the case of the homogeneous sphere. We conclude in Section 6 by discussing the different ways this work can be used in practical simulations of CCSNe explosions.
2. Radiative transfer equation and asymptotic limits. In this section, we consider the monochromatic (independent of energy) radiative transfer equation in spherical symmetry. We assume the background to be static and exclude inelastic and anisotropic scattering processes. The radiative transport equation in natural units (c = 1) reads where f (t, r, µ) is the distribution function of the transported particles that depends on the time t, radius r and cosine of the angle between the propagation direction and the radial direction µ; κ a (r) and κ s (r) are the absorption and scattering opacities, respectively; b(r, µ) is the equilibrium distribution and J(t, r) is the angular average of f (t, r, µ) defined in Eq. (4). We also introduce the linear transport operator L f of f . The first two moment equations are obtained by multiplying Eq. (1) by 1 and µ respectively and integrating over µ. This gives L J (J, H) := ∂ t J + 1 r 2 ∂ r (r 2 H) = κ a (B − J), and L H (J, H, K) := ∂ t H + 1 r 2 ∂ r (r 2 K) + The two moments equations (2) and (3) basically express the conservation of energy and momentum. The three angular moments J, H and K of the radiation field are proportional to the radiation energy density, flux and pressure. We also used the integrated equilibrium distribution B = 1 2 1 −1 b dµ and defined the transport operators L J and L H for J and H, respectively. Note that these operators also depends on higher order moments and we obtain in general a hierarchy of equations that needs to be closed by a prescription of the form for a one-moment or a two-moment closure, respectively. A classical two-moment closure has been introduced in [5] and other two-moment closures are discussed in [9] in relation with flux-limited diffusion. It is also convenient to define the following flux factors.
Definition 2.1 (Flux ratio and variable Eddington factor). The flux ratio h(t, r) is defined by and the variable Eddington factor k(t, r) is defined by k(t, r) := K(t, r) J(t, r) .
In order to enable a complete discussion of the IDSA, we need to discuss the asymptotic behavior of the radiative transfer equation because the IDSA is based on the coupling of three different asymptotic limits of this equation as shown in [2]. We discuss the diffusion limit, the reaction limit and the free streaming limit in terms of closure relations similar to Eq. (5). 2.1. Diffusion limit. Following [8, pp. 350-353], we can use the diffusion approximation of (1) by imposing the closure conditions where κ := κ a + κ s is the total opacity. These conditions lead to the diffusion equation where we also defined the diffusive transport operator L diff J for J. In this limit, the variable Eddington factor satisfies The diffusion limit is valid in regions where the mean free path κ −1 is short compared with the typical length scale of the system. Such a region will be refered to as opaque. This limit can also be obtained by asymptotic expansions, see for example [2].
where the free streaming flux ratio can be estimated by a geometrical relation under the assumption that the streaming particles are emitted from an infinitely opaque homogeneous sphere of radius R. In this case, we have [3,6] Assumption 1 (Scattering sphere approximation). The scattering sphere onemoment closure relation is given by the flux ratio The radius R is called the neutrinosphere radius.
This flux ratio corresponds, as we will rederive in Section 4, to the flux ratio of an infinitely opaque sphere of radius R. In Section 4 we also compute the corresponding variable Eddington factor. The neutrinosphere radius R is usually computed using the optical depth τ and the neutrinosphere radius satisfies τ (R) = 2 3 . As for the diffusion and the reaction limit, the stationary state free streaming limit can be derived by asymptotic expansions, see [2]. This approximation is valid where the mean free path κ −1 is large compared with the typical length scale of the system. In these regions, the particles propagate freely and we therefore refer to them as transparent regions. In term of flux factor, the free streaming limit satisfies lim r→∞ h free (r, R) = 1, lim r→∞ k free (r, R) = 1.
These limits come from the fact that in this case, the distribution function is strongly outward peaked and accumulates on µ = 1. THE IDSA AND THE HOMOGENEOUS SPHERE   5 3. Isotropic Diffusion Source Approximation (IDSA). As mentioned in Section 1, in practical calculations the direct solve of Eq. (1) is too costly and one need to rely on approximations. The Isotropic Diffusion Source Approximation (IDSA) has been developed by Liebendörfer et al. in [6] in the context of core-collapse supernovae (CCSNe) modeling. This is an approximation of the radiative transfer of neutrinos in this kind of stellar explosion. It is based on the observation that in a CCSN event, the inner core of the star becomes dense enough to be opaque to neutrinos, whereas the outer layers of the star are transparent. Following the discussion in Section 2, the radiative transfer of neutrinos in the inner core should be well described by the diffusion or the reaction limit, whereas in the outer layers it should be well described by the free streaming limit. This observation has led to the idea of using an additive decomposition of the distribution function into two components: a trapped component describing the inner core neutrinos and a free streaming component describing the outer layers neutrinos. Liebendörfer and colleagues also noted that it was enough, for their simulations, to compute the angular averaged J distribution function. The IDSA is therefore an heterogeneous approximation of the first angular moment J of the neutrino distribution function f . We now recall the hypotheses of the IDSA and state it for our model problem.
3.1. Ansatz: Decomposition into trapped and streaming neutrinos. We assume a decomposition of on the whole domain into distribution functions f t and f s supposed to account for trapped and for streaming neutrinos, respectively. This ansatz is also valid for all the moments by linearity of integration where we used the notation introduced in Eq. (4). The IDSA is a two components one-moment approximation of Eq. (2) that we recall here Assumption 2. Eq. (21) can be well approximated by a system of equation of the form where Σ is the diffusion source, and J t + J s ≈ J.
The IDSA is completed by defining the operators L t J and L s J together with the closure relations for the trapped and the streaming components.

Assumption 3.
We assume that f t , κ a and κ s are isotropic and that f s is in the stationary state free streaming limit. Hence, the IDSA uses the following closure relations , The IDSA is defined by Eqs (22)-(24) together with the closure relations given in Eq. (25) and the operator definitions given in Eq. (26).
The diffusion source Σ is defined in a way to switch between three different regimes to which we now come.

3.2.
The three regimes of the IDSA. The IDSA has three distinct regimes corresponding to the different values that the diffusion source Σ can take: 1. Reaction regime: Σ = 0; 2. Diffusion regime: 3. Free streaming regime: Σ = κ a B.
We now discuss these three regimes.
Remark 1. In this work, when speaking about the behavior of the IDSA, we will speak about regimes. But when speaking about the asymptotic behavior of the real distribution, we will speak about limits. Hence, it is important to distinguish clearly between these two terms.
3.2.1. Reaction regime. In the reaction regime, the IDSA equations (27) and (28) become which shows that the trapped component is in the reaction limit, see Eq. (12) and J s reac is solution of Eq. (31), that is Because 1 2 ≤ h free (r) ≤ 1, J s reac (r) > C e −r 2 r 2 which diverges in 0 in an unphysical way. Therefore, we have C = 0 and and the solution of the IDSA in the reaction regime is in the reaction limit.
or equivalently that is the trapped component is in the diffusion limit up to a term of absorption of streaming particles. The first equation is obtained by passing the diffusion term on the left hand side, comparing with Eq. (9). The second is obtained by multiplying Eq. (35) by r 2 and integrating over r. This modified diffusion limit can also be obtained by asymptotic expansions, this has been shown in [2].
Remark 2. Note that in order to have a positive distribution of streaming particles J s diff ≥ 0, it is necessary to have a decreasing distribution function of trapped ∂ r J t diff ≤ 0. That is, any non-decreasing solution of Eq. (36) is not physical. Substituing (37) into (36), one can decouple the two equations and the new equation for the trapped particles is given by where we explicit the diffusion transport operator L diff J . 3.2.3. Free streaming regime. In the free streaming regime, the IDSA equations (27) and (28) become which shows that the streaming component is in the stationary free streaming limit, see Eq. (13). The trapped equation can be solve analytically, we obtain That is the trapped component is exponentially decreasing and vanishes in the stationary state.

Mathematical and numerical issues.
We have discussed in Subsection 3.2 that the different regimes of the IDSA correspond to the asymptotic limit of the radiative transfer equation. The main problem of the IDSA is its treatment of transient regions where neither the diffusion, the free streaming nor the reaction limit apply. Then, nevertheless, one of the above regimes occurs. In the case of CCSN modeling, a discussion of this topic by numerical results comparing the IDSA with the full radiative transfer equation in the spherically symmetric case is given in [6]. In this paper, it is shown that the IDSA solution approximates pretty well the reference solution computed by solving the full radiative transfer equation. It is also shown that the biggest deviations occur in the "semi-transparent" transient region around the neutrinosphere. Similar results have been founded in a numerical example for a model problem in [1]. The limiters for the coupling term in (24) can be physically and mathematically motivated, see for example [2,6], but there does not seem to be a rigorous explanation of these limitations by means of the full radiative transfer equation.
In the following, we take a closer look at some properties of the diffusion source Σ and show that this coupling can lead to unphysical solutions and numerical instabilities.
3.3.1. Spurious trapped. The first problem that we discuss here is referred to as the Spurious Trapped problem. This is a problem in modeling that occurs in regions with small but non zero opacities (κ small). The problem originates in the definition of the coupling term Σ. Such a region should be transparent and therefore described by the free streaming limit. Hence, the IDSA limiter should always choose the free streaming regime in this case. We now construct an example in which the definition of diffusion source Σ leads to a wrong regime. Proposition 1. In a transparent region where κ a = O(ε) and where the distribution J t is locally constant and J s < B, we have that is the IDSA does not choose the free streaming regime in this case.
Proof. The condition on κ a is needed in order to have a transparent region, for more details about this asymptotic limit, see [2,7]. We now compute the value of the diffusion source in this case.
where we used the fact that J t is locally constant. We now observe that κ a > 0 and that J s < B by hypothesis, this leads to This concludes the proof.
We now discuss the relevance of the hypotheses used in Proposition 1. As we have mentioned in the proof, the assumption on κ a just reflects the fact that we are in a transparent region. In such a region, Eq. (41) implies that the trapped distribution vanishes. Finally, the condition on J s is natural in radiative transfer problems. In fact, the equilibrium distribution is mainly driven by emission and absorption processes which are small in a transparent region. The transport processes therefore tend to reduce the distribution which leads to J s < B.
Now that we know that the hypotheses can be realized, we look at the resulting dynamics of the trapped particles Because B − J s > 0, we can assume that there exists B such that B − J s ≥ B > 0 and that is, the trapped particle distribution is growing.       Remark 3. Note that if κ a = 0 in the free streaming regime, a non-vanishing trapped distribution will not be eliminated because its evolution becomes and the distribution, whatever it is at first, is stationary. Now that we have seen that the spurious trapped problem may occur, we show a numerical illustration of this problem. For the discretization, we used the discretization proposed by Liebendörfer in [6] described in more detail in [1,7]. For the grid parameters, we choose N r = 50 and ∆t = 0.1. For the parameters of the radiative transfer equation, we choose B = 1, κ s = 0 and (48) The results are displayed in Figure 1. In this Figure, we also show the proportions h t = J t J t +J s and h s = J s J t +J s of trapped and streaming particles on the right panels. This illustrates the fact that the trapped component is growing. This problem develops on a slow timescale, and becomes visible only on a long timescale. When ε goes to zero, the time needed to create spurious trapped particles tends to infinity and this problem disappears. In Figure 2, we display the time needed to reach a stationary distribution dominated by spurious trapped particles in the streaming region. The results show that the time needed to reach a stationary state grows roughly with an order of O(ε −1 ). Numerically, a linear regression on the results presented in Figure 2, where we excluded the 5 first points because for too big ε the free streaming limit is no longer a good approximation, gives an order of O(ε −0.90 ). It is reasonable to think that inside a CCSN, the time interval over which the solution is computed is too short for this problem to develop.

3.3.2.
Instability of the coupling. The second problem we discuss is also linked with the diffusion source Σ. It takes the form of an instability that develops at a boundary between an opaque and a transparent region. It originates in the fact that the diffusion term in Σ can become large for two reasons: the first is because κ is small and the domain is transparent, and this triggers a correct transition to the streaming limit; the second is because the concavity becomes large. At the boundary between a very opaque and a transparent region, the second case may occur leading to a switch to the streaming regime in a region where the opacity κ is large. That is, where the diffusion limit is valid and the diffusion regime of the IDSA should apply. As a consequence, the IDSA virtual boundary between the opaque and the transparent domains is shifted inward. This shift is helped by the creation of     limit is recovered near the real boundary and starts to grow again as expected, but at the virtual boundary of the IDSA, the trapped distribution continues to locally decrease leading to a locally radially growing distribution of trapped particles which is not physical, see Remark 2. As an effect, an instability of the underlying model is created and propagated inward.
As a numerical illustration of such an instability, we apply the IDSA to the homogeneous sphere test case defined by κ s (r) = 0, κ a (r) = 1 r<6 and B = 1. For the grid parameters, we set N r = 10 4 and ∆t = 0.1. Figure 3 displays the formation of the instability near the boundary of the homogeneous sphere. We see for example in the first subfigure that the virtual boundary of the IDSA, where the trapped distribution vanishes, stands at a radius r < R. This illustrates the shift of the interface discussed above. The next subfigures show the early developement of the instability and its inward transport due to the feedback of the streaming component into the trapped component dynamics, see Eq. (38).
The behavior of the solution on longer times is displayed in Figure 4. This figure shows the evolution of the instability. The solution becomes more complicated as time passes. Note that the solution remains bounded and never explodes like a usual numerical instability due to numerical problems such as a CFL condition. This instability in the solution comes from the modeling and the definition of the coupling term Σ.
Remark 4. The instability is linked with the size of the diffusion term. Note that this quantity also depends on the discretization as the second derivative will be of the order of O((∆r) −2 ). This means that if the discretization is coarse enough, then the instability does not occur. This explains the fact that we used a coarse grid in the spurious trapped numerical experiment. In fact, if we refine the grid, the spurious trapped numerical experiment will also show the numerical instability that we just described.
We have shown with these two examples that the coupling in the IDSA equations can lead to wrong solutions that are due to the application of a wrong asymptotic regime in parts of the domain. This can lead to the creation of spurious trapped particles or to the creation of an instability in the solution. We have also discussed the influence of the advection term coming from the coupling with streaming particles in the trapped diffusion equation. It seems that this term helps the developement of the instability.
A natural question to ask is: why did the IDSA work in the case of the CCSN models performed in [6]? We answer this question below.

3.4.
Why does the IDSA work in practice? As we mentioned in Section 1, in the numerical tests performed by Liebendörfer et al. in [6], the IDSA is quite accurate and good enough to be used in simulations. We just discussed two important mathematical issues of this approximation. How is it possible that the IDSA works in practice when it fails on simple examples? Part of the answer has been given in Remark 4. If the discretization is coarse enough, the instability can not develop. In the application of the IDSA to the CCSN modeling, the number of grid points is N r = 1000 and the repartition of points is such that the instability of the underlying model does not develop. However, it is probable that a refinement of the radial grid will lead to the apparition of an instability in the IDSA. Note that in a real CCSN simulation, the reactions in the center become very large and as a result, the instability is quickly damped. In our example for the instability, we used an artificially small opacity κ = 1. For larger opacities, the instability is damped and vanishes in the limit of the infinitely opaque homogeneous sphere. Another possible explanation is that in the CCSN numerical experiments, the radiative transfer equation is coupled to the background matter. This coupling is done through the trapped particles and their distribution is reconstructed at every time step with a smoothing effect that may help to get rid of the instability. Regarding the spurious trapped particles, we have observed that the growing of the spurious trapped particles component is slow and that it takes a long time to dominate the dynamics, if ever. In fact, in a CCSN simulation, the time window on which the neutrino radiative transfer needs to be computed is short and the spurious trapped particles do not contribute significantly to the dynamics.
We now study in full detail the case of the homogeneous sphere where an analytical solution is known. This idealized example will shed light on the coupling mechanism and will help us to partly correct or avoid the problems we have described. 4. The case of the stationary homogeneous sphere. In this example we want to solve the monochromatic equation (1) with κ s = 0, that is We assume that the equilibrium distribution B is constant. The homogeneous sphere setting is completed by requiring a specific form for κ a . We set where κ is a constant and 1 is the characteristic function. This example has been chosen for two main reasons. The first is because it has a closed form solution that is for example derived in Chapter 2 of Duderstadt and Martin [4]. The second reason is that this example naturally has two different regions where the limit dynamics correspond to the limits used in the derivation of the IDSA. It is therefore a very good and natural case study for the IDSA. The steady state radiative transfer equation (49) has an analytical solution given by where and G(r, µ) : This solution is easily obtained from the formal solution of the transport equation; see e.g. [4].
The first three angular moments of the analytical exact solution are given by and and and can be computed by direct integration from Equation (51) using the parity of the function G(r, µ) with respect to µ. We now give some particular values of the analytical solutions that will be useful for the analysis of the IDSA. (60) Proof. The proof is done by direct integration.
Before applying the IDSA to the homogeneous sphere, we discuss the limit of the homogeneous sphere when the opacity κ → ∞. We call this limit the infinitely opaque homogeneous sphere.
For the stationary distribution function (51), we have the following limit For the first three angular moments, we have and and These results directly come from Eqs (54)-(56) by noting that the dependence on κ is contained only in the integral terms that vanish when κ → ∞ as shown in the following proposition.
Proof. The proof for I 3 is clear as the only dependence on κ is in the exponential. The proof for I 1 and I 2 relies on the fact that the integrants satisfy µ sinh(κrµ)e −κRG(r,µ) = O(e κ(rµ−RG(r,µ)) ), because the behavior of cosh and sinh is dominated by their positive exponential parts. For I 1 and I 2 to vanish in the limit κ → ∞, it is sufficient to check that rµ − RG(r, µ) < 0 for r < R, or equivalently, because all these quantities are positive, that (rµ) 2 − (RG(r, µ)) 2 < 0 for r < R. We obtain This concludes the proof.
In this limit, we can compute the variable Eddington factor k = K J and the flux factor h = H J defined in Definition 2.1.

Proposition 4.
In the limit κ → ∞ the flux factor h R and the variable Eddington factor k R of the homogeneous sphere of radius R are given by and (72) Proof. The proof is done by direct computation from equations (62)-(64).
With these results on the behavior of the analytical solution of the homogeneous sphere, we can now study the IDSA in this particular case.

4.1.
Applying the IDSA to the homogeneous sphere. In order to apply the IDSA to the homogeneous sphere test case and to overcome the numerical problems illustrated in Subsection 3.3, we need to reformulate the IDSA. This is done by introducing the following assumption.
Assumption 4. We assume that the radius R of the homogeneous sphere corresponds to the neutrinosphere radius. We also assume that 1. The domain of validity of the diffusion and the reaction regimes is characterized by r < R. 2. The domain of validity of the free streaming regime is characterized by r ≥ R. This is a very important assumption because it allows us to avoid the problematic diffusion source Σ and to realize the coupling directly and explicitly. This way, we eliminate both the spurious trapped and the instability problem discussed above and we will have access to the modeling error of the IDSA.
Remark 5. Assumption 4 allows us to simplify the Eqs (31) and (32) for the streaming particles and (35) and (37) for the trapped particles, recalling (13), because the function h free = g idsa = 1 2 in these cases.

Remark 6.
We remark here that Assumption 4 does not match the usual definition of the neutrinosphere based on the optical depth, see Eq. (15). In fact, in the case of the homogeneous sphere, one can compute the neutrinosphere radius R ν where we used the definition of κ a , see Eq. (50); that is, the usual definition of the neutrinosphere radius only corresponds to the homogeneous sphere radius in the limit κ → ∞ as R − R ν = 2 3κ .

JÉRÔME MICHAUD
Remark 7 (IDSA in the limit κ → ∞). In the limit κ → ∞, the diffusion limit is irrelevant and we have the following solution If we set C(t) := R 2 B 4 , the solution of this coupling is equivalent to the analytical solution (62) of the infinitely opaque homogeneous sphere. This condition can be ensured by requiring the boundary condition J s free (R, t) = B 2 . Using Assumption 4 and Remark 5, one can rewrite the IDSA equations (37)-(40) without reaction and with κ = 0 for r ≥ R as with g(r) = 1 2 1 + 1 − R r 2 the geometrical factor appearing in Eq. (14).
It is possible to reconstruct the first and second angular moments H and K by defining flux factors for the trapped and the free streaming components. We set and Remark 8. We note that the definition of h s R corresponds to the definition of h free = g idsa , see Eq. (14). In fact, this definition can be used as a justification for the definition of g idsa in the derivation of the IDSA. The definition of these flux factors extend the stationary free streaming limit to the first order moment solution.
Proposition 5. With Equations (79) and (80), we can reconstruct H idsa (r) and K idsa (r) through the relations This reconstruction has the following properties: 1. IDSA closure in the diffusion regime: 2. Free streaming flux ratio: 3. Diffusion variable Eddington factor:

Free streaming variable Eddington factor:
where h R and k R are the flux factors for the infinitely opaque homogeneous sphere and h idsa R := H idsa J t +J s and k idsa R := K idsa J t +J s are the IDSA flux factors.
Proof. We prove separatly the four points.
1. The first statement can be proved by inserting Eqs (79) and (78) into Eq. (81). 2. In the region r ≥ R, we have J t = 0. As a consequence, h s R = h idsa R = h R by Eqs (79) and (71). 3. In the region r < R, we have k t R = k s R = 1/3. Therefore, This concludes the proof.
The proposition 5 shows that the flux factors of the IDSA correspond to the flux factors of the infinitely opaque homogeneous sphere. However, the flux ratio of the IDSA in the diffusion region r < R does not match the diffusion limit relation (8) even if the form is similar. In fact, where we used Eq. (78) 2 to eliminate the streaming component in the last equality. This means that the diffusion limit is not well represented by this formulation. Moreover, the diffusion equation (38) used for the trapped particles does not correspond to the usual diffusion approximation equation (9). Another possible problem of this formulation is that the amount of streaming particles at the homogeneous sphere radius is not necessarily accurate. If this value is inaccurate, then the error will be propagated by the free streaming equation in the whole domain where r ≥ R. This can create large approximation errors.
In order to deal with the problem of the advection term and the possible inaccuracy of the amount of streaming particles at the homogeneous sphere radius, we propose in the next subsection some possible modifications of the approximation.

4.2.
Improvements of the scheme. We now propose some improvements of the IDSA. We propose to use the diffusion equation (9) instead of Eq. (78) 1 for the evolution of the trapped particles and to modify the reconstruction of the streaming particles in the diffusion regime to correctly predict the amount of streaming particles at the homogeneous sphere radius. This is done through the following assumption Assumption 5. In the diffusion region, we assume that the free streaming component can be reconstructed from the trapped component by a relation of the form The constant C is chosen in such a way that J s (R) = J(R) , see Proposition 2.
The scheme obtained can be written as where g(r) = The scheme ensures that the trapped particles are evolved by a diffusion equation and that the amount of streaming particles is correct in the whole domain where r ≥ R. Another advantage of this formulation is that it has a stationary state closed form solution.
In the stationary state limit, we can compute the closed form solution of system (85) given in the following proposition. Proposition 6. The closed form solution of the stationary state limit of (85) is given by where Proof. To prove this proposition, it is enough to substitute the proposed solution into system (85) and checking that it is a solution. This comparison leads to the following result.
Proposition 7. The relative error at r = 0 of the New IDSA approximation (86) is given by .
The value of Err 0 can be used as a measure of the validity of the New IDSA approximation (85). Figure 5 shows that the error vanishes when κR is large. This figure also shows that the approximation fails when κR is small, that is when the diffusion approximation is not valid, as expected. 5. Numerical experiments. For the numerical experiments, we apply the two versions of the IDSA designed before to the homogeneous sphere test case. We refer to Eq. (78) as the Old IDSA approximation and to Eq. (85) as the New IDSA approximation. We use a backward Euler scheme with a finite difference scheme in space with a very fine grid, N r = 10 6 . Such a fine grid eliminates discretization errors and the remaining discrepancy can be attributed to modeling error. We use 0 initial conditions and boundary conditions are discussed in [7,Chap.4]. The test case we choose is an homogeneous sphere of radius R = 6. The value of κ is varied in the different tests.
As a first example, we apply the Old IDSA and the New IDSA to the homogeneous sphere of radius R = 6 with κ = 1. Figure 6 shows the relative approximation error of these two schemes. We notice that the Old IDSA (in blue) has a considerably larger error in the free streaming region where r ≥ R than the New IDSA (in red). These results show that even if the flux factors are accurate and the diffusion relation given in Proposition 5 is satisfied by the Old IDSA approximation, it does not mean that the approximation is good. The modifications given is Section 4.2 greatly improve the situation. The error still has a peak at the homogeneous sphere radius, but this is expected as a consequence of the definition of approximations (78) and (85). A better view of this example is given in Figure 7 where we also give the value of the flux factors h and k which can be computed from the value of J t and J s from Proposition 5.
The results show that the Old IDSA overestimates the value of streaming particles at the radius of the homogeneous sphere; this leads to an overestimation of the streaming particles in the whole region where r ≥ R, explaining the results of Figure  6. We also notice that, compared to the New IDSA, the Old IDSA underestimates the proportion of the trapped particles and overestimates the proportion of the streaming particles in the diffusive region (r < R). In a CCSN experiment, the overestimation of streaming particles can lead to larger heating rates behind the shock and an artificially increased likeliness for explosions. From this viewpoint, one would say that the Old IDSA overestimates the neutrino heating, see [6]. Regarding the flux factors, we see in Figure 7 that the flux ratio of the Old IDSA is worse that the flux ratio of the New IDSA; that is the correction done in the J distribution reflects in the flux ratio. There are, however, no differences in the approximation of the variable Eddington factors. This is a consequence of the reconstruction of H idsa and K idsa given in Proposition 5 from which the flux factors are computed.
To complement the study, we let the values of κ vary from 1 to 100 and compute the relative L 2 errors. In Figure 8 convergence of these three quantities is O(κ −1/2 ). This rate of convergence has only been obtained numerically and we want to prove this fact in future work.
The first panel of Figure 8 shows that the convergence of J accelerates for large values of κ, but this behavior has yet to be understood. 6. Conclusion and outlook. In this paper, we have discussed the behavior of the IDSA on the homogeneous sphere test case. We have shown that the current implementation of this approximation has modeling issues coming from the definition of the coupling mechanism, that is from the diffusion source Σ. The two main problems are the Spurious Trapped and Instability problems. In the case of the homogeneous sphere, we have changed the coupling mechanism assuming the different regions where each regime apply to be known. Doing so, the two modeling problems are avoided allowing us to access the modeling error of the IDSA. We show with numerical experiments that the IDSA is overestimating the streaming component both in the diffusive and in the free streaming part of the domain. To address the other issue, we proposed a new version of the IDSA based on the analytical solution of the homogeneous sphere. This very particular model has been shown to converge to the solution of the infinitely opaque homogeneous sphere case. We also used the analytical solution of this case to obtain the flux ratio and the variable Eddington factor for the stationary state free streaming limit used in the IDSA. The new version of the IDSA presented has also been shown to converge for the first and second order moments with these flux factors. The order of convergence has been numerically estimated to O(κ −1/2 ) and the proof of this fact is still missing. We have shown in this paper that simplifying the model can shed some light on the behavior of the approximation. To go further in the analysis, we need to go back to more complicated radiative transfer problems. The analysis presented in this work is very specific to the homogeneous sphere example. However, many astrophysical problems do have a very similar structure, a dense body surrounded by a transparent atmosphere. This is the case for example in CCSN simulations, stellar winds, planetary disks and many other astrophysical problems. A possible direction in which our results can be extended is to use the neutrinosphere radius as the limiting radius for trapped particles. Note that Remark 6 shows that these two radius are usually not equal, but close when κ is large. This allows the computation of trapped particles with the New IDSA scheme. Then, one needs to find a way to reconstruct the streaming particle distribution inside the neutrinosphere in order to obtain a boundary condition for the transport of streaming neutrinos outside the neutrinosphere. Once such a condition is obtained, one can use any radiative transfer code to compute the outside distribution.