© J. Perdigon et al. 2021

1. 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¹ (MATISSE), operating in the mid-IR, or the Atacama Large Millimeter Array² (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 intensity 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, 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 frequency-dependent 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.

2. 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 by $\widehat{n}$, and the frequency by the subscript ν. The transport equation for the specific intensity $I_{\nu }(\mathit{r},\widehat{n},t)$ at the position r in the $\widehat{n}$ direction with isotropic and coherent scattering is

$$\begin{array}{c}\hfill \frac{1}{c}{\partial }_t{I}_{\nu }+\widehat{n}.\mathbf{\nabla }{I}_{\nu }=-{\kappa }_{\nu }^{\mathrm{ext}}{I}_{\nu }+{\kappa }_{\nu }^{\mathrm{abs}}{B}_{\nu }+{\kappa }_{\nu }^{\mathrm{sca}}{J}_{\nu }.\end{array}$$(1)

J_ν = J_ν(r,t) is the mean specific intensity, B_ν = B_ν(T(r,t)) is the Planck function and ${\kappa }_{\nu }^{\text{ext}}$, ${\kappa }_{\nu }^{\text{abs}}$ and ${\kappa }_{\nu }^{\text{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

$$\begin{array}{c}\hfill J_{\nu }=\frac{1}{4\pi }{\displaystyle \underset{4\pi }{\int }{I}_{\nu }{\mathrm{d}}^2\widehat{n},{\mathit{H}}_{\nu }=\frac{1}{4\pi }\underset{4\pi }{\int }{I}_{\nu }\widehat{n}{\mathrm{d}}^2\widehat{n},}\end{array}$$(2)

where the integration is performed over all directions. These quantities are linked by the zeroth moment of Eq. (1),

$$\begin{array}{c}\hfill \frac{1}{c}{\partial }_t{J}_{\nu }+\mathbf{\nabla }.{\mathit{H}}_{\nu }={\kappa }_{\nu }^{\mathrm{abs}}(B_{\nu }-J_{\nu }).\end{array}$$(3)

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_ν,

$$\begin{array}{c}\hfill I_{\nu }=J_{\nu }{\psi }_{\nu }(\mathit{r},\widehat{n},t),{\mathit{H}}_{\nu }=J_{\nu }{\mathit{h}}_{\nu }(\mathit{r},t),\frac{1}{4\pi }{\displaystyle \underset{4\pi }{\int }{\psi }_{\nu }{\mathrm{d}}^2\widehat{n}=1,}\end{array}$$(4)

where h_ν is the normalised flux and is expressed as

$$\begin{array}{c}\hfill {\mathit{h}}_{\nu }(\mathit{r},t)=\frac{1}{4\pi }{\displaystyle \underset{4\pi }{\int }\widehat{n}{\psi }_{\nu }(\mathit{r},\widehat{n},t){\mathrm{d}}^2\widehat{n}.}\end{array}$$(5)

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 ψ_ν,

$$\begin{array}{c}\hfill {\psi }_{\nu }(\mathit{r},\widehat{n},t)=\frac{1}{1+({\mathit{h}}_{\nu }-\widehat{n}).{\mathit{R}}_{\nu }},\end{array}$$(6)

with

$$\begin{array}{c}\hfill {\mathit{R}}_{\nu }(\mathit{r},t)=\frac{-\mathbf{\nabla }{J}_{\nu }}{{\omega }_{\nu }{\kappa }_{\nu }^{\mathrm{ext}}{J}_{\nu }},{\omega }_{\nu }=\frac{{\kappa }_{\nu }^{\mathrm{abs}}{B}_{\nu }+{\kappa }_{\nu }^{\mathrm{sca}}{J}_{\nu }}{{\kappa }_{\nu }^{\mathrm{ext}}{J}_{\nu }}·\end{array}$$(7)

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 (${\kappa }_{\nu }^{\text{abs}}$ = 0). L&P showed that, using Eq. (6) in Eq. (5), h_ν is proportional to R_ν,

$$\begin{array}{c}\hfill {\mathit{h}}_{\nu }=\lambda \left(R_{\nu }\right){\mathit{R}}_{\nu },\lambda \left(R_{\nu }\right)=\frac{1}{R_{\nu }}(\frac{1}{tanh{R}_{\nu }}-\frac{1}{R_{\nu }}),\end{array}$$(8)

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

$$\begin{array}{c}\hfill \frac{1}{c}{\partial }_t{J}_{\nu }-\mathbf{\nabla }.\left(D_{\nu }\mathbf{\nabla }{J}_{\nu }\right)={\kappa }_{\nu }^{\mathrm{abs}}(B_{\nu }-J_{\nu }),\end{array}$$(9)

with the non-linear diffusion coefficient

$$\begin{array}{c}\hfill D_{\nu }=\frac{\lambda \left(R_{\nu }\right)}{{\omega }_{\nu }{\kappa }_{\nu }^{\mathrm{ext}}}.\end{array}$$(10)

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.

3. 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_ν.

3.1. Inner spherical cavity

3.1.1. Expression of the incident flux in the FLD formalism

One example of physical significance is to write a condition on the flux $F_{\nu }^{\text{in}}$ entering a boundary surface at the position r_S,

$$\begin{array}{c}\hfill F_{\nu }^{\mathrm{in}}({\mathit{r}}_{\mathrm{S}})=J_{\nu }({\mathit{r}}_{\mathrm{S}}){\displaystyle \underset{\widehat{s}.\widehat{n}\le 0}{\int }\widehat{s}.\widehat{n}{\psi }_{\nu }({\mathit{r}}_{\mathrm{S}},\widehat{n}){\mathrm{d}}^2\widehat{n}.}\end{array}$$(11)

Here, $\widehat{s}$ is the outward normal vector to the surface of the envelope ($\widehat{s}=-\widehat{r}$ at the inner edge). The right-hand side (RHS) of Eq. (11) is the expression of an incoming flux in the FLD formalism. Equation (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 (${\mathit{R}}_{\nu }=\pm R_{\nu }\widehat{r}$), which is exact for spherically symmetric problems. The equation can then be rewritten as a non-linear mixed BC,

$$\begin{array}{c}\hfill F_{\nu }^{\mathrm{in}}({\mathit{r}}_{\mathrm{S}})=\pi ({\alpha }_{\nu }{J}_{\nu }({\mathit{r}}_{\mathrm{S}})+{\beta }_{\nu }\widehat{s}.\mathbf{\nabla }{J}_{\nu }{|}_{\mathit{r}={\mathit{r}}_{\mathrm{S}}})\end{array}$$(12)

with

$$\begin{array}{c}\hfill {\alpha }_{\nu }={2\frac{ln(cosh{R}_{\nu })}{R_{\nu }tanh{R}_{\nu }}|}_{\mathit{r}={\mathit{r}}_{\mathrm{S}}},{\beta }_{\nu }=2D_{\nu }{|}_{\mathit{r}={\mathit{r}}_{\mathrm{S}}}.\end{array}$$(13)

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 ${\kappa }_{\nu }^{\text{ext}}$, ω_ν and a particular incident intensity distribution Γ(r_S, μ) of the form

$$\begin{array}{c}\hfill \mathrm{\Gamma }({\mathit{r}}_{\mathrm{S}},\mu )=\frac{1}{coth(R_{\nu }({\mathit{r}}_{\mathrm{S}}))-\mu }=R_{\nu }({\mathit{r}}_{\mathrm{S}}){\psi }_{\nu }({\mathit{r}}_{\mathrm{S}},\mu ).\end{array}$$(14)

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

$$\begin{array}{c}\hfill \gamma =\frac{{\alpha }_{\nu }-2\lambda (R_{\nu })\widehat{s}.{\mathit{R}}_{\nu }}{2-4\lambda (R_{\nu })\widehat{s}.{\mathit{R}}_{\nu }}.\end{array}$$(15)

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 ${\kappa }_{\nu }^{\text{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.

3.1.2. 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_⋆ (${\mathit{r}}_S=R_{\text{in}}\widehat{r}$) and surface temperature T_⋆. For this, we need to specify the flux $F_{\nu }^{\text{in}}$. We have two contributions,

$$\begin{array}{c}\hfill F_{\nu }^{\mathrm{in}}={\displaystyle \underset{\widehat{s}.\widehat{n}\le 0}{\int }}\widehat{s}.\widehat{n}[B_{\nu }(T_{\star })+I_{\nu }^{\mathrm{e}}(\mathit{r},-\widehat{n})]{\mathrm{d}}^2\widehat{n}.\end{array}$$(16)

The first contribution, denoted by B_ν(T_⋆) comes from the star and the other from the inner boundary itself and is expressed as $I_{\nu }^{\text{e}}(\mathit{r},-\widehat{n})=J_{\nu }\left(\mathit{r}\right){\psi }_{\nu }(\mathit{r},-\widehat{n})$. As shown by Fig. 1, the vector $\mathit{r}={\mathit{r}}_{\text{S}}-2({\mathit{r}}_{\text{S}}.\widehat{n})\widehat{n}$ corresponds to the opposite point at the inner boundary, along $\widehat{n}$. 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 ${\psi }_{\nu }(\mathit{r},-\widehat{n})={\psi }_{\nu }({\mathit{r}}_{\text{S}},-\widehat{n})$. To perform the angular integration, we aligned the n_z axis with the unitary vector $\widehat{r}$. For the star, the integration on μ = cos(θ) (θ being the angle between n_z and $\widehat{n}$) spans from ${\mu }_0=\sqrt{1-{(R_{\star }/R_{\text{in}})}^2}$ to 1, and for the inner cavity, it spans from 0 to μ₀,

$$\begin{array}{c}\hfill F_{\nu }^{\mathrm{in}}=2\pi ({\displaystyle \underset{{\mu }_0}{\overset{1}{\int }}}\mu B_{\nu }(T_{\star })d\mu +J_{\nu }({\mathit{r}}_S){\displaystyle \underset{0}{\overset{{\mu }_0}{\int }}}\mu {\psi }_{\nu }({\mathit{r}}_S,-\mu )\mathrm{d}\mu ).\end{array}$$(17)

Fig. 1. Geometry of an inner spherical cavity (solid black line) illuminated by a central star (grey disc).

The incident flux $F_{\nu }^{\text{in}}$ is then expressed as,

$$\begin{array}{c}\hfill F_{\nu }^{\mathrm{in}}=\pi [{\left(\frac{R_{\star }}{R_{\mathrm{in}}}\right)}^2{B}_{\nu }(T_{\star })+{\gamma }_{\nu }{J}_{\nu }\left({\mathit{r}}_S\right)],\end{array}$$(18)

with,

$$\begin{array}{c}\hfill {\gamma }_{\nu }={\frac{2[{\mu }_{0}tanh{R}_{\nu }-ln(1+{\mu }_{0}tanh{R}_{\nu })]}{R_{\nu }tanh{R}_{\nu }}|}_{\mathit{r}={\mathit{r}}_S}.\end{array}$$(19)

The inner BC for the FLD equation in the case of an inner spherical cavity enclosing a star is

$$\begin{array}{c}\hfill ({\alpha }_{\nu }-{\gamma }_{\nu })J_{\nu }\left(R_{\mathrm{in}}\right)+{\beta }_{\nu }\widehat{s}.\mathbf{\nabla }{J}_{\nu }{|}_{r=R_{\mathit{in}}}={\left(\frac{R_{\star }}{R_{\mathrm{in}}}\right)}^2{B}_{\nu }(T_{\star }),\end{array}$$(20)

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,

$$\begin{array}{c}\hfill J_{\nu }=\frac{B_{\nu }(T_{\star })}{2}(1-\sqrt{1-{\left(\frac{R_{\star }}{r}\right)}^2}).\end{array}$$(21)

In the optically thin limit, the FLD Eq. (9) and the BC Eq. (20) reduce to

$$\begin{array}{c}\hfill r^2{J}_{\nu }=\mathrm{const}.=R_{\mathrm{in}}^2{J}_{\nu }\left(R_{\mathrm{in}}\right),J_{\nu }\left(R_{\mathrm{in}}\right)=\frac{1}{4}{\left(\frac{R_{\star }}{R_{\mathrm{in}}}\right)}^2{B}_{\nu }(T_{\star }),\end{array}$$(22)

respectively. The solution of the FLD equation, in this limit is then

$$\begin{array}{c}\hfill J_{\nu }=\frac{1}{4}{\left(\frac{R_{\star }}{r}\right)}^2{B}_{\nu }(T_{\star }),\end{array}$$(23)

in agreement with Eq. (21) when (R_⋆/r)² ≪ 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%).

3.2. Outer boundary

For the inner cavity, we would like to impose that no incoming radiation enters the external shell of the envelope ($F_{\nu }^{\text{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,

$$\begin{array}{c}\hfill J_{\nu }(R_{\mathrm{out}})-{\zeta }_{\nu }\widehat{s}.{\mathit{H}}_{\mathit{\nu }}(R_{\mathrm{out}})=0,\end{array}$$(24)

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

$$\begin{array}{c}\hfill {\zeta }_{\nu }=\frac{{\int }_0^1{I}_{\nu }\mathrm{d}\mu }{{\int }_0^1\mu I_{\nu }\mathrm{d}\mu }.\end{array}$$(25)

It depends on the anisotropy of the emergent radiation field. In the optically thin limit, the radiation field is along the $\widehat{s}$ direction (spherical symmetry), thus I_ν ∝ δ(μ − 1) and hence ${\zeta }_{\nu }\to 1$. In the diffusion regime where the emergent field is isotropic, ${\zeta }_{\nu }\to 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 ζ_ν,

$$\begin{array}{c}\hfill {\zeta }_{\nu }=\zeta (R_{\nu })=\frac{2+\alpha (R_{\nu })tanh(R_{\nu })}{\alpha (R_{\nu })+2\lambda (R_{\nu })R_{\nu }}.\end{array}$$(26)

In the two limits (R_ν ≫ 1 and R_ν ≪ 1), we recover ${\zeta }_{\nu }\to 1$ and ${\zeta }_{\nu }\to 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 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_{\nu }(R_{\text{out}})=R_{\nu }^{\text{e}}$ (see Eq. (41)). Because we prescribe a value for R_ν at the external edge, $R_{\nu }^{\text{e}}$ is also used to compute the non-linear diffusion coefficient $D_{\nu }(R_{\text{out}})=D_{\nu }(R_{\nu }^{\text{e}})=\lambda (R_{\nu }^{\text{e}})/{\omega }_{\nu }{\kappa }_{\nu }^{\text{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

$$\begin{array}{c}\hfill J_{\nu }(R_{\mathrm{out}})+\zeta (R_{\nu }^{\mathrm{e}})D_{\nu }(R_{\nu }^{\mathrm{e}})\widehat{s}.{\mathbf{\nabla }{J}_{\nu }|}_{r=R_{\mathrm{out}}}=0,\end{array}$$(27)

where we have expressed ${\mathit{h}}_{\nu }(R_{\text{out}})=-D_{\nu }(R_{\nu }^{\text{e}})\nabla {J_{\nu }|}_{r=R_{\text{out}}}$. In Sect. 5.2 we compare this BC in an astrophysical application with respect to the original BCs of L&P.

3.3. Approximation for the emergent flux

We used an extrapolation of the non-linear diffusion coefficient $D_{\nu }(R_{\nu }^{\text{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

$$\begin{array}{c}\hfill {\mathit{F}}_{\nu }^{\mathrm{out}}=4\pi {\mathit{H}}_{\nu }(R_{\mathrm{out}})=-4\pi D_{\nu }(R_{\nu }^{\mathrm{e}}){\mathbf{\nabla }{J}_{\nu }|}_{r=R_{\mathrm{out}}}.\end{array}$$(28)

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 Figs. 2 and 3.

Fig. 2. Non-grey case: normalised temperature profiles (upper panels) and SEDs (lower panels) for four different opacities τν0 = 1, 10, 100, and 1000 (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 Ivezic et al. (1997).

Fig. 3. Grey case: normalised temperature profiles (left panel) and emerging fluxes (right panel) for four different opacities τ = 0.01, 1, 10, and 100 (blue, orange, green and red, respectively) with a constant density profile (p = 0). The solid lines represent the FLD curves, and the black dots with the error bars σ indicate the MC profiles from Niccolini & Alcolea (2006).

3.4. 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,

$$\begin{array}{cc}\hfill F_{\nu }^{\mathrm{fall}}(R_{\mathrm{in}})& =2\pi J_{\nu }(R_{\mathrm{in}}){\displaystyle \underset{-1}{\overset{-{\mu }_0}{\int }}}\mu {\psi }_{\nu }(R_{\mathrm{in}},\mu )\mathrm{d}\mu \hfill \\ \hfill & =\frac{2\pi J_{\nu }(R_{\mathrm{in}})}{R_{\nu }tanh{R}_{\nu }}[ln\left(\frac{1+tanh{R}_{\nu }}{1+{\mu }_{0}tanh{R}_{\nu }}\right)-(1-{\mu }_0)tanh{R}_{\nu }],\hfill \end{array}$$(29)

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

$$\begin{array}{c}\hfill \sigma T_{\mathrm{eff}}^4=\sigma T_{\star }^4+{\left(\frac{R_{\mathrm{in}}}{R_{\star }}\right)}^2{\displaystyle \underset{0}{\overset{\infty }{\int }}}{F}_{\nu }^{\mathrm{fall}}(R_{\mathrm{in}})d\nu .\end{array}$$(30)

We imposed a fixed value for T_eff, and the temperature of the star T_⋆ was updated accordingly.

4. 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,

$$\begin{array}{c}\hfill \frac{1}{r^2}{\partial }_r\left(r^2{D}_{\nu }{\partial }_r{J}_{\nu }\right)={\kappa }_{\nu }^{\mathrm{abs}}(J_{\nu }-B_{\nu }).\end{array}$$(31)

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

$$\begin{array}{c}\hfill {\displaystyle \underset{0}{\overset{\infty }{\int }}}{\kappa }_{\nu }^{\mathrm{abs}}{B}_{\nu }\mathrm{d}\nu ={\displaystyle \underset{0}{\overset{\infty }{\int }}}{\kappa }_{\nu }^{\mathrm{abs}}{J}_{\nu }\mathrm{d}\nu .\end{array}$$(32)

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_⋆.

4.1. Numerical scheme

4.1.1. 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:

$$\begin{array}{c}\hfill A_{k,i+\frac{1}{2}}{J}_{k,i+1}-A_{k,i}{J}_{k,i}+A_{k,i-\frac{1}{2}}{J}_{k,i-1}=-b_{k,i}{B}_k(T_i).\end{array}$$(33)

The non-linear coefficients A are given by

$$\begin{array}{c}\hfill \begin{array}{c}\hfill A_{k,i\pm \frac{1}{2}}=r_{i\pm \frac{1}{2}}^2{D}_{k,i\pm \frac{1}{2}}{\frac{\mathrm{d}x}{\mathrm{d}r}|}_{i\pm \frac{1}{2}},A_{k,i}=A_{k,i+\frac{1}{2}}+b_{k,i}+A_{k,i-\frac{1}{2}},\\ \hfill b_{k,i}=\mathrm{\Delta }{x}^2{\left({\frac{\mathrm{d}x}{\mathrm{d}r}|}_i\right)}^{-1}{r}_i^2{\kappa }_{k,i}^{\mathrm{abs}}.\end{array}\end{array}$$(34)

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

$$\begin{array}{c}\hfill J_{k,i+\frac{1}{2}}=\frac{1}{2}(J_{k,i+1}+J_{k,i}),{\mathrm{\nabla }{J}_{\nu }|}_{k,i+\frac{1}{2}}={\frac{\mathrm{d}x}{\mathrm{d}r}|}_{i+\frac{1}{2}}\frac{J_{k,i+1}-J_{k,i}}{\mathrm{\Delta }x}\widehat{r}.\end{array}$$(35)

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

$$\begin{array}{c}\hfill J_{k,i}^{n+1}=\frac{b_{k,i}{B}_k(T_i^n)+A_{k,i+\frac{1}{2}}^n{J}_{k,i+1}^n+A_{k,i-\frac{1}{2}}^n{J}_{k,i-1}^{n+1}}{A_{k,i}^n},\end{array}$$(36)

and the temperature is updated after only one Gauss-Seidel spatial sweep at each frequency through Eq. (32), which we rewrite as

$$\begin{array}{c}\hfill {\displaystyle \sum _{k=0}^{n_{\nu }-1}}{W}_k{\kappa }_{k,i}^{\mathrm{abs}}{B}_k(T_i^{n+1})={\displaystyle \sum _{k=0}^{n_{\nu }-1}}{W}_k{\kappa }_{k,i}^{\mathrm{abs}}{J}_{k,i}^{n+1}.\end{array}$$(37)

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.

4.2. 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,

$$\begin{array}{c}\hfill T_{\star }^{n+1}={(T_{\mathrm{eff}}^4-\frac{1}{\sigma }{\left(\frac{R_{\mathrm{in}}}{R_{\star }}\right)}^2{\displaystyle \sum _{k=0}^{n_{\nu }-1}}{W}_k{F}_k^{\mathrm{fall},n+1})}^{\frac{1}{4}}.\end{array}$$(38)

Here again, we replaced the frequency integration by a quadrature formula with the associated weights W_k. $F_k^{\text{fall},n+1}$ is given by Eq. (29) and computed with the freshly updated values of $J_{k,\frac{1}{2}}^{n+1}$ and $R_{k,\frac{1}{2}}^{n+1}$.

4.3. Boundary conditions

We used two ghosts cells (one at each boundary of the domain) in order to simplify the implementation of our BCs. In doing so, the inner BC was imposed at the wall between the first (i = 0) and second cell (i = 1), and at the outer BC between the last (i = n_x − 1) and penultimate cell (i = n_x − 2). As indicated by Eqs. (20) and (27), we need to specify J_ν and ∇J_ν at these interfaces. For this, we used Eq. (35) and write

$$\begin{array}{cc}\hfill J_{\nu }(R_{\mathrm{in}}|R_{\mathrm{out}})& \approx \frac{J_{k,0|n_x-2}+J_{k,1|n_x-1}}{2},\hfill \\ \hfill \widehat{s}.\mathbf{\nabla }{J}_{\nu }{|}_{R_{\mathrm{in}}|R_{\mathrm{out}}}& \approx {\frac{\mathit{dx}}{\mathit{dr}}|}_{\frac{1}{2}|n_x-\frac{3}{2}}\frac{J_{k,1|n_x-1}-J_{k,0|n_x-2}}{\mathrm{\Delta }x}.\hfill \end{array}$$(39)

Accordingly, the values of J_ν in the ghosts cells were updated immediately after one Gauss-Seidel sweep Eq. (36) to ensure the BCs,

$$\begin{array}{cc}\hfill J_{k,0}& =\frac{2\mathrm{\Delta }x{\left(\frac{R_{\star }}{R_{\mathrm{in}}}\right)}^2{B}_k(T_{\star })-[\mathrm{\Delta }x({\alpha }_k-{\gamma }_k)-2\frac{dx}{dr}{|}_{\frac{1}{2}}{\beta }_k]J_{k,1}}{\mathrm{\Delta }x({\alpha }_k-{\gamma }_k)+2\frac{\mathrm{d}x}{\mathrm{d}r}{|}_{\frac{1}{2}}{\beta }_k},\hfill \\ \hfill J_{k,n_x-1}& =\frac{[2\frac{\mathrm{d}x}{\mathrm{d}r}{|}_{n_x-\frac{3}{2}}\zeta (R_k^{\mathrm{e}})D_k(R_k^{\mathrm{e}})-\mathrm{\Delta }x]J_{k,n_x-2}}{2\frac{\mathrm{d}x}{\mathrm{d}r}{|}_{n_x-\frac{3}{2}}\zeta (R_k^{\mathrm{e}})D_k(R_k^{\mathrm{e}})+\mathrm{\Delta }x},\hfill \end{array}$$(40)

with α_k, β_k and γ_k being defined by Eqs. (13) and (19). For the extrapolated value $R_k^{\text{e}}$ in $\zeta (R_k^{\text{e}})$ Eq. (26) and $D_k(R_k^{\text{e}})$, we used a second-order Lagrange extrapolation,

$$\begin{array}{c}\hfill R_k^{\mathrm{e}}=3(R_{k,n_x-\frac{5}{2}}-R_{k,n_x-\frac{7}{2}})+R_{k,n_x-\frac{9}{2}}.\end{array}$$(41)

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

$$\begin{array}{c}\hfill J_{k,i}^0=\frac{1}{4}{\left(\frac{R_{\star }}{r_i}\right)}^2{B}_k(T_{\star }),\end{array}$$(42)

from which we deduce the corresponding temperature profile T⁰ with the help of Eq. (37).

5. Numerical tests: spherically symmetric envelopes

5.1. Benchmarks from Ivezic et al. (1997)

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 Ivezic et al. (1997). 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_⋆ = 2500 K. This envelope extends from the inner radius R_in to the outer radius R_out = 1000 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 a power law of the form n(r) = n₀(R_in/r)^p. The radial optical depth τ_ν of the envelope is linked to the density profile by

$$\begin{array}{c}\hfill {\tau }_{\nu }={\displaystyle \underset{R_{\mathrm{in}}}{\overset{R_{\mathrm{out}}}{\int }}}{\kappa }_{\nu }^{\mathrm{ext}}\mathrm{d}r={\displaystyle \underset{R_{\mathrm{in}}}{\overset{R_{\mathrm{out}}}{\int }}}{C}_{\nu }^{\mathrm{ext}}n(r)\mathrm{d}r,\end{array}$$(43)

where $C_{\nu }^{\text{ext}}$ is the extinction cross-section coefficient. It is defined by

$$\begin{array}{c}\hfill \begin{array}{cc}{C}_{\nu }^{\mathrm{abs}}=C_{{\nu }_0}^{\mathrm{abs}},C_{\nu }^{\mathrm{sca}}=C_{{\nu }_0}^{\mathrm{sca}}\hfill & \mathrm{if}\nu \ge {\nu }_0,\hfill \\ C_{\nu }^{\mathrm{abs}}=C_{{\nu }_0}^{\mathrm{abs}}\left(\frac{\nu }{{\nu }_0}\right),C_{\nu }^{\mathrm{sca}}=C_{{\nu }_0}^{\mathrm{sca}}{\left(\frac{\nu }{{\nu }_0}\right)}^4\hfill & \mathrm{if}\nu \le {\nu }_0,\hfill \\ C_{\nu }^{\mathrm{ext}}=C_{\nu }^{\mathrm{sca}}+C_{\nu }^{\mathrm{abs}},\hfill \end{array}\end{array}$$(44)

with $C_{{\nu }_0}^{\text{abs}}=(1-\eta )C_{{\nu }_0}^{\text{ext}}$, $C_{{\nu }_0}^{\text{sca}}=\eta C_{{\nu }_0}^{\text{ext}}$, ν₀ 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 = 1, 10, 100,  and 1000. This created eight different cases to test the accuracy of our code. The coefficient n₀ in the density profile is derived with the help of τ_ν0 and p,

$$\begin{array}{c}\hfill n_0=\frac{(p-1){\tau }_{{\nu }_0}}{C_{{\nu }_0}^{\mathrm{ext}}{R}_{\mathrm{out}}}{\left(\frac{R_{\mathrm{out}}}{R_{\mathrm{in}}}\right)}^p{[{\left(\frac{R_{\mathrm{out}}}{R_{\mathrm{in}}}\right)}^{p-1}-1]}^{-1}.\end{array}$$(45)

The normalised temperature profile T/T_in and the normalised SED λF_λ/F ($F={\int }_0^{\infty }{F}_{\lambda } \text{d}\lambda$) of the envelope are shown in Fig. 2 for each case.

The Ivezic benchmarks were produced with version 2 of DUSTY³ (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⁻⁶ 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 = 1000). 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%.

5.2. Grey spherical shell with the Monte Carlo code from Niccolini & Alcolea (2006)

We compared the FLD code with a 3D MC radiative transfer code (Niccolini & Alcolea 2006). We wished to test the new BCs 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 temperature 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 differences 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_λ),

$$\begin{array}{c}\hfill \text{mean}(ϵ(T))=\frac{\sum _{i=0}^{N_x}{W}_iϵ(T_i)}{{\displaystyle \sum _{i=0}^{N_x}}{W}_i},\text{mean}(ϵ(F_{\lambda }))=\frac{\sum _{k=0}^{N_{\lambda }}{W}_kϵ(F_k)}{{\displaystyle \sum _{k=0}^{N_{\lambda }}}{W}_k}\end{array}$$(46)

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

$$\begin{array}{c}\hfill W_i={\left(\frac{\sigma (T_i)}{T_i}\right)}^{-2},W_k={\left(\frac{\sigma (F_k)}{F_k}\right)}^{-2}\end{array}$$(47)

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,

$$\begin{array}{c}\hfill J_{\nu }(R_{\mathrm{out}})+2D_{\nu }\widehat{s}.\mathbf{\nabla }{J}_{\nu }{|}_{r=R_{\mathrm{out}}}=0.\end{array}$$(48)

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

$$\begin{array}{c}\hfill \frac{{\int }_{{\mu }_0}^1{B}_{\nu }(T_{\star })\mathrm{d}\mu }{{\int }_{{\mu }_0}^1\mu B_{\nu }(T_{\star })\mathrm{d}\mu }=\frac{2(1-{\mu }_0)}{1-{\mu }_0^2}=\frac{2}{1+{\mu }_0}\end{array}$$(49)

with ${\mu }_0=\sqrt{1-{\left(\frac{R_{\star }}{R_{\text{out}}}\right)}^2}$, the cosine of the stellar angular size at R_out. As μ₀ → 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 (μ₀ → 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.

Fig. 4. Comparison of the outer boundary conditions presented here (solid lines) and from 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).

6. 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 ray-tracing module, which is a significant gain in computational time.

The next step, which will be our main concern for our future work, is the generalisation of these boundary conditions to non-spherically symmetric media, in particular for circumstellar discs, where the BCs will take a more complex form. In this regard, the extension to 2D will require numerical optimisations. Several directions of improvements are already being studied, such as the use of multi-frequency adaptive spatial grids, or acceleration procedures for efficiently solving the FLD equation under the radiative equilibrium condition. The last point is also crucial for solving extreme optically thick envelopes.

Acknowledgments

This study has been supported by the Lagrange laboratory of Astrophysics and funded under a 3-year PhD grant from Ecole Doctorale Sciences Fondamentales et Appliquées (EDSFA) of the Université Côte-d’Azur (UCA). The benchmark profiles were generated with the radiative transfer code DUSTY, available at http://faculty.washington.edu/ivezic/dusty_web. Part of the computations were carried out with the help of OPAL-Meso computing facilities. The authors are grateful to the OPAL infrastructure from Observatoire de la Côte d’Azur (CRIMSON) for providing resources and support.

References

All Tables

All Figures

	Fig. 1. Geometry of an inner spherical cavity (solid black line) illuminated by a central star (grey disc).
In the text

	Fig. 2. Non-grey case: normalised temperature profiles (upper panels) and SEDs (lower panels) for four different opacities τν0 = 1, 10, 100, and 1000 (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 Ivezic et al. (1997).
In the text

	Fig. 3. Grey case: normalised temperature profiles (left panel) and emerging fluxes (right panel) for four different opacities τ = 0.01, 1, 10, and 100 (blue, orange, green and red, respectively) with a constant density profile (p = 0). The solid lines represent the FLD curves, and the black dots with the error bars σ indicate the MC profiles from Niccolini & Alcolea (2006).
In the text

Fig. 4. Comparison of the outer boundary conditions presented here (solid lines) and from 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).

In the text