A new formal solution of the radiative transfer in arbitrary velocity fields
aa r X i v : . [ a s t r o - ph . S R ] F e b Astronomy&Astrophysicsmanuscript no. paper c (cid:13)
ESO 2018October 29, 2018
A new formal solution of the radiative transfer in arbitrary velocityfields.(Research Note)
Sebastian Knop , Peter H. Hauschildt , and Edward Baron , Hamburger Sternwarte, Gojenbergsweg 112, 21029 Hamburg, Germanye-mail: [sknop,yeti]@hs.uni-hamburg.de Homer L. Dodge Department of Physics and Astronomy, University of Oklahoma, 440 W Brooks, Rm 100, Norman, OK 73019-2061 USAe-mail: [email protected]
ABSTRACT
Aims.
We present a new formal solution of the Lagrangian equation of radiative transfer that is useful in solving the equation ofradiative transfer in the presence of arbitrary velocity fields.
Methods.
Normally a term due to the inclusion of the wavelength derivative in the Lagrangian equation of radiative transfer isassociated with a generalised opacity. In non-monotonic velocity fields, this generalised opacity may become negative. To ensure thatthe opacity remains positive, this term of the derivative is included in the formal solution of the radiative transfer problem.
Results.
The new definition of the generalised opacity allows for a new solution of the equation of radiative transfer in the presence ofvelocity fields. It is especially useful for arbitrary velocity fields, where it e ff ectively prevents the occurrences of negative generalisedopacities and still allows the explicit construction of the Λ -operator of the system needed for an accelerated Λ -iteration. We performedtest calculations, where the results of old, established solutions were compared with the new solution. The relative deviations neverexceeded 1 % and so the new solution is indeed suitable for use in radiative-transfer modelling.Non-monotonic velocity fields along photon paths frequently occur in three-dimensional hydrodynamical models of astrophysicalatmospheres. Therefore, the formal solution will be of use for multidimensional radiative transfer and has immediate applications inthe modelling of pulsating stars and astrophysical shock fronts. Key words. radiative transfer
1. Introduction
The Lagrangian equation of radiative transfer for moving atmo-spheres (Mihalas 1980) has been successfully solved for di ff er-ent astrophysical atmospheres, in particular supernova and novaatmospheres, as well as stellar atmospheres with winds.The inclusion of velocity fields in the equation of transfergives rise to two terms, one which is proportional to the spe-cific intensity. This term is of the same form as the term aris-ing from the opacity. Thus, the term from the velocity field canbe included in the definition of a generalised opacity (see forinstance Hauschildt 1992). For positive, velocity gradients, thecontribution of the velocity field to the generalised opacity isalways positive. In the case of negative, velocity gradients, thevelocity field contribution may become dominant and produce anegative, generalised opacity. Thus, the method of solution mustavoid exponentially increasing terms in the formal solution be-cause such terms make the numerical computation unstable.We show how to avoid this problem by not introducing ageneralised opacity but instead including the velocity field inthe formal solution of the radiative transfer. Negative opacitywill occur in non-monotonic velocity fields, thus any numeri-cal scheme for the formal solution must be able to handle bothpositive and negative velocity gradients. Numerically, the for-mal solution of the radiative transfer problem is then no longeran initial-value problem but a boundary-value problem. A formalsolution of this situation was described by Baron & Hauschildt (2004) and the revised formal solution in this paper is derivedfrom that solution. In addition, the formal solution incorporatesthe two discretisation schemes for the wavelength derivative andtheir mixing described by Hauschildt & Baron (2004). It it im-portant to note that the new formal solution must be consistentwith the solution to the scattering problem, e.g. an Accelerated Λ -iteration (ALI) (Scharmer 1981) as used in this work.An immediate application of this new solution is to bothone-dimensional model atmospheres of pulsating giant stars andshock fronts of all kinds, such as the interaction of a supernovaenvelope with the interstellar medium producing both forwardand reverse shocks. In the course of multidimensional modellingof astrophysical flows, the use of these solutions will become un-avoidable since arbitrary velocity fields occur frequently, due toe.g. Rayleigh-Taylor instabilities. In fact it is impossible to per-form comoving radiative transfer in arbitrary velocity fields bymeans of a characteristic method without the use of a formal so-lution that totally avoids the occurrence of negative generalisedopacities.We derive the revised formal solution with a di ff erent defini-tion of the opacity in Sect. 2. In Sect. 3, we demonstrate that thissolution produces results consistent with established solutionsand present brief conclusions in Sect. 4. Knop et al.: A new formal solution of the radiative transfer in arbitrary velocity fields. (RN)
2. The construction of the formal solution
We solve the equation of radiative transfer by following the so-lution along characteristic rays. In this form, the transport of thespecific intensity I is described along a photon path (with pathlength s ) and is independent of the dimensionality of the under-lying atmosphere:d I λ d s = η λ − χ λ I λ − a λ I λ − a λ ∂ ( λ I λ ) ∂λ , (1)where η is the emissivity and χ is the total extinction opacity.The dependence on the comoving wavelength has been explic-itly separated from the di ff erential dd s λ giving rise to the 4 a λ I λ term that can be incorporated into a generalised opacity. Theterm a λ = ∂λ∂ s is due to the wavelength coupling and is, for in-stance, given in spherical symmetry by: a λ = γ γ µ ( µ + β ) ∂β∂ r + β (1 − µ ) r ! , (2)where β = v / c , γ = / p − β , and r and µ = cos ϑ arethe phase-space coordinates (wavelength is of course a phase-space coordinate). We note that it is the sign of the wavelength-coupling term, and not the sign of the velocity gradient, that de-termines the upwind sense of the wavelength derivative.The form of Eq. 1 is the same for all geometries and dimen-sionalities. However, the form of the di ff erential dd s and the cou-pling term a λ will di ff er.To formally solve Eq. 1, the wavelength derivative ∂ ( λ I ) ∂λ mustbe discretised. We adopt the notation of Hauschildt & Baron(2004), where two possible discretisations of the wavelengthderivative are combined in a Crank-Nicholson-like scheme.Since we allow for arbitrary velocity fields, the formal so-lution must allow for an arbitrary sense of wavelength deriva-tive that is, the upwind direction may correspond either to longeror shorter wavelengths. To insure the stability of the discretisa-tion, local upwind schemes are introduced (Baron & Hauschildt2004) depending on the coupling term a λ . We assume a sortedwavelength grid λ l − < λ l < λ l + , and the wavelength depen-dence is represented by the wavelength index. The wavelengthderivative can then be written as: ∂ ( λ l I l ) ∂λ = λ l I l − λ l − I l − λ l − λ l − : a λ l ≥ λ l I l − λ l + I l + λ l − λ l + : a λ l < = h p − l I l − + p | l I l + p + l I l + i where the p • l are defined as: p − l = − λ l − λ l − λ l − p | l = λ l λ l − λ l − p + l = a λ l ≥ , p − l = p | l = λ l λ l − λ l + p + l = − λ l + λ l − λ l + a λ l < p • l depend not only on wavelength butalso on spatial position since the sign of the coupling term a λ may change along the characteristic.For mixing parameters of ξ ∈ [0 ,
1] for the two di ff erentdiscretisations, the equation of radiative transfer then reads:d I l d s = η l − χ l I l − a l (cid:16) + ξ p | l (cid:17) I l − ξ a l (cid:16) p − l I l − + p + l I l + (cid:17) − [1 − ξ ] a l (cid:16) p − l I l − + p | l I l + p + l I l + (cid:17) (4) See Chen et al. (2007) for a more complete description.
To solve Eq. 4 for monotonic velocity fields, it is customaryto define a generalised opacity:ˆ χ l = χ l + a l (cid:16) + ξ p | l (cid:17) (5)This approach in general, however, often fails because thegeneralised opacity in Eq. 5 may become negative if: − a l > χ l + a l ξ p | l (6)(Note that the p | l coe ffi cient always has the same sign as a l ). Forstrong wavelength couplings — as for instance in shock flows— the condition 6 is easily fulfilled. To avoid negative opticaldepths along the characteristics, the generalised opacity must bedefined di ff erently. In the ξ = p | l vanish. In principle,one could consider a method of correcting the formal solution atthese points, but the corrections must then also be applicable tothe explicit construction of the Λ -operator.A positive generalised opacity is assured by the followingdefinition:ˆ χ l = χ l + ξ a l p | l , (7)which only incorporates the physical opacity and a purely pos-itive contribution from the velocity field. The remaining term4 a l I l must then be incorporated into the formal solution of theequation of radiative transfer. This term then has the same e ff ectas if it was included in the opacity definition. For positive a l , thecontribution to the formal solution is negative and therefore de-creases the intensity along the ray in a way similar to an increasein opacity provides. For negative a l , the contribution is positivewhich behaves like a negative opacity. This new solution indeedprovides identical results to the previous discretisation method(see Sect. 3).After integrating Eq. 4 along a photon path from path length s to s , the formal solution between two spatial points on a char-acteristic can then be written in terms of optical depth, since theoptical depth τ λ along the characteristic relates to the path lengthby means of d τ l = ˆ χ l d s . I l ( τ ) = I l ( τ ) e τ − τ + Z τ τ ˆ S l e τ − τ d τ − Z τ τ ˜ S l e τ − τ d τ (8)where τ i = τ l ( s i ) andˆ S l = χ l ˆ χ l S l − ξ a l χ l (cid:16) p − l I l − + p + l I l + (cid:17)! (9)˜ S l = a l ˆ χ l (cid:16) [1 − ξ ] (cid:16) p − l I l − + p + l I l + (cid:17) + h + [1 − ξ ] p | l i I l (cid:17) (10)In a discrete atmosphere model, the integrals in Eq. 8 can beevaluated by means of piecewise linear- or parabolic interpola-tion of the auxiliary source functions: Z τ i τ i − ˆ S l e τ − τ i d τ ≈ α i , l ˆ S i − , l + β i , l ˆ S i , l + γ i , l ˆ S i + , l (11) Z τ i τ i − ˜ S l e τ − τ i d τ ≈ α i , l ˜ S i − , l + β i , l ˜ S i , l (12)where the coe ffi cients α , β , and γ are described inOlson & Kunasz (1987) and Hauschildt (1992). Note that thecontribution from the explicit wavelength derivative has beenonly linearly interpolated, whereas the implicit part can also beparabolically interpolated. nop et al.: A new formal solution of the radiative transfer in arbitrary velocity fields. (RN) Fig. 1: Comparison of the comoving spectra for the cases ξ = ξ = km / s around a fictitious line center.The line is omitted and the continuum is purely absorptive.The construction of the matrix equation for the formal solu-tion as well as the analytic construction scheme for a Λ opera-tor described in Baron & Hauschildt (2004) can be used withoutmodification in the new formal solution. These steps will there-fore not be repeated here.
3. The test calculations
To compare the new solution with other solutions, we used asimplified test setup (see also Hauschildt & Baron (2004) andKnop et al. (2007)). The opacity consists of a grey continuumand a single atomic line. The continuum opacity χ κ is grey andvaries with the radial structure according to χ κ ∝ r − . The ex-pression χ κ = κ κ + σ κ includes absorption κ κ as well as scattering σ κ , which are related by a thermalisation parameter ǫ κ = κ κ /χ κ .The corresponding opacities κ line , σ line , and thermalisation ǫ line are used for the — assumed to be Gaussian shaped — spectralline of the two-level-atom, where its strength is parameterisedrelative to the continuum opacity R = χ line /χ κ . The atmospherehas an extension of r min < r <
101 , where r min = cm. Aradial optical-depth scale is mapped onto the radial grid with thecontinuum opacity and ranges from 10 − < τ < . To includewavelength couplings in the atmosphere, a velocity field is im-posed on the radial structure.At first, we checked that, for an absent velocity field, the oldand new formal solutions provide identical results. Then we usedmonotonic velocity fields, in order to compare the new formalsolutions with the well tested solution described in Hauschildt(1992). For the velocity field, a linear law v ( r ) = v max r max / r with v max = / s was assumed.In the upper panel of Fig. 1a, the comoving frame spectrumfor a pure absorptive continuum without a line — σ κ = R = ξ =
0. The new solution is plotted witha dashed black line, while the old solution is represented by asolid grey line. The lower panel of Fig. 1a illustrates the dif-ference between the two solutions, implying that both solutionsreproduce the expected flat continuum and di ff er about 1%. Forcomparison, the computations in Fig. 1a were repeated for thecase ξ = ξ = ffi culties in re-producing a flat continuum because to its di ff usive behaviour, asalso demonstrated by Hauschildt & Baron (2004). This is also (a) (b) Fig. 2: Same as in Fig. 1, apart from the fact that the spectral lineis omitted and the continuum has a thermalisation parameter of ǫ κ = − . (a) (b) Fig. 3: Same as in Fig. 1, apart from the fact that the spectralline is strongly scattering ( ǫ line = − ) and the continuum has athermalisation parameter of ǫ κ = − .true for the new solution, as can be seen in the upper panel ofFig. 1b; the new formal solution closely follows the old solu-tion, as is also evident in the lower panel of Fig. 1b. The relativedi ff erence never exceeds 1 %, which is less than the deviation ofthe continuum from flatness.Since no scattering was included in the opacity, there is noneed for an ALI and a single Λ -step will solve directly the sys-tem. The comparisons between the comoving spectra then di-rectly yield the di ff erences in the formal solutions.In the next test, a scattering continuum opacity — of ǫ κ = − — was added, but the line was still omitted. The results forthe case of ξ = ξ = ff erences are shown in thebottom panels. The di ff erences between the scattering continuaare of the order of only 10 − %. The di ff usive behaviour of thefully implicit solutions appears again in the deviations from aflat continuum.In further tests, a strong scattering line — with R =
100 and ǫ line = − — was added to the continuum opacity which is alsoa ff ected by additional scattering — ǫ κ = − . The results forthe ξ = Knop et al.: A new formal solution of the radiative transfer in arbitrary velocity fields. (RN)
Fig. 4: Alternating velocity field in terms of optical depth.profile is a typical scattering line in an expanding atmosphere.The apparent excellent match of the two solutions is verified bythe plot of the relative di ff erences in the lower panel of Fig. 3a.The maximal di ff erence is of the order of 10 − %.The same calculation has been performed for the ξ = ff usive behaviour of this wavelength discretisation.However, again the match between the new and the old solutionis excellent.The excellent agreement between the old and new comovingspectra in the presence of scattering — with relative di ff erencesof the order of 10 − % — may at first seem unusual becausethe di ff erences in the formal solution itself were of the orderonly of one percent (see Figs. 1a and 1b). However, in thesecases, no ALI was needed due to the lack of scattering. In thescattering cases, an ALI was performed. Since only exact ele-ments from the Λ -operator have been used in the constructionof the Λ ∗ -operator, the operators derived from the old and newformal solutions also di ff er. From the excellent agreement of thefinal results, it is obvious that these di ff erences cancel each otherand enable far closer agreement for the moments of the radiationfield in the ALI than the formal solution for a given source func-tion.In the preceding tests, there was no need to employ the newsolution because the velocity fields were monotonic. To testthe new solution for alternating wavelength couplings, a non-monotonic, alternating velocity field was used.To demonstrate the robustness of the new solution, we per-formed calculations for an alternating velocity profile — shownin Fig. 4 — that penetrates into the deep layers of the atmo-sphere. This run of the velocity profile was not even remotelypossible within the old framework, and hence no comparisonwith the old solution is shown. The comoving spectrum is shownin Fig. 5. The peculiar features and broadness of the line resultfrom the extreme velocity fields in the atmosphere. The ξ = ξ = a l I l term from the opacityand its addition as an explicit source term has apparently littleinfluence on the solution. In all of our test calculations, all dif-ferences with respect to established solutions were minor. In thecase of scattering in the line as well as in the continuum the dif-ferences become almost negligible. In comparison, the di ff erent Fig. 5: Co-moving spectrum for the presence of the non-monotonic velocity field shown in Fig. 4 for ξ =
1. The wave-length scale is given in km / s around a fictitious line center. Theline is strongly scattering ( ǫ line = − ) and the continuum has athermalisation parameter of ǫ κ = − .discretisations had a far more dominant e ff ect on the solutionthan the addition of the explicit term, as can be seen in the caseof the ξ = ff usive behaviour is tracedout. We should note further that the speed of convergence of theformal solution did not change significantly.The new solution method appears to be well suited to generalradiative-transfer modelling involving non-monotonic velocityfields. This has immediate applications to 1D spherical symmet-ric atmosphere models, e.g. for pulsating atmospheres of vari-able stars as well as shock fronts in supernova flows.It remains to be seen whether the di ff erent numerical natureof the solution, due to an additional explicit term, will be trou-blesome in some cases. Our test computations did not identifyan example of numerically unstable behaviour during the explo-ration of the parameter space.
4. Conclusions
We have presented a new formal solution of the equation of ra-diative transfer. Its advantage with respect to established solu-tions lies in the avoidance of negative opacities that occur by de-sign in non-monotonic velocity fields. It allows the computationof the arbitrarily wavelength-coupled, radiative-transfer equa-tion that arises with complete avoidance of negative, generalisedopacities that render the solution unstable. The new approach isrequired for any modelling that involves alternating wavelengthcouplings. Since one considers radiative transfer in multidimen-sional situations, non-monotonic velocity fields do naturally oc-cur.The relative di ff erences between the new and old solutionswas not found to exceed 1 % in all test calculations and werein most cases several orders of magnitude smaller. Therefore,the new solution is a viable replacement of the old solutions anduseful in modelling the atmospheres of variable stars or radiativetransfer in hydrodynamically obtained structures. Acknowledgements.
This work was supported in part by SFB 676 from theDFG, NASA grant NAG5-12127, NSF grant AST-0707704, and US DOE GrantDE-FG02-07ER41517. This research used resources of the National EnergyResearch Scientific Computing Center (NERSC), which is supported by theO ffi ce of Science of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098; and the H¨ochstleistungs Rechenzentrum Nord (HLRN). Wethank all these institutions for a generous allocation of computer time. References
Baron, E. & Hauschildt, P. H. 2004, A&A, 427, 987 nop et al.: A new formal solution of the radiative transfer in arbitrary velocity fields. (RN)5