New boundary conditions for the approximate flux-limited diffusion radiative transfer in circumstellar environments: test case study for spherically symmetric envelopes

In order to constrain the models describing circumstellar environments, it is necessary to solve the radiative transfer equation in the presence of absorption and scattering, coupled with the equation for radiative equilibrium. However, solving this problem requires much CPU time, which makes the use of automatic minimisation procedures to characterise these environments challenging. The use of approximate methods is then of primary interest. One promising candidate method is the flux-limited diffusion (FLD), which recasts the radiative transfer problem into a non-linear diffusion equation. One important aspect for the accuracy of the method lies in the implementation of appropriate boundary conditions (BCs). We present new BCs for the FLD approximation in circumstellar environments that we apply here to spherically symmetric envelopes. At the inner boundary, the entering flux (coming from the star and from the envelope itself) may be written in the FLD formalism and provides us with an adequate BC. At the free outer boundary, we used the FLD formalism to constrain the ratio of the mean radiation intensity over the emerging flux. In both cases we derived non-linear mixed BCs relating the surface values of the mean specific intensity and its gradient. We implemented these conditions and compared the results with previous benchmarks and the results of a Monte Carlo radiative transfer code. A comparison with results derived from BCs that were previously proposed in other contexts is presented as well. For all the tested cases, the average relative difference with the benchmark results is below 2\% for the temperature profile and below 6\% for the corresponding spectral energy distribution or the emerging flux. We point out that the FLD method together with the new outer BC also allows us to derive an approximation for the emerging flux.


Introduction
The study of circumstellar environments at different stages of stellar evolution is of crucial importance. These environments reflect the physical processes in action, from the star formation with the presence of accretion discs to late stages in the evolution, in which strong stellar winds shape the circumstellar envelopes. Observations at high angular resolution allow to probe and characterise the circumstellar material by determining densities, temperatures, abundances, velocity fields, etc. The exploitation of instruments such as the Multi AperTure mid-Infrared SpectroScopic Experiment 1 (MATISSE), operating in the mid-IR, or the Atacama Large Millimeter Array 2 (ALMA) in the sub-millimetric domain offer complementary views of these environments, that give access to regions close to the central star up to the outer regions using a multi-wavelength approach.
Circumstellar matter is generally composed of a mixture of gas and dust particles that absorbs and scatters the incident stellar radiation. The envelope is then heated by the radiation, and a radiative equilibrium can be reached in which the envelope also emits radiation in the infrared domain.
In order to constrain the models describing circumstellar environments, it is necessary to solve the radiative transfer equation under the assumption of radiative equilibrium. Several numerical techniques exist to solve this problem in one, two, and three dimensions. The Monte Carlo method is popular because it can be adapted to any geometry and can handle many physical processes (see Steinacker et al. 2013 for a thorough review). However, solving the radiative transfer problem requires much CPU time, which makes the use of any automatic minimisation procedure to characterise these environments challenging.
In this context, the use of approximate methods is of primary interest. One promising candidate is the flux-limited diffusion (FLD), introduced by Levermore & Pomraning (1981) (L&P hereafter). This description numerically simplifies the problem by recasting the radiative transfer equation into a non-linear diffusion equation for the mean specific intensity of the radiation field (see Sect. 2).
Physically, the boundary conditions (BCs) for the radiative transfer equation are obtained from the known specific intensity incident on the surface of the object. However, the outgoing in-Article number, page 1 of 10 arXiv:2107.13993v1 [astro-ph.SR] 29 Jul 2021 A&A proofs: manuscript no. main tensity is not known a priori and has to be obtained from the solution of the radiative transfer problem. It is thus not obvious to find a BC for the mean specific intensity at the surface of the object. A consequent theoretical work has already been done in finding satisfying BCs for the FLD method (Pomraning 1986(Pomraning , 1988. These BCs were derived with the assumption that a boundary layer could be defined, which might not be true in astrophysical applications where the media may not be seen as an infinite half-space by the radiation for some frequencies. Furthermore, as far as we know, they have never been numerically tested in an astrophysical context. The FLD approximation has already been implemented in several astrophysical applications. Sonnhalter et al. (1995); Yorke & Sonnhalter (2002) used the FLD to solve the frequencydependent radiative transfer to model protostellar discs and massive star formation, respectively. In these studies, the central star was treated as an additional source and the specific mean intensity J ν at the outer edge of the media was set to be equal to the Planck function B ν (T out ), with a prescribed temperature T out at the boundary. Some improvements were made later, in the context of the radiation hydrodynamics problem for massive star formation (see Kuiper et al. 2010;Mignon-Risse et al. 2020). These more sophisticated hybrid codes split the radiation field into two components, the stellar and the dust component, where the FLD method only solves for the latest part. In this treatment, the Dirichlet boundary condition at the outer edge only applies to the dust component. This relies on the assumption that the dust temperature is known at the interface with the interstellar medium. In the problem we consider, the temperature at the outer boundary is not known a priori and must be derived as part of the solution to the radiative transfer problem coupled with the radiative equilibrium equation. In a non-grey problem, we need BCs that can properly handle several regimes of optical thicknesses, for different frequencies.
We present new BCs for the FLD theory in circumstellar environments that we tested and implemented in the case of spherically symmetric envelopes. The condition is derived from the prescription of (i) the incident flux, derived from an extended stellar source and the self-heating of the envelope at the inner boundary and from (ii) the ratio of the mean specific intensity over the radiative flux at the free outer boundary. We show that they both may be written as mixed BCs relating the mean specific intensity and its gradient at the surfaces of the envelope. They consequently lead to a more realistic description of the radiation field (compared to simple Dirichlet or von Neumann boundary conditions) while still remaining sufficiently easy to implement. As a by-product of our investigations, we also derived an approximate expression for the emergent flux at the free outer surface.
The paper is organised as follows: in Sect. 2 we recall the bases of the FLD theory. In Sect. 3 we present the new boundary conditions and in Sect. 4 their numerical implementations. In Sect. 5 we test the accuracy of our results by comparing the temperature profile in the envelope and the spectral energy distribution (SED) of the outgoing flux with the results of two radiative transfer codes, namely (i) DUSTY (Ivezic & Elitzur 1997), which numerically solves the integral equation for the energy density, and (ii) a Monte Carlo (MC) radiative transfer code (Niccolini & Alcolea 2006). Additionally, we compare the derived boundary conditions with the original boundary conditions of Levermore & Pomraning (1981). Finally, in Sect. 6, we conclude and present some perspective for our future work.

The Flux Limited Diffusion Theory
In the following, we present the original work of L&P and introduce the relevant background for the derivation of the BCs in Sect. 3. The position is denoted by r, the direction of propagation byn, and the frequency by the subscript ν. The transport equation for the specific intensity I ν (r,n, t) at the position r in then direction with isotropic and coherent scattering is (1) is the Planck function and κ ext ν , κ abs ν and κ sca ν are the extinction, absorption and scattering coefficients, respectively. The zeroth and first moments of the specific intensity, namely J ν and H ν , are defined as where the integration is performed over all directions. These quantities are linked by the zeroth moment of Eq. (1), The FLD approximation is a closure of the system of the moment equations by expressing H ν as a function of J ν . This is done by expressing the specific intensity I ν as a function of the mean specific intensity J ν , where h ν is the normalised flux and is expressed as The ψ ν function is called the normalised intensity and quantifies the anisotropy of the radiation field I ν . In the optically thin and thick limits, this function reduces to a Dirac distribution and a constant, respectively. The FLD approximation consists of assuming that the anisotropy of the radiation field is a conserved quantity, yielding the expression for ψ ν , with The quantity denoted by R ν plays a key role in the description of the local radiation field in the medium. It expresses the ratio of the (effective) mean free path over the characteristic length of the variation of the mean specific intensity. Consequently, the limits R ν 1 and R ν 1 correspond to the optical thin and thick regimes, respectively. The quantity ω ν is called the effective albedo. It is equal to unity in the absence of true absorption (κ abs ν = 0). L&P showed that, using Eq. (6) in Eq. (5), h ν is proportional to R ν , Article number, page 2 of 10 where R ν is the norm of R ν and λ (R ν ) is the 'flux-limiting' parameter. Finally, using Eq. (8) in Eq. (4) allows Eq.
(3) to be rewritten as a non-linear diffusion equation for the mean specific intensity with the non-linear diffusion coefficient We note that in the FLD approach, the radiative net flux H ν is related to the gradient of the mean specific intensity J ν by H ν = −D ν ∇J ν . It shows some similarities with the Fick law, which applies in the stellar interior, but here the diffusion coefficient depends on a non-linear way on the mean specific intensity. In the optically thick regime, the FLD approximation reduces to a linear diffusion equation, whereas in the optically thin regime, it reduces to an advection equation, as expected.

Boundary conditions
In this section, we specify the time-independent BCs that were implemented for the FLD Eq. (9). They are defined at the specific boundaries of a circumstellar shell, namely an inner spherical cavity that is illuminated by an enclosed star, and an outer boundary with no incoming radiation. The BCs for the radiative transfer equation Eq. (1) would be obtained by specifying the value of the incident specific intensity I ν on the considered surface. In the framework of the FLD approximation, we face two problems: (i) The FLD equation applies to the mean specific intensity, whereas the physical constraint is on the ingoing specific intensity at the boundaries, and (ii) we also explained in Sect. 2 that the FLD approach implies a specific angular dependence of I ν , given by the function ψ ν (see Eq. (6)). This specific dependence would not be consistent with any arbitrary value of the incident radiation field at the surface. Consequently, the actual BCs for the radiative transfer equation are in general incompatible with the FLD solution. This is expected because the FLD approximation is in principle not valid close to the boundaries of objects. Pomraning (1988) derived the BCs for the FLD equation from the decomposition of the radiative transfer problem into an interior problem described by the FLD equation and an additional boundary layer term. The match of the interior and boundary layer solutions yields the BCs for the FLD equation. A relation of the surface values of the mean specific intensity and its gradient is obtained, but with rather involved coefficients depending on integrals of the Chandrasekhar H function and on the surface value of R ν . In the literature, Dirichlet BCs are often used with a prescribed value for the temperature at the outer boundary of the medium. As already pointed out, the surface temperature is not known a priori and must be derived from the solution of the radiative transfer problem.
Here, we propose to impose BCs for the zeroth and the first angular moments of the radiation field, in a form that is compatible with the FLD approach, in order to ensure a smooth match with the interior FLD solution. We obtained mixed Robin-type BCs that relate the surface values of the mean specific intensity and its gradient, but the coefficients are quite simple analytical functions of the surface value of R ν .

Expression of the incident flux in the FLD formalism
One example of physical significance is to write a condition on the flux F in ν entering a boundary surface at the position r S , .n ψ ν (r S ,n) d 2n .
Here,ŝ is the outward normal vector to the surface of the envelope (ŝ = −r at the inner edge). The right-hand side (RHS) of Eq. (11) is the expression of an incoming flux in the FLD formalism. Eq. (11) is the general form of the inner BC without any assumption on the geometry of the problem. The integral can be performed analytically using Eq. (6) for ψ ν if we assume that the vector R ν is normal to the surface (R ν = ± R νr ), which is exact for spherically symmetric problems. The equation can then be rewritten as a non-linear mixed BC, This BC may be regarded as an implicit relation of the surface values of J ν and R ν that yields an FLD solution compatible with the given incident flux. It is interesting to note that although Eq. (12) is different in its coefficients and construction from the original BC Eqs. (56) and (66) in L&P, they are analytically equivalent. In L&P, the coefficients were derived to give an exact transport result for the case of a source-free (B ν = 0) half-space media, with a constant κ ext ν , ω ν and a particular incident intensity distribution Γ(r S , µ) of the form This specific form is proportional to the angular dependence in the FLD formalism. The correspondence between α ν in Eq. (13) and Eq. (56) in L&P is Hence, specifying the FLD flux at the boundaries of a spherically symmetric domain will actually give the same BC as is obtained by solving the exact transport result of a source-free (B ν = 0) half-space media, with a constant κ ext ν , ω ν and a particular incident intensity distribution Γ(r S , µ), given by Eq. (14). We note that this equivalence only applies for the spherical and planar symmetric systems and no longer holds when we specify Eq. (11) for other configurations, where the radiative flux is not orthogonal to the boundary surface.

Incident flux from an extended source and the envelope
We want to write the BC at an inner spherical cavity located at a distance R in from the centre of a star of radius R (r S = R inr ) and surface temperature T . For this, we need to specify the flux F in ν . We have two contributions, Article number, page 3 of 10 A&A proofs: manuscript no. main The first contribution, denoted by B ν (T ) comes from the star and the other from the inner boundary itself and is expressed as I e ν (r, −n) = J ν (r) ψ ν (r, −n). As shown by Fig.1, the vector r = r S − 2(r S .n)n corresponds to the opposite point at the inner boundary, alongn. Because of this dependence, this BC is no longer local by nature and cannot be expressed in a closed form, except in spherical symmetry where J ν (r) = J ν (r S ) and ψ ν (r, −n) = ψ ν (r S , −n). To perform the angular integration, we aligned the n z axis with the unitary vectorr. For the star, the integration on µ = cos(θ) (θ being the angle between n z andn) spans from µ 0 = 1 − (R /R in ) 2 to 1, and for the inner cavity, it spans from 0 to µ 0 , The incident flux F in ν is then expressed as, with, The inner BC for the FLD equation in the case of an inner spherical cavity enclosing a star is with α ν and β ν , the quantities defined by Eq. (13). We can analytically write the FLD solution with this inner boundary condition in the limit where the envelope is optically thin (R ν 1) and compare it with the known analytical solution for the dilution of the mean specific intensity in free space, In the optically thin limit, the FLD Eq. (9) and the BC Eq. (20) reduce to respectively. The solution of the FLD equation, in this limit is then in agreement with Eq. (21) when (R /r) 2 1. The relative difference in temperature between Eq. (21) and the FLD solution Eq. (23) is ≈ 15 % at the star surface, 1 % at r ≈ 2.5 R . For r > ∼ 10 R (values for our test cases presented in Sect. 5), this difference becomes negligible ( < ∼ 0.06 %).

Outer boundary
For the inner cavity, we would like to impose that no incoming radiation enters the external shell of the envelope (F in ν = 0). However, in the FLD formalism, the angular dependence of the radiation is given by the specific form of the ψ ν function (Eq. (6)). The incoming radiation vanishes only when R ν becomes infinite, which also results in a sharp-peaked distribution for the emerging specific intensity. This is not physically realistic, and this inconsistency is expected because the FLD method is rigorously not valid close to the surface of the object. Another approach is then required to describe the behaviour of the radiation on the external edge. Inspired by L&P and Pomraning (1986), we seek a BC in the form of a closure relation between the mean specific intensity and the radiation flux at the surface, that is, with ζ ν , a coefficient we need to determine. At the outer edge and without incoming radiation, ζ ν has to be understood as the ratio of the energy density over the emerging flux. This ratio can be expressed in spherical symmetry as It depends on the anisotropy of the emergent radiation field. In the optically thin limit, the radiation field is along theŝ direction (spherical symmetry), thus I ν ∝ δ(µ − 1) and hence ζ ν − → 1. In the diffusion regime where the emergent field is isotropic, ζ ν − → 2. In the original study of L&P, ζ ν was chosen to be equal to 2, which means that it correctly describes the optically thick cases where R ν 1. If the BC is to correctly describe different optical regimes, we need ζ ν to be a function of R ν . In the framework of the FLD approximation, I ν = J ν ψ ν and the specific angular dependence of I ν , given by ψ ν , is used to compute the surface value of ζ ν , In the two limits (R ν 1 and R ν 1), we recover ζ ν − → 1 and ζ ν − → 2. We note that the latter limit reduces to the boundary condition Eq. (56) in L&P. Now we need to specify the value of R ν at the external boundary of the envelope. The behaviour of the radiation at this interface is only dictated by the interior solution Article number, page 4 of 10 because there is no incoming radiation. Consequently, to ensure a smooth match of the BC and the interior FLD solution, we used a second-order extrapolated value for R ν (R out ) = R e ν (see Eq. 41). Because we prescribe a value for R ν at the external edge, R e ν is also used to compute the non-linear diffusion coefficient D ν (R out ) = D ν (R e ν ) = λ(R e ν )/ω ν κ ext ν and the coefficients of the numerical scheme Eq.(34) on the external edge. The outer BC without incoming radiation that we implemented is then where we have expressed H ν (R out ) = −D ν (R e ν ) ∇ J ν | r=R out . In Sect. 5.2 we compare this BC in an astrophysical application with respect to the original BCs of L&P.

Approximation for the emergent flux
We used an extrapolation of the non-linear diffusion coefficient D ν (R e ν ) to relate the flux at the external edge to the gradient of the mean specific intensity. As there is no incident flux, this provides us with an approximate expression for the emergent flux, given by This approximation is tested in Sect. 5. It reproduces the results of different radiative transfer codes that solve the full transfer equation very well, as shown in Fig. 2 and Fig. 3.

Radiative equilibrium and warming of the stellar surface
Because of the geometric extension of the star (see Fig. 1), there is part of the radiation that emerges from the envelope that falls back onto the star, with R ν being evaluated at r = R in . We used the same conventions as in Sect. 3.1 to perform the angle integration. This part of the flux is hidden from the rest of the envelope, and the radiative equilibrium inside the cavity leads to a warming of the stellar surface (Niccolini et al. 2003). This effect can be quite dramatic, and reach up to 30 % of the total stellar luminosity that is obscured for optically-thick grey shell, as in the test case presented in Sect. 5.2. This has to be taken into account to properly ensure the radiative equilibrium condition throughout the full space from the stellar surface to the outer boundary of the envelope. To fulfil the radiative equilibrium condition at the stellar surface we write We imposed a fixed value for T eff , and the temperature of the star T was updated accordingly.

Numerical implementation
The FLD Eq. (9) is a non-linear diffusion equation that has to be solved numerically for each point of space, time, and frequency.
In the following, we limit ourselves to the 1D time-independent FLD equation, Here, r denotes the radial variable from the centre of the envelope. An additional constraint to Eq. (31) is given by the equation of the radiative equilibrium These equations are solved for a spherically symmetric envelope of inner radius R in and outer radius R out , surrounding a star of radius R .

Finite-difference approximation
Following a finite-difference procedure, we discretise and sample in a logarithmic way the frequency domain into n ν points, denoted by the subscript k. Space is discretised into n x cells and denoted by the subscript i. J ν is computed at the cell centres and the vector H ν on the walls. The differential operator is approximated with a second-order finite differences operator. The equation is solved with respect to a new variable x = f (r) on a regular grid of constant step ∆x = (x(R out ) − x(R in )) / (n x − 2) to allow r to be non-uniformly sampled. We make use of one ghost cell for each grid border to ensure the BCs. We obtain the following system of equations: The non-linear coefficients A are given by The non-linear nature of the equation arises from the expression of the coefficients A and RHS in Eq. (33). They implicitly depend on J ν , through the diffusion coefficient D ν and the radiative equilibrium Eq. (32), respectively. The coefficients A require an estimation of J ν and its gradient (see Eqs. (10) and (7)) at the cell walls, given by

Iterative scheme
Several strategies are possible in order to solve the FLD Eq. (33) coupled with Eq. (32). The simplest approach is to use an iterative method to fully solve Eq. (33) and to update the temperature through Eq. (32). Iterating between these two processes until convergence yields the solution of the problem. This procedure, commonly called the Λ-iteration in the literature (see Yorke & Sonnhalter 2002), becomes very slow and does not converge for large optical depths. In analogy with the usual accelerated Λ-iteration (ALI) methods, Yorke & Sonnhalter (2002) found an improved convergence behaviour by splitting the solution of equation Eq. (31) in Eq. (32).
We found a simple method, inspired by the Gauss-Seidel approach, to solve Eqs. (33) and (32) simultaneously instead of repetitively. If we denote by n the iteration of the method, we have and the temperature is updated after only one Gauss-Seidel spatial sweep at each frequency through Eq. (32), which we rewrite as We replaced the frequency integration by a quadrature formula with the associated weights W k . The left-hand side (LHS) of Eq. (37) was pre-computed and stored in a table, for a wide range of temperatures, allowing the RHS to be linearly interpolated in this table (in the logarithm of the integral for better accuracy) . By doing so, we avoided using a Newton-Raphson procedure to determine the new temperature, which reduces the computational time.
Our procedure consisted of repeatedly updating J ν and T with the help of Eqs. (36) and (37) until we reached convergence. The coefficients A Eq. (34) and the BCs Eq. (40) were immediately updated for each frequency k after one Gauss-Seidel spatial sweep. We note that this procedure is different from the usual Λ-iteration presented above because the temperature is updated simultaneously with J ν , within the same iteration n.

Update of the stellar temperature
The radiative equilibrium inside the inner cavity requires updating the stellar temperature (see Sect. 3.4). Following Eq. (30), we updated the stellar temperature at the end of each iteration n of the numerical scheme, Here again, we replaced the frequency integration by a quadrature formula with the associated weights W k . F fall,n+1 k is given by Eq. (29)
Accordingly, the values of J ν in the ghosts cells were updated immediately after one Gauss-Seidel sweep Eq. (36) to ensure the BCs, with α k , β k and γ k being defined by Eqs. (13) and (19). For the extrapolated value R e k in ζ(R e k ) Eq. (26) and D k (R e k ), we used a second-order Lagrange extrapolation,

Initial conditions
Initial conditions for both J ν and T must be provided in order to solve Eq. (33). It is clear that the overall convergence speed strongly depends on the initial setup of the solution, but there is also a trade-off with the stability, that is, the ability of the solution to converge, for a wide variety of cases. As an initial guess, we used the analytic solution of the FLD in the optically thin limit Eq. (23), and we write from which we deduce the corresponding temperature profile T 0 with the help of Eq. (37).

Numerical tests: spherically symmetric envelopes
5.1. Benchmarks from  We tested the accuracy of our FLD code with our Robin-type mixed boundary conditions in a general and realistic case, by comparing it with the 1D benchmark problems realised by . We recall the conditions of the test and refer to the original paper for further information. A point source surrounded by a spherically symmetric envelope of matter at radiative equilibrium irradiates as a black body at the temperature T = 2 500 K. This envelope extends from the inner radius R in to the outer radius R out = 1 000 R in . The inner radius is set so that the temperature at the inner radius is always T in = T (R in ) = 800 K. The density profile n(r) is assumed to be  Fig. 2. Non-grey case : Normalised temperature profiles (upper panels) and SEDs (lower panels) for four different opacities τ ν 0 = 1, 10, 100, and 1 000 (blue, orange, green and red, respectively) and two density power laws: p = 0 (left panels) and p = 2 (right panels). The solid lines represent the FLD curves, and the black dots indicate the benchmark profiles from  a power law of the form n(r) = n 0 (R in /r) p . The radial optical depth τ ν of the envelope is linked to the density profile by where C ext ν is the extinction cross-section coefficient. It is defined by with C abs ν 0 = (1 − η) C ext ν 0 , C sca ν 0 = η C ext ν 0 , ν 0 the frequency corresponding to λ = 1 µm and η the albedo, set to 1/2 for these tests. The benchmark problems are thus completely defined by two parameters: (i) the exponent in the density power law p = 0, 2, and (ii) the radial optical depth of the envelope at ν 0 , τ ν 0 = 1, 10, 100, and 1 000. This created eight different cases to test the accuracy of our code. The coefficient n 0 in the density profile is derived with the help of τ ν 0 and p, The normalised temperature profile T/T in and the normalised SED λF λ /F (F = ∞ 0 F λ dλ) of the envelope are shown in Fig. 2 for each case.
The Ivezic benchmarks were produced with version 2 of DUSTY 3 (Ivezic & Elitzur 1997). Because our code is of different nature, the spatial and frequency grids are different. We then compared the results by linearly interpolating our profiles (in log − log scale) on the DUSTY grids. We used 128 points for space and frequency, with a logarithmic sampling. The corresponding relative differences are displayed in Table 1. We also point out that we restricted the comparison, for the normalised SEDs, to the frequency domain where λF λ /F ≥ 10 −6 because of non-physical results of the DUSTY code below this threshold, for the smallest wavelengths.
The FLD results and the benchmarks agree well. The average relative differences in the temperature profiles mean( (T )) is of the order of 1 %, with a maximum value of approximately 4 %, achieved by the most optically thick envelopes (τ ν 0 = 1 000). The average of the relatives differences in the normalised SEDs mean( (λF λ /F)) always stays below 3 %, with the exception of the optically thick envelope with constant density profile, where this difference reaches 6 %. Niccolini & Alcolea (2006) We compared the FLD code with a 3D MC radiative transfer code (Niccolini & Alcolea 2006 Levermore & Pomraning (1981) (dashed lines) for the test case presented in Sect. 5.2. The colour code is the same as in Fig.3. The emerging fluxes are displayed in semi-log scale to highlight the differences between the BCs. The lower panels show the relative differences profiles with respect to the MC code from Niccolini & Alcolea (2006) in temperature (left) and in emerging flux (right).

Grey spherical shell with the Monte Carlo code from
in the less realistic but more extreme case of a spatially small spherically symmetric grey envelope. We expect the boundary effects to play a major role for this type of problems. The inner radius was set to R in = 10 R and the outer radius R out = 20 R . We assumed a constant density profile (p = 0) in the envelope. Our test cases consisted of determining the normalised tempera-ture profiles T/T in and the emerging fluxes F λ for several cases, ranging from optically thin (τ = 0.01) up to the optically thick envelopes (τ = 100). The corresponding profiles are shown in Fig. 3.
Because the codes are different, we interpolated our results linearly (in log − log scale) on the MC grids. The relative differ-  (T ) (F λ ) 0.01 0/0/0 0/3/17 1 1/1/5 0/3/17 10 1/1/7 1/8/33 100 0/0/3 1/2/13 Notes. Relative differences (in %) for the temperature profile (T ) and for emerging flux (F λ ) shown in Fig. 3. The differences are presented in the form mean( ) / std( ) / max( ) and rounded to the closest percent. mean( ) is given by Eq. (46). ences between the two codes, are displayed in Table 2. As an additional feature, the MC code also provides an estimation of the errors on the temperature σ(T ) and on the emerging flux σ(F λ ), computed from the MC noise (Niccolini & Alcolea 2006). We used this information to compute a more relevant mean value for (T ) and (F λ ), where N x (N λ ) is the number of spatial (wavelength) points of the MC grid, (T i ) ( (F k )) is the relative error (in %) on the temperature (emerging flux) between our results and the MC results, and W i (W k ) is the inverse square of the MC relative errors, defined as The two results agree well. The average of the relative differences of the temperature profile mean( (T )), remains of the order of 1 % for all the cases we tested. The largest differences are reached for the intermediate cases (τ = 1, 10) and are located on the external edge of the envelope. This is expected because the FLD approximation is known to perform well in the optically thin and thick regimes, but it is less well suited to describe these intermediate cases. Nevertheless, the temperature profile is still quite well reproduced and the emerging flux is not affected by the small errors on the temperature close to the outer edge. We point out that the BCs derived in this paper allow us to successfully reproduce the correct behaviour of the temperature profile for the optically thick envelope where it shows a quite steep decrease at the outer surface, as shown by Fig. 3. The emerging fluxes F λ agree within 1 % on average, except for the optically thick case (τ = 100), where mean( (F λ )) reaches about 5 %. We conclude this section with a comparison of the outer BC from this study Eq. (27) and the original BC Eq. (56) from L&P, given by, We previously mentioned in Sect. 3.1 that the inner BC Eq. (20) is analytically identical to that of L&P in spherical symmetry, so we restrict the comparison to the outer edge of the envelope. It is important to notice that the BC of L&P was not originally intended to describe this class of problems, but the comparison still remains instructive for studying the importance of the BCs for the accuracy of the solution. In Eq. (48), the factor 2 (Eq. 58 in L&P) was originally used to give the correct ratio of the energy density over the emerging radiative flux in plane-parallel geometry for an optically thin slab illuminated by an isotropic incident radiation field. However, for the special case of a spherical envelope surrounding a black-body star, this ratio becomes , the cosine of the stellar angular size at R out . As µ 0 → 0 or equivalently R /R out → 1, this ratio increases to 2, as expressed by L&P, because we recover the case of a plane-parallel geometry with an isotropic incident radiation field. Far from the star (µ 0 → 1 or R /R out → 0), the ratio tends to 1, associated with an incoming sharp-peaked radiation. Hence, the L&P BC that set the ratio to 2 will strongly deviate from the analytic limit 1, in the optically thin regime. In Fig. 4 we display the relative differences in the temperature profiles (T ) and in the emerging fluxes (F λ ) for the same test case as presented at the beginning of this section. We note that this test case is not realistic, however, it allows to compare different optical regimes, in contrast to a more realistic problem in which the radiation is free in the external regions most of the time, such as for the test case presented in Sect. 5.1. We recall that the BC from L&P is a limiting case of the BC derived in this study, in the case R 1, where the emerging radiation is almost isotropic (see Sect. 3.2), hence it is not surprising that the results converge to the same profile, for a optically thick grey envelope (τ = 100). On the other-hand, in the optically thin case (τ = 0.01) where we would expect the ζ coefficient Eq. (26) to be close to unity, the BC from L&P performs poorly as expected. For intermediates regimes (τ = 1, 10), the entanglement of the error of the BC and FLD method itself makes any comparison very hard. We note that although the temperature profile is closer to the benchmark result for τ = 10 with the BC of L&P, this is not the case for the associated emerging flux. We also note that the BC, although defined locally, can have a global effect on the whole solution, as shown by the temperature profile of the test case τ = 1. To conclude, we also point out that we tried to implement the L&P BC in the non-grey cases (Sect. 5.1) but we were unable to reach a satisfying convergence of the computations.

Conclusion
The FLD approximation together with the new BCs, yields promising results in correctly describing the radiation transport inside spherically symmetric circumstellar envelopes. These conditions, derived in Sect. 3 from physically consistent constraints on the behaviour of the radiation field at the inner and outer surfaces, allow us to compute with a good accuracy the temperature profile and the SED for a wide range of configurations, from very small to very large optical thicknesses. As shown in Sect. 5, it reproduces the correct temperature profile within ≤ 2%, and the SED or emerging flux at less than ≤ 6% on average, with respect to the solution of the full radiative transfer equation under radiative equilibrium.
The numerical solution of the 1D non-linear diffusion equation Eq. (31) coupled with the radiative equilibrium Eq. (32) was performed with a Gauss-Seidel method-based iterative scheme, in which the temperature is updated at each iteration step. Furthermore, we point out that the FLD approximation implemented with the proposed outer-boundary condition provides a simple approximation for the flux emitted by the envelope. This allows computing the SED or emerging flux without using any