Issue 
A&A
Volume 653, September 2021



Article Number  A139  
Number of page(s)  10  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202040236  
Published online  24 September 2021 
New boundary conditions for the approximate fluxlimited diffusion radiative transfer in circumstellar environments
Test case study for spherically symmetric envelopes
Université Côte d’Azur (UCA), Observatoire de la Côte d’Azur (OCA), CNRS, Laboratoire Lagrange, Lagrange, 06304 Nice Cedex 4, France
email: jeremy.perdigon@oca.eu
Received:
24
December
2020
Accepted:
2
July
2021
Context. 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.
Aims. In this context, the use of approximate methods is of primary interest. One promising candidate method is the fluxlimited diffusion (FLD), which recasts the radiative transfer problem into a nonlinear 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.
Methods. 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 nonlinear 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.
Results. 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. This feature avoids additional formal solutions for the radiative transfer equation in a set of rays (raytracing computations).
Conclusions. The FLD approximation together with the proposed new BCs performs well and captures the main physical properties of the radiative equilibrium in spherical circumstellar envelopes.
Key words: radiative transfer / methods: numerical / circumstellar matter
© J. Perdigon et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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 midInfrared SpectroScopic Experiment^{1} (MATISSE), operating in the midIR, or the Atacama Large Millimeter Array^{2} (ALMA) in the submillimetric domain offer complementary views of these environments, that give access to regions close to the central star up to the outer regions using a multiwavelength 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 fluxlimited diffusion (FLD), introduced by Levermore & Pomraning (1981) (L&P hereafter). This description numerically simplifies the problem by recasting the radiative transfer equation into a nonlinear 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 halfspace 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; MignonRisse 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 nongrey 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 selfheating 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 byproduct 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 , and the frequency by the subscript ν. The transport equation for the specific intensity at the position r in the direction with isotropic and coherent scattering is
J_{ν} = J_{ν}(r,t) is the mean specific intensity, B_{ν} = B_{ν}(T(r,t)) is the Planck function and , and 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 ( = 0). L&P showed that, using Eq. (6) in Eq. (5), h_{ν} is proportional to R_{ν},
where R_{ν} is the norm of R_{ν} and λ(R_{ν}) is the ‘fluxlimiting’ parameter. Finally, using Eq. (8) in Eq. (4) allows Eq. (3) to be rewritten as a nonlinear diffusion equation for the mean specific intensity
with the nonlinear 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 nonlinear 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 timeindependent 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 Robintype 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 entering a boundary surface at the position r_{S},
Here, is the outward normal vector to the surface of the envelope ( at the inner edge). The righthand 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 (), which is exact for spherically symmetric problems. The equation can then be rewritten as a nonlinear mixed BC,
with
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 sourcefree (B_{ν} = 0) halfspace media, with a constant , ω_{ν} 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 sourcefree (B_{ν} = 0) halfspace media, with a constant , ω_{ν} 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_{⋆} () and surface temperature T_{⋆}. For this, we need to specify the flux . We have two contributions,
The first contribution, denoted by B_{ν}(T_{⋆}) comes from the star and the other from the inner boundary itself and is expressed as . As shown by Fig. 1, the vector corresponds to the opposite point at the inner boundary, along . 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 . To perform the angular integration, we aligned the n_{z} axis with the unitary vector . For the star, the integration on μ = cos(θ) (θ being the angle between n_{z} and ) spans from to 1, and for the inner cavity, it spans from 0 to μ_{0},
Fig. 1.
Geometry of an inner spherical cavity (solid black line) illuminated by a central star (grey disc). 
The incident flux 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%).
3.2. Outer boundary
For the inner cavity, we would like to impose that no incoming radiation enters the external shell of the envelope (). 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 sharppeaked 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 . In the diffusion regime where the emergent field is isotropic, . 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 and . 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 secondorder extrapolated value for (see Eq. (41)). Because we prescribe a value for R_{ν} at the external edge, is also used to compute the nonlinear diffusion coefficient 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 . 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 nonlinear diffusion coefficient 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 Figs. 2 and 3.
Fig. 2.
Nongrey 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,
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 opticallythick 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.
4. Numerical implementation
The FLD Eq. (9) is a nonlinear 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 timeindependent 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_{⋆}.
4.1. Numerical scheme
4.1.1. Finitedifference approximation
Following a finitedifference 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 secondorder 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 nonuniformly sampled. We make use of one ghost cell for each grid border to ensure the BCs. We obtain the following system of equations:
The nonlinear coefficients A are given by
The nonlinear 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
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 GaussSeidel 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 GaussSeidel 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 lefthand side (LHS) of Eq. (37) was precomputed 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 NewtonRaphson 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 GaussSeidel 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,
Here again, we replaced the frequency integration by a quadrature formula with the associated weights W_{k}. is given by Eq. (29) and computed with the freshly updated values of and .
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
Accordingly, the values of J_{ν} in the ghosts cells were updated immediately after one GaussSeidel sweep Eq. (36) to ensure the BCs,
with α_{k}, β_{k} and γ_{k} being defined by Eqs. (13) and (19). For the extrapolated value in Eq. (26) and , we used a secondorder Lagrange extrapolation,
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 tradeoff 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).
5. Numerical tests: spherically symmetric envelopes
5.1. Benchmarks from Ivezic et al. (1997)
We tested the accuracy of our FLD code with our Robintype 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_{0}(R_{in}/r)^{p}. The radial optical depth τ_{ν} of the envelope is linked to the density profile by
where is the extinction crosssection coefficient. It is defined by
with , , ν_{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 1000. 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 () 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 nonphysical results of the DUSTY code below this threshold, for the smallest wavelengths.
Results from the comparison with DUSTY.
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_{λ}),
Results from the comparison with the MC code.
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 planeparallel geometry for an optically thin slab illuminated by an isotropic incident radiation field. However, for the special case of a spherical envelope surrounding a blackbody star, this ratio becomes
with , 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 planeparallel 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 sharppeaked 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 otherhand, 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 nongrey 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 semilog 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 nonlinear diffusion equation Eq. (31) coupled with the radiative equilibrium Eq. (32) was performed with a GaussSeidel methodbased iterative scheme, in which the temperature is updated at each iteration step. Furthermore, we point out that the FLD approximation implemented with the proposed outerboundary condition provides a simple approximation for the flux emitted by the envelope. This allows computing the SED or emerging flux without using any raytracing 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 nonspherically 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 multifrequency 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.
Available at http://faculty.washington.edu/ivezic/dusty_web/
Acknowledgments
This study has been supported by the Lagrange laboratory of Astrophysics and funded under a 3year PhD grant from Ecole Doctorale Sciences Fondamentales et Appliquées (EDSFA) of the Université Côted’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 OPALMeso 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
 Ivezic, Z., & Elitzur, M. 1997, MNRAS, 287, 799 [Google Scholar]
 Ivezic, Z., Groenewegen, M. A. T., Men’shchikov, A., & Szczerba, R. 1997, MNRAS, 291, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Kuiper, R., Klahr, H., Dullemond, C., Kley, W., & Henning, T. 2010, A&A, 511, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321 [NASA ADS] [CrossRef] [Google Scholar]
 MignonRisse, R., González, M., Commerçon, B., & Rosdahl, J. 2020, A&A, 635, A42 [CrossRef] [EDP Sciences] [Google Scholar]
 Niccolini, G., & Alcolea, J. 2006, A&A, 456, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Niccolini, G., Woitke, P., & Lopez, B. 2003, A&A, 399, 703 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pomraning, G. C. 1986, J. Quant. Spectr. Rad. Transf., 36, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Pomraning, G. C. 1988, J. Comput. Phys., 75, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Sonnhalter, C., Preibisch, T., & Yorke, H. W. 1995, A&A, 299, 545 [NASA ADS] [Google Scholar]
 Steinacker, J., Baes, M., & Gordon, K. D. 2013, ARA&A, 51, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Yorke, H. W., & Sonnhalter, C. 2002, ApJ, 569, 846 [NASA ADS] [CrossRef] [Google Scholar]
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.
Nongrey 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 semilog 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.