Graziani F. (editor) Computational Methods in Transport

Подождите немного. Документ загружается.

Radiative Transfer in Astrophysical Applications 17

The point is that χ and η generally depend on radiation intensity as well,

moreover in a complex non-linear and non-local fashion. We shall return to

this point in the next section.

We also note that it is common in astrophysical applications to introduce

the optical depth along a given direction as dτ = χds,andthesource function

S ≡

. (4)

3 Absorption, Emission and Scattering Coeﬃcients

The absorption and emission coeﬃcients are introduced phenomenologically,

and are deﬁned analogously to the speciﬁc intensity, namely as the energy

removed or added to a beam of radiation at unit frequency range, solid angle,

area, and time. It is customary to split them into two parts,

χ = κ

+ κ

, (5)

and

η = η

+ η

, (6)

where κ

represents the true absorption (a process where a photon is ab-

sorbed and subsequently destroyed by a thermal process such as a collisional

de-excitation); κ

is the scattering coeﬃcient (a photon is absorbed and sub-

sequently re-emitted). Similarly, η

is the coeﬃcient of thermal emission (a

photon is created at the expense of the thermal pool); and η

is the scat-

tering part of the emission coeﬃcient, corresponding to emission of a photon

following previous absorption.

Scattering absorption coeﬃcient may be written as

(ν)=



(ν) , (7)

where the summation extends over all species α that do scatter radiation, n

is the number density, and σ

the corresponding cross-section. For instance,

for hot star atmospheres (spectral types O and B), the dominant, and usu-

ally the only scattering process that is considered is the electron (Thomson)

scattering, in which case n

(ν)=n

; n

being the electron density, and

the Thomson cross-section. For even hotter objects (e.g., hot accretion

disks, T ≈ 10

K and higher), the dominant scattering process is the Comp-

ton scattering. For cooler stars (A, F, G), the Rayleigh scattering on atomic

hydrogen (and helium) becomes important; for even cooler stars Rayleigh

scattering on the hydrogen molecule H

and other molecules is important,

and for coolest stars and substellar mass objects (brown dwarfs and giant

planets), scattering on cloud (dust) particles becomes dominant.

The scattering part of the emission coeﬃcient may be written as

18 I. Hubeny

(ν)=





dν





dΩ



I(ν



, n



(ν



, n



,ν,n) , (8)

where R is the redistribution function, which has the meaning of the prob-

ability density that a photon is absorbed at frequency ν



and direction n



(within appropriate unit intervals around these values), and re-emitted at ν

and n.

The scattering part of the emission coeﬃcient is the culprit behind most

of the numerical complexity of the transfer problem. It makes the problem

pathologically implicit because in order to solve the transfer equation we

need to know the emission eﬃcient, but it (or at least its scattering part) is

a function of the speciﬁc intensity; moreover it couples many frequencies and

directions.

However, although η

is the most troubling part of the problem, we shall

show that all the other coeﬃcients do also generally depend on radiation ﬁeld,

and usually in a complex non-local and non-linear way. From the point of

view of microscopic physics, one can write down the absorption and emission

coeﬃcients as follows (we only indicate the dependence on frequency, and

omit an indication of their dependence on position and time)

(ν)=



j>i

− (g

] σ

(ν)+





− n

∗

−hν/kT



iκ

(ν)



κκ

(ν, T )



1 − e

−hν/kT



, (9)

where the three terms represent, respectively, the contributions of bound-

bound transitions (i.e. spectral lines), bound-free transitions (continua), and

free-free absorption (also called bremsstrahlung).

Here, n

is the occupation number (population) of an atom in the energy

level labeled i, g

the corresponding statistical weight, and n

∗

denotes an equi-

librium population of level i corresponding to temperature T, and density ρ.

σ(ν) are the corresponding cross-sections; subscript κ denotes the “contin-

uum,” and n

the ion number density. The relation between the bound-bound

cross section σ

(ν) and the well-known Einstein coeﬃcients for the for the

photo-excitation is σ

(ν)=(hν

/4π)B

φ(ν); φ(ν) is the so-called absorp-

tion proﬁle coeﬃcient, normalized to unity,



φ(ν) dν =1.Itrepresentsthe

conditional probability density that if a photon is absorbed in the transition

i → j, it is absorbed in the frequency range (ν, ν + dν).

Analogously, the thermal emission coeﬃcient is given by

(ν)=

2hν

⎡

⎣



j>i

)σ

(ν)+



∗

iκ

(ν)e

−hν/kT



κκ

(ν, T )e

−hν/kT



. (10)

Radiative Transfer in Astrophysical Applications 19

The three terms again describe the bound-bound, bound-free, and free-free

emission processes, respectively.

The absorption and emission coeﬃcients are thus described through the

corresponding cross sections—given by the atomic physics, and possibly also

by the local thermodynamic parameters, temperature and density; and the

atomic level populations for all the levels involved in the microscopic processes

that give rise to an absorption and emission at frequency ν.Theirnumber

may be enormous. Here is the essential complication. One can either assume

the Local Thermodynamic Equilibrium (LTE), in which case the level pop-

ulations are given through the Saha-Boltzmann formula, i.e. n

= n

∗

(T,ρ),

so that the thermal absorption and emission coeﬃcient are functions of only

temperature and density. However, in many cases of interest, the LTE approx-

imation breaks down, and one has to determine the atomic level populations

by a set of statistical equilibrium (rate) equations. For a given level i,therate

equation reads (for simplicity, we consider here the time-independent form;

the general form would contain additional terms (1/c)(∂n

/∂t)+∇·(vn

)).



j=i

+ C



j=i

+ C

) , (11)

where R

is the rate of radiative transition from level i to level j, C

the analogous collisional rate. Equation (11) thus expresses the principle of

statistical equilibrium, namely that the total number of transitions out of

level i is equal to the total number of transitions into level i (we note that

if this is satisﬁed for all microscopic processes i ↔ j we recover LTE). The

essential complication arises due to the fact that while the collisional rates

are functions of temperature and perturber density only, the radiative rates

are given by



dν



dΩ

4π

hν

(ν) I(ν, n) , (12)

(and analogously for the downward rate); that is, the level populations that

are needed for evaluating even the thermal absorption and emission coeﬃ-

cients are determined through the (yet unknown) radiation ﬁeld, moreover in

essentially all frequencies and directions. An approach where at least some

important level populations are determined by a self-consistent solution of the

statistical equilibrium equation and the radiative transfer equation is called

in astrophysics the non-LTE (or NLTE) approach, and its elaboration was in

fact the main topic of the ﬁeld of stellar atmospheres theory in the last three

decades.

Finally, I stress that even in the case of LTE, one still has to deal with

with complex non-local and non-linear coupling of radiation ﬁeld and the

state parameters. This is because the local temperature is determined through

some form of energy balance (radiative equilibrium for hot star atmospheres;

radiative + convective equilibrium in cooler stars; plus additional dissipation

of mechanical energy in case of stellar chromospheres or accretion disks, etc.).

20 I. Hubeny

In the simplest case of pure radiative equilibrium, which reads



∞



[χ(ν, T )I(ν, n) − η(ν, n,T)] dν dΩ =0, (13)

we see that the temperature structure depends on the radiation ﬁeld in a

complex way.

4 Hierarchies of Approximations

There is a wide range of diﬀerent radiative transfer problems considered in

astrophysics, depending on a degree of sophistication in treating the inter-

action of radiation and matter. I will brieﬂy summarize the basic types of

problems below, although the list is by no means exhaustive.

4.1 Treatment of Absorption and Emission

• Formal solution

This is the simplest case, in which the emission and absorption coeﬃcients

(or the optical depth and the source function) are fully known. In this case

the transfer equation is a linear diﬀerential equation, and its solution is

easy. The only concern is to ﬁnd a suﬃciently fast and eﬃcient numerical

procedure. The topic was recently reviewed for instance in [5] and [6].

The basic use of the formal solution is as an intermediate step of a global

iterative method, as we shall discuss later on.

• Monochromatic(coherent scattering) problems

In this case, χ and S may depend on frequency ν, but there is no coupling

to other frequencies ν



. In other words, all the frequencies are separate,

one one can solve for one frequency at a time. There is a subdivision

according to various possibilities of treating the angular coupling:

– anisotropic scattering

Here, the source function contains a term proportional to

(n) ∝



σ(n) p(n



, n) I(n



) dΩ



where p is the phase function.

– isotropic scattering, where p = 1, and thus the scattering part of the

source function is independent of direction and is given by

(n)=S

∝



I(n



) dΩ



≡ J,

where J is the zero-order moment of the speciﬁc intensity, called in

astrophysics the mean intensity of radiation (it is called the scalar ﬂux,

φ, in the neutron transport theory).

Radiative Transfer in Astrophysical Applications 21

• Non-coherent scattering problems

In this case, the scattering source function contains the term

(ν, n) ∝



R(ν



, n



,ν,n)I(ν



, n



) dν



dΩ



, (14)

where R is the redistribution function. Typical examples of the situations

where a complex frequency and angular redistribution of radiation is im-

portant are for instance Compton scattering in high energy astrophysics,

resonance scattering in spectral lines, or Mie scattering on cloud particles

for coldest stars and giant planets. Simpliﬁed cases include:

– angular and frequency redistribution are separate,

R(ν



, n



,ν,n)=

R(ν



,ν) p(n



, n) ,

which is a very good approximation in most cases, and which simpliﬁes

a numerical solution considerably.

– complete redistribution, in which cases the scattering is viewed as two

uncorrelated process of absorption and subsequent re-emission, with

the same probability densities,

R(ν



,ν)=φ(ν



) φ(ν) ,

where φ(ν) is the absorption proﬁle coeﬃcient. This is usually a good

approximation for scattering in spectral lines except strong resonance

lines when one has to deal with more general redistribution function.

In astrophysical terminology, the latter situation is called partial re-

distribution, in contrast to the complete redistribution.

– coherent scattering, which we mentioned earlier. corresponds to

R(ν



,ν)=φ(ν



) δ(ν



− ν) ,

where δ is the Dirac δ-function.

4.2 Treatment of Spatial Dimensions

• 1-D plane-parallel geometry

In this case one assumes that the medium is composed of plane-parallel,

horizontally-homogeneous layers. The problem is thus 1-dimensional, in

Cartesian geometry. This approach is usually used for stellar atmospheres.

Although it may seem very crude, it actually provides an excellent ap-

proximation for atmospheres whose thickness is much smaller that the

stellar radius, which is the case for most stellar types except giants and

supergiants with extended atmospheres. Also, such an approach is war-

ranted because some other essential pieces of physics, such as departures

from LTE, and/or a radiative/convective balance, are important, which

22 I. Hubeny

leads to a necessity of treating a large number of frequencies and possibly

directions. This approach can be used even in the case of surface inho-

mogeneities, as long as the radial properties vary much faster than the

horizontal properties (which, again, is usually the case). 1-D plane-parallel

approach is also often being used for geometrically thin accretion disks.

There are several codes for this type of problems; the most widely used are

a universal stellar atmosphere and accretion disk program TLUSTY [7–9],

and the package PRO2 [10, 11].

• 1-D spherical geometry

It is used for extended atmospheres, stellar winds, supernovae, H II re-

gions, and others. Although it is still a one-dimensional problem, it is

numerically much more demanding than the previous 1-D plane-parallel

problem because one has to take into account many more angles (di-

rections), because the radiation intensity is becoming more and more

forward-peaked when going outward. The most sophisticated codes of

this sort are CMFGEN [12], and PHOENIX [13].

• 2-D spherical geometry

It is used for rapidly rotating stars, core-collapse supernovae, etc.

• 2-D cylindrical geometry

It is mostly used for accretion disks and core-collapse supernova simula-

tions. A well-described code of this sort is ALTAIR [14, 15].

• 3-D Cartesian geometry

So far it has been used for snapshots of 3-D hydrodynamic simulations “in

the box” for solar convection, modeling viscous dissipation in accretion

disks, mass transfer in close binary systems, simulating large-scale cosmo-

logical structures, and others. Obviously, the problem is computationally

extremely demanding, so the 3-D transfer simulation use a number of

various approximations (e.g., [16]).

4.3 Treatment of frequency

Ideally, one has to treat as many frequency points as to faithfully reproduce

the frequency dependence of the absorption (and emission) coeﬃcient. While

the bound-free and free-free processes have relatively smooth cross-sections

(although, strictly speaking, the photoionization cross-section do sometimes

contain a number of resonances, so they do not be smooth after all) and

thus can be reproduced by a relatively small number of discretized frequency

points, the bound-bound processes (spectral lines) have very sharply peaked

cross-sections. Moreover, a typical stellar spectrum contains literally milions

of spectral lines, which should all be consedred explicitly in order to construct

an accurate model atmosphere (essentially because the rate equations and the

energy balance equation contain integrals of the radiation intensity over the

whole frequency range). Model stellar atmospheres (and other astronomical

radiation-dominated objects) either ignore the eﬀects of the majority of spec-

tral lines and treat only the most important ones; or do treat, with a varying

Radiative Transfer in Astrophysical Applications 23

degree of approximation, essentially “all” spectral lines. Model atmospheres

of the latter kind are called line-blanketed models. I list below several op-

tions for treating the frequency dependence of relavant quantities, ordered

by increasing accuracy, and by the number of frequency points required.

• Gray problem

In this case, there is only one frequency; better speaking, one either

assumes that the absorption and emission coeﬃcients are frequency-

independent, or one works in terms of frequency-integrated quantities.

In the context of stellar atmospheres, such an approach is being used as

an interesting pedagogical tool for understanding the temperature struc-

ture, and for providing a starting solution for an iteration scheme to solve

a general stellar atmosphere problem.

• Multi-group/multi-gray

is a generalization of the gray model. One introduces several, O(1) –

O(10), frequency bins (groups), and averages all the relevant frequency-

dependent quantities over the individual bins/groups. Again, this ap-

proach is no longer being used in the context of stellar atmosphere models;

although [17] presented a generalization this approach for NLTE model

atmospheres with a treatment of many spectral lines; however, the num-

ber of bins was signiﬁcantly higher – O(10

)toO(10

). The multi-group

approach is however still being used in radiation hydrodynamic studies (or

any problem which is already so computationally heavy that one cannot

aﬀord to deal with many frequencies), and, in particular, for treating neu-

trino transport in the core-collapse supernova simulations. In the latter

case, the relevant cross-sections are relatively smooth functions of fre-

quency/energy, so that such an approach is a very good approximation

anyway.

• Opacity Distribution Functions

One introduces a number of frequency intervals [typically O(100) - O(10

)];

instead of simply averaging over the intervals, one resamples the opacity

in the individual intervals into monotonic functions of frequency, which

can then be treated by means of a small number of frequency points [typ-

ically O(10) per interval]. One can thus represent the whole frequency

space in roughly O(10

) frequency points.

• Opacity Sampling

This is essentially a Monte-Carlo treatment of frequency dependence. One

selects a set of random frequency points that statistically reproduce a

complex frequency behavior of the absorption and emission coeﬃcient.

Obviously, the larger the number of frequency points, the higher the ac-

curacy of the solution. In typical (stellar atmosphere) applications, one

thus deal has to deal with O(10

)tofewtimesO(10

) frequency points.

• Full frequency resolution

Such an item would be missing in a similar review only a few years ago

because it would have been viewed as too computationally demanding.

24 I. Hubeny

However, such approaches are possible now, at least in 1-D. Recent grids

of NLTE line-blanketed model stellar atmospheres, such as the grid of

O star models [18] that uses Opacity Sampling with a resolution of 0.75

ﬁducial Doppler widths, can be viewed as more or less “full” frequency

resolution. This approach needs 2–3 ×10

frequency points; cooler star

models would need up to O(10

) frequency points, but this is becoming

possible with current computers.

4.4 Treatment of angular dimensions

Brieﬂy, there are two types of approaches. One can either treat an angular

dependence of speciﬁc intensity explicitly, or to use either variable Eddington

factors (or Eddington tensor in more than 1 spatial dimension), or some other

closure relations, and to solve the transfer equation only for the moments of

speciﬁc intensity, such as the mean intensity J. The latter approach is very

advantageous in the case of isotropic scattering. For the problems of general

anisotropic scattering, one usually adopts the former approach and solves the

transfer problem for the speciﬁc intensities.

In 1-D geometries one typically considers one angular variable (the polar

angle measured with respect to the normal to the surface), and assumes az-

imuthal symmetry of the radiation ﬁeld. In that case, one works in terms of

azimuthally-averaged speciﬁc intensity. In 1-D plane-parallel geometry, it is

usually suﬃcient to consider a very small number of polar angles (3 is a typ-

ical value in stellar atmosphere models). However, model atmospheres with

strong irradiation from outside (for instance for very close-in extrasolar giant

planets), and moreover with strongly anisotropic scattering phase function,

require at least several hundred angles (e.g., [19]).

1-D spherical models do typically require more angles that 1-D plane-

parallel models, usually there are roughly as many polar angles as the number

of radial zones. In multi-D transfer problems the number of angles depend on

actual situation; some workers in radiation hydrodynamic models for solar

atmosphere use only O(10) angles (see, e.g. [16]); but most transfer studies

in multi-D consider O(100) to O(10

) angles.

5 General Problem

5.1 State Vector

The state of a radiation-dominated medium, as discussed above, is fully de-

scribed by the following general state vector

ψ = {I(t, r,ν,n),T(t, r),ρ(t, r), v(t, r), [n

(t, r),i =1,NL]}, (15)

where the state parameters are, respectively, the speciﬁc intensity of radi-

ation, temperature, density, velocity vector, and the atomic (possibly also

Radiative Transfer in Astrophysical Applications 25

molecular) level populations (occupation numbers); here NLis the total num-

ber of energy levels for which the statistical equilibrium equation is solved,

and for which one thus allows departures from LTE. In some simpliﬁed situ-

ations, certain components may be missing. For instance, for static models,

v ≡ 0, so that the velocity is skipped from the state vector. Similarly, in LTE,

the level populations are known functions of temperature and density, and

thus are skipped from the state vector. In we consider the so-called restricted

NLTE problem, that is assuming that the temperature, density (and velocity

in the case of dynamic models) are known, one is left only with radiation

intensity and level populations to be solved for. Finally, if we assume LTE

for level populations, and assuming ﬁxed temperature, density, and velocity,

the only unknown quantity is the radiation intensity. This is actually the case

in many problems in the neutron transport theory and the Earth atmosphere

studies.

Let us ﬁrst summarize typical numbers of the components of the state

vector, in order to appreciate the numerical complexity of the problem.

• t – many time steps, depending on the problem. However, the vast ma-

jority of astrophysical radiation transport problem is time-independent;

• r – typically O(100

sp.dim.

) grid points where N

sp.dim.

is the number of

spatial dimensions;

• n – O(10

)toO(10

) grid points, depending on problem (see above);

• ν – O(10

)toO(10

) points, depending on problem (see above);

• NL – O(10

)toO(10

) levels, depending on problem.

Let us now give some examples of actual numbers used in the current

state of the art numerical simulations in astrophysics. We will only count the

number of individual components of discretized speciﬁc intensity; the number

of other components of the state vector is negligible. The total number is a

product of number of spatial positions, number of angles, and number of

frequencies. In the subsequent examples, we will follow this order.

• 1-D static plane-parallel NLTE line-blanketed model stellar atmospheres

× few × (few × 10

) ≈ 10

elements;

• 1-D close-in extrasolar giant planets (with anisotropic scattering)

× 10

) ≈ 10

elements;

• 2-D core-collapse supernovae (neutrino transport)

(few ×10

) ×10

×(few ×10

) ≈ 10

to 10

elements, but per time step;

there are O(10

) time steps.

These are already existing simulations. In the (hopefully near) future, one

hopes to achieve for instance the following:

• 3-D core-collapse supernovae (neutrino transport)

(few ×10

) × 10

× (few × 10

) ≈ 10

to 10

but per time step; ≈ 10

altogether;

26 I. Hubeny

• 3-D static plane-parallel NLTE line-blanketed model stellar atmospheres

× 10

≈ 10

elements.

5.2 Methods of Solution

There are several basic possibilities, which I outline below. Again, the list is

far from being complete, and is meant for illustration only. The approaches

are categorized by means of a degree of sophistication to which the interaction

between radiation and matter is treated numerically.

• not solving, really

This is obviously the simplest, and therefore most popular, approach to

the problem. There are several possibilities, including:

– optically thin approximation—χ is essentially zero; consequently no

interaction of radiation with matter.

– diﬀusion approximation—an almost opposite approximation, valid for

large optical depths. In this case, one can express the moments of

speciﬁc intensity through the gradient of temperature; radiation in-

tensity is thus eliminated from the state vector, but the results are

approximate, the more so in the optically thin regime.

– ﬂux-limited diﬀusion—an improvement of the diﬀusion approximation,

but still inaccurate at small optical depths;

– escape probability methods—essentially assume that a photon after its

emission is either re-absorbed on the spot, or escape freely from the

medium; the branching ratio between these two is a function of optical

distance form the nearest boundary. Again, the radiation intensity is

eﬀectively eliminated from the state vector.

• Monte-Carlo methods

These are powerful methods, in particular in multi-D geometries, but

there are problems when treating NLTE problem (but see some contribu-

tions in this Volume for eﬃcient generalizations of Monte-Carlo methods

to NLTE situations).

• “Exact” numerical solutions. This is the most important class of solutions,

so the next Section will be fully devoted to it.

6 Exact Numerical Solution

In astrophysical applications, there are essentially two types of problems:

i) solution of a full set of structural equations, generally with a large number

of frequency points and a large number of atomic level populations to be

solved for—the problem is thus highly non-linear; and

ii) solution of a (generally anisotropic) scattering problem with a known

structure (temperature, density, and velocity are speciﬁed)—a linear problem.