Radiative transfer

A. Jerkstrand, U. Noebauer, C. Vogl + several others

February 2, 2018

1 The radiation field

The radiation field is in general a 7-dimensional field quantity, varying over 3 spatial dimensions,
2 angle dimensions, frequency, and time. If the spatial resolution is Ns , the angular Nµ , frequency
Nν , and time Nt , just storing the radiation field takes
 3  2   
11 Ns Nµ Nν Nt
N = 10 (1)
100 10 1000 1
i.e. 100 billion numbers, orders of magnitude more than what can be stored in a modern computer
(∼ 1 billion numbers).
The vast majority of RT work is therefore limited to spherical symmetry, where 2 space and 1
angle variables are dropped, reducing the dimensionality from 7 to 4. The simplest problems we can
conceive of drop also the time dependence (steady-state) and frequency variable (gray transport),
reducing dimensionality to 2.
If the chief objective is to obtain the dynamic effects of the radiation field, much simplification
and reduced dimensionality and resolution is often motivated. On the other hand, if the objective
is to predict spectra in high detail, larger problems need to be considered.

1.1 Spatial resolution Ns

In principle, the radiation field can change over length scales corresponding to the mean-free-path
(the average length a photon travels before interaction, see more below). Equivalently, an optical
depth of 1000 requires 1000 bins.
An immediate problem makes itself known; individual lines as well as bound-free continua can
reach enormous optical depths of order τ = 1010 or more. Clearly, such situations cannot be
resolved by any kind of finite difference scheme.

1.2 Frequency resolution Nν

The needed frequency resolution is determined by the frequency scale over which emission and
absorption coefficients vary. The limiting factor in photon transport becomes bound-bound tran-
sitions where the line widths are determined by the thermal broadening.
 1/2
p T
v= 2kT /m = 12 km s−1 Ā−1/2 (2)
104 K
A line at say 5000 Å then has a thermal width of λ0 v/c = 0.2 Å, and to resolve it we would need
frequency spacing of ≪ 0.2 Å. To model the whole optical spectrum (3000-10,000 Å), we would
then need ≫ 35, 000 frequency points.
For continuous emission and absorption, a significantly smaller number of frequency points
would be sufficient.

2 Radiation field quantities

As particles can be treated as ensembles in statistical mechanics and thermodynamics, photons
and other radiation particles can be treated in a macroscopic field framework. For this to work,
over any scale in length, angle, or frequency over which physical conditions change, the number

1
of photons must be ≫ 1. This is fulfilled in the vast majority of situations encountered, but
exceptions exist (e.g. single-photon heating of dust grains).
The most commonly used description of the radiation field in the macroscopic picture is the
specific intensity Iν (x, y, z, Θ, Φ, t). This describes the energy flux per unit area at position
x = (x, y, z) = (r, θ, φ), at frequency ν, into direction n = (Θ, Φ) = “ω ′′ (per steradian). In
general we will take Θ to be the angle relative to the radial direction, and Φ to be the orthonormal
(azimuthal) angle. The unit of Iν is erg cm−2 s−1 Hz−1 ster−1 . An important property of Iν is
that in the absence of sources or sinks, it stays constant along a ray.
The specific intensity related to the electromagnetic field E as (Kanschat 2009)
Z t+∆t/2
1
Iν = E 2 (x, t′ )dt′ (3)
∆t t−∆t/2

where a time-scale ∆t long enough to cover many oscillations of E but short enough to have
neglegible evolution of the amplitude is assumed.
An equivalent quantity is the photon number density ψν :
1
ψν = Iν (4)
chν
which has unit cm−3 Hz−1 ster−1 .
A third is the photon distribution function fR :
c2
fR = Iν (5)
h4 ν 3
which has units cm−4 g−1 s ster−1 and describes the number density in phase space (space and
momentum). A description of the 6-D phase space distribution is of course equivalent to the dis-
tribution over space, frequency, and angles. fR is a relativistically invariant distribution function.
3
This means that to transform Iν between frames, there is a (ν/ν ′ ) conversion factor.

2.1 Zeroth moments

The mean intensity Jν is the angle-average of Iν :
Z
1
Jν = Iν dω (6)
4π
which
R is also called the zeroth moment of the radiation field. It has the same units as Iν . Here
dω denotes integration over all angles.
An equivalent quantity is the monochromatic radiation energy density Eν :
Z
Eν ≡ hν ψν dω (7)

which has units erg cm−3 Hz−1 . Combining the last two equations gives Eν = 4π/cJν .

2.2 First moments

The monochromatic radiation flux Fν is a vector defined by
Z
Fν = Iν (ω)n(ω)dω (8)

which is also called the first moment of the radiation field. It has unit erg s−1 cm−2 Hz−1 , and
thus respresents the net flux of energy through a unit surface (so the luminosity at radius r is
4πr2 Fν (r) in spherical symmetry).
Sometimes an alternative first moment Hν , or the Eddington flux, is used:
1
Hν ≡ Fν (9)
4π
In spherical symmetry, Fν (and Hν ) are always in the radial direction and their magnitudes Fν and
Hν are used for notation, then called the flux. Note that by definition of the above moments the

2
constraint |Fν | ≤ cEν holds, which expresses that radiation energy cannot be transported faster
than with the speed of light.
A third equivalent quantity is the monochromatic radiation momentum density gν :

gν = c−2 Fν (10)

which has units erg s cm−4 Hz−1 .

2.3 Second moments

The (3 × 3) radiation stress tensor Pij
ν is
Z



Pij
ν = ψν (ω)hν ni · n(ω) nj · n(ω) dω (11)

and is also called the second moment of the radiation field. It represents the net rate of transport
of momentum in direction i across a surface with normal in direction j. It has dimension dyne
cm−2 Hz−1 , and is symmetric. From the definitions of Pij ν and Eν directly follows the general
relation X
Pii
ν = Eν (12)
i=1,2,3

between the trace of the second-moment tensor and the zeroth moment.
In spherical symmetry Pijν is a diagonal matrix. If we let the three direction axes be (θ̂, φ̂, r̂),
then
 
(Eν − Pν )/2 0 0
Pij
ν =
 0 (Eν − Pν )/2 0  (13)
0 0 Pν
where Z Z 1
2π
Pν ≡ Pνr̂r̂ = ψν (ω)hν (r̂ · n(ω)) (r̂ · n(ω)) dω = Iν µ2 dµ (14)
c −1

where we have introduced µ = r̂ · n. For isotropic radiation fields (Iν (µ) = Iν = Jν ), we get
Pν = 4π/3cJν = Eν /3.
Sometimes alternative quantities Kij
ν , Kν (called ’the K-integral’) are used:

c ij
Kij
ν = P (15)
4π ν
Z 1
c 1
Kν ≡ Pν = Iν µ2 dµ (last for spherical symmetry) (16)
4π 2 −1

The mean radiation pressure P¯ν is

1 X
P̄ν ≡ P ii = Eν /3 (17)
3 i=1,2,3 ν

Note that in general P̄ν equals Pν only if Pν = Eν (which results in the free-streaming case).

2.4 Other quantities

Sometimes the normalized moments are better suited to characterize the radiation field than
the ordinary moments. The normalized 1st moment Fν /(cEν ) = Hν /Jν is called flux factor and
it measures how much net transport of radiation energy takes place. Its absolute value is bounded
by 0 and 1. Values close to zero are indicative of conditions close to isotropy (in which nearly no net
energy is transported in any direction), while values close to 1 are resulting in the free-streaming
regime (in which all radiation particles move along the same direction and therefore transport all
their energy in that direction).
The variable Eddington factor fν is the ratio of the second and zeroth moments in the
radial direction:
Pν Kν
fν ≡ = (18)
Eν Jν

3
Its value must lie between 0 and 1. Similarly to the flux factor the Eddington factor measures the
degree of anisotropy of the radiation field, while going from isotropic diffusion to free-streaming
conditions, fν runs from 1/3 to 1 (while values smaller than 1/3 are possible). As a geometric
measure, it can usually be obtain to satisfactory accuracy by solving the time-independent transfer
equation (page 343).

Quantity General Radial component

Radiation field Iν , ψν , fR -
Zeroth moment Jν , Eν -
First moment Fν ,Hν , gν Fν , Hν , gν
Second moment Pijν ,Kν
ij
Pν , Kν

Table 1: Radiation field quantities.

Some methods deploy even higher order integrals of Iν , but most use some combination of those
described here.

3 Transfer equation quantities

To set up the transfer equation, we need terms specifying the creation and destruction of radiation;
these are the emission and absorption coefficients.
The emission coefficient (or emissivity) is denoted ην (x, y, z, Θ, Φ, t) and is in general a
7-D quantity, with unit erg s−1 cm−3 Hz−1 ster−1 . The product ην ds gives the amount of specific
intensity gained per unit cm. In the fluid rest frame (comoving frame), η is in general still 7-D as
scattering emissivity can be angle-dependent. If there is no scattering or if the scattering isotropic,
it reduces to an isotropic 5-D quantity.
The absorption coefficient is denoted χν (x, y, z, Θ, Φ, t) and is also a 7-D quantity, with unit
cm−1 . The product Iν χν gives the amount of specific intensity dIν lost per cm of path. The
absorption coefficient is also called opacity in Mihalas. In the fluid rest frame, χ is isotropic so
reduces to a 5-D quantity.
The mean-free path λν is
1
λν ≡ (19)
χν
and has unit cm; it is how far a photon typically travels before it interacts with the matter.
The optical depth between two points is the integral of the absorption coefficient between the
points Z
τν (x, x′ ) = χν (x + ns)ds (20)

Absorption can sometimes be divided into the processes of scattering and thermalization. By
scattering we mean that the photon is reemitted with the same energy in the fluid frame, but in a
new direction. This happens for example in Thomson scattering, and in resonance line scattering.
By thermalization we mean that the photon energy goes to thermal kinetic energy of the gas (and
possibly also to potential energy). This happens for instance when a line absorption is followed
by collisional deexcitation, or when photoionization (here also potential energy is created along
with thermal energy). One should note that in the general case, one does not know a-priori what
happens following a line absorption for instance; the atom can do a myriad of things. Thus, this
division between scattering and thermalization does not really exist in the most general problems
(but for many simplified ones used in practise).
Mihalas represents thermal absorption by κ and scattering by σ, so
χν = κν + σν (21)
the emissivity can also be divided into thermal and scattering components
ην = ηνt + ηνs (22)
The source function is the ratio between emissivity and absorption coefficient
ην
Sν = (23)
χν
it has the same units as Iν .

4
3.1 Scattering details
In general, the scattering emissivity function is (e.g. Eq 1-23 in Mihalas 1978)
Z ∞Z
s 1
ην = Iν ′ (dω ′ )R(ν ′ , n′ , ν, n)dν ′ dω ′ (24)
4π 0

where R is the redistribution function, a 10-D quantity (3 space, 2 incoming angle, 2 outgoing
angle, 1 incoming frequency, 1 outgoing frequency, 1 time). Its high dimensionality almost always
requires simplifying assumptions.

3.1.1 Complete coherence limit (page 13 in Mihalas 1978)

In this limit we assume that ν = ν ′ (in the comoving frame), so we have ’perfect scattering’. This
is usually the treatment for continuum scattering, and sometimes line scattering (Mihalas 1978).
Then
R = g(n′ , n)σν ψ(ν ′ )δ(ν − ν ′ ) (25)
where ψ is the normalized line profile function. If we also assume isotropy, then g = 1, and
Z
s 1
ην = σν Iν dω ′ = σν Jν (26)
4π
2
For dipole scattering g = 3/4(1 + cos (n · n′ ) ).

3.1.2 Complete incoherence (complete redistribution) limit (page 12 in Mihalas

1978)
In this limit we instead assume that there is no correlation at all between ν ′ and ν. Then..

4 The transfer equation

The transfer equation corresponds to the Boltzmann equation for photons. The non-relativistic
transfer equation in the observer frame is (Eq 76.3 in Mihalas & Mihalas 1984), along a path s
 
1 ∂ ∂
+ Iν (x, n, t) = ην (x, n, t) − χν (x, n, t)Iν (x, n, t) (27)
c ∂t ∂s

Depending on the problem/algorithm, ην and χν may be explicitly known, or they may be functions
of Iν (in which case the problem is much more difficult). In the first case we have to solve a PDE
in two variables (t, s), in the second case an integro-PDE equation.
The ∂/∂s differential corresponds in generic notation to n · ∇:

∂I ∂I ∂x ∂I
= + ... = nx + ... = (n · ∇) I (28)
∂s ∂x ∂s ∂x
In curvilinear coordinates, moving on a straight line path ds corresponds to continuous rotation
of the basis vectors (e.g., in spherical symmetry the angle to the normal Θ changes along a ray
unless Θ = 0.) In spherical symmetry,

∂ ∂r ∂ ∂Θ ∂
= + (29)
∂s ∂s ∂r ∂s ∂Θ
where we have used dr = cos Θds and rdΘ = − sin Θds.
Eq 27 then becomes
 
1 ∂ ∂ 1 − µ2 ∂
+µ + Iν (x, n, t) = ην (x, n, t) − χν (x, n, t)Iν (x, n, t) (30)
c ∂t ∂r r ∂µ

5
4.1 Zeroth moment of the transfer equation
Integrating the transfer equation over angle, we get (Eq. 78-2 in M84)
Z
1 ∂Jν 1
+ ∇ · Hν = ην − χν Iν dω (31)
c ∂t 4π
Here the n in the second term has moved inside the integral, giving a first-moment of the radiation
field term (H). R
In the comoving frame, χν is isotropic, which means that the RHS will contain only Iν dω = Jν ,
not Iν .
By substituting the expression relating Jν to Eν , one can show that the first moment equation
expresses energy conservation : the time-change in energy density is the net creation minus the
net outflow (see Eq 78.3 in Mihalas & Mihalas 1984).
Radiative equilibrium means that the gas emits exactly as much as it absorbs
Z Z
ην − χν Iν dω = 0 (32)
ν ω

This must hold in any static situation. Since such a situation also has ∂Jν /∂t = 0, it follows from
Eq. 31 that ∇ · H = 0, which gives L = constant. Thus, in a stellar atmosphere for example, the
luminosity is the same at each depth. In the interior of the star radiative equilibrium does not
hold because of energy generation by fusion, and a gradient in L is established.
Discuss situations when radiative equilibrium can be assumed. In SNe if gamma rays from
radioactivity are included in the radiation budget?

4.2 First moment of the transfer equation

Multiply the transfer equation 27 by n and integrate over angle. This gives
 Z Z
1 ∂ 1
+n·∇ Iν n(ω)dω = (ην − χν Iν ) n(ω)dω (33)
c ∂t 4π

This gives (derive this)

Z Z
1 ∂Hν 2 1
+∇· Iν n dω = (ην − χν Iν ) n(ω)dω (34)
c ∂t 4π

Or (Mihalas & Mihalas (1984, Eq 78.9))

Z
1 ∂Hν 1
+ ∇ · Kij
ν = (ην − χν Iν ) n(ω)dω (35)
c ∂t 4π
This equation expresses momentum conservation of the radiation field. The time rate of change
of momentum is the transport of momentum across the boundary surface (term 2 on LHS), plus
addition of momentum by emission (term 1 on RHS), minus removal of momentum by absorption.
Then, the radiative force due to absorption must equal this last term
Z Z
1
fR = χν Iν n(ω)dω (36)
c
Derive Eddington luminosity from this expression.
In spherical symmetry
Z
1 ∂Fν ∂Pν 3Pν − Eν 2π
+ + = (ην − χν Iν ) µdµ (37)
c2 ∂t ∂r r c
Again, in the comoving frame the RHS will contain not Iν but Hν .

6
4.3 Closure
The moment equations, at least in the comoving frame (or static case), do not contain Iν , but only
angle-integrated quantities. Thus, the dimensionality of the problem is reduced by 2 in 3D/2D
problems, and 1 in spherical symmetry. This is the main motivation for seeking solutions to them
instead of to the full transfer equation.
However, we always have at one unknown moment (possibly non-scalar) more than the number
of equations. This is called the closure problem. Of course, we cannot arbitrarily simplify our
problem just by integrating the equations. All the angular details of Iν need still to be known to
define an exact closure relation.
A closure relation can be obtained by variable Eddington factor methods (often used in
non-hydro), where an iteration with a formal solver of the normal transfer equation (with source
function fixed) is done. Note there is some subtility here : if we anyway have to solve the transfer
equation many times, why don’t we just take J etc from its solution directly? The answer is that
convergence is faster if we determine J from the moment equations (more later).
Even more simply some codes use analytic formulae from the geometry (often used in radi-
ation hydro).
Note that in multi-D, we have 4 equations but 13 (Jν ,Fν ,Pij ν ) unknowns. Thus we need not
one closure relation but 9 in each cell. In spherical symmetry we have 2 equations and 3 unknowns,
so a scalar closure is sufficient.

4.4 Time-dependence or not

A central question for many problems is whether the time-derivative term is needed or not. A
good discussion of this is given in Mihalas & Mihalas (Section 6.5 in 1984).
Define a optically thin radiation flow timescale tR :
l
tR = (38)
c
where l is some length scale.
If the distance l is optically thick (l/λ = τ ≫ 1), it takes instead a diffusion time scale td

l2
td = = τ tR (39)
cλ
for a photon to cross the distance l. (τ 2 = l2 /λ2 scatterings, each taking a time λ/c).
Define also a fluid flow timescale
l
tf = (40)
v
where v is a fluid velocity scale.
Regimes for the transfer:
1. (Effectively) optically thin regime. If the gas is optically thin, tR /tf = v/c.
Mihalas & Mihalas (1984): "If v/c ≪ 1, the radiation field at any position adjusts essentially
instantaneously to changes in physical conditions”.
For example, if the radiation field changes because the density is getting lower in an optically
thin outflow with v/c ≪ 1, we can drop the time-derivative.
If we care about the effect of radiation on the fluid flow, we need to resolve the fluid flow
timescale. If v/c ≪ 1, tR /tf ≪ 1, and a change in fluid conditions changes the radiation field
on a time-scale much shorter than the fluid flow time scale.
Another way to see this is to write (?)
1 ∂ ∂ 1 ∂ ∂
+ = + (41)
c ∂t ∂s c ∂t v∂t
So the first and second terms have relative size v/c.
Note that this line of reasoning assumes that changes in emission and absorption are related
to changes in flow, i.e. occur on time scale tf .
Mihalas & Mihalas (1984) defines this regime as quasi-static.

7
 Note that is the region is optically thick, but τ is small enough and/or v small enough that
td /tf = τ v/c < 1 (i.e. 1 ≪ τ ≪ c/v). Then we come to the same regime as the optically
thin case.
2. Optically thick, frozen diffusion regime. Now assume the region is optically thick and
td /tf ≫ 1. Then the formal diffusion time is long compared to the flow time, and it seems
we would need to solve the diffusion time-dependently.
However, thermalization can save the day here. If the optical depth is high, for many sit-
uations the radiation will interact with the matter and that ’resets’ conditions. Then, the
’effective’ time that photons travel without being destroyed is not td but something smaller.
3. Optically thick, dynamic diffusion regime. Assume τ is large enough that td /tf =
τ v/c > 1, but not so large as to td /tf ≫ 1, i.e. td ∼ tf .
Mihalas & Mihalas (1984) states that here we must include the time-derivative term, and
solve on the fluid time scale.
4. External power source. If some time-varying external power source, such as a neutrino
source or radioactivity, governs the radiation field, the radiation field can change on the time-
scale ts of that source, unrelated to any flow time scale tf . Thus, if we take tR = R/c = v/ct,
we need tR ≪ ts for a time-independent solution, or
c
t≪ ts (42)
v
For the optically thick case dynamic diffusion situation, we have to solve time-dependently
on the shorter time-scale of ts and tf .
Discussion : come up with examples and identify regimes Explicit vs implicit solutions?

5 Reference Frames
Photons propagate with the speed of light. Consequently, special relativistic effects always have
to be considered to a certain extent when addressing radiative transfer problems. This implies
particularly that one has to think about the proper reference frame in which the radiation field is
described and its evolution solved.
There are three frames which are of fundamental importance in radiative transfer theory (e.g.
Mihalas1986, page 6)
1. Atom frame. Frame in which the individual atom undergoing an interaction has zero
velocity.
2. Comoving frame (CMF). Frame at (x,t) in which the integral of velocity vectors in a
small volume around x is zero. Identifical to the Lagrangian frame. Since the matter is at
rest in this frame, the opacity and often also the emissivity are isotropic (they are, however,
in general still complex functions of frequency).
3. Observer frame (lab-frame, LF). Frame in which observer is at rest.
When considering LF and CMF it is important to realise that while the LF is an inertial frame,
the CMF is in general not. Since the ambient material moves at a velocity which is a function
of time and space, the CMF cannot be defined globally. Instead, CMF should be understood as
a sequence of inertial frames each moving with the instantaneous velocity of the fluid element
under consideration. Lorentz transformations can then be used to transform between these inertial
frames and the LF1 .
Combining energy and momentum gives a four-vector and Lorentz transforming it yields two
important special relativistic effects in radiative transfer, namely the Doppler shift

ν0 = γν (1 − n · v/c) , (43)
ν = γν0 (1 + n0 · v/c) (44)
1 Mihalas & Mihalas (1984) points out that if one applies this concept rigorously, a self-consistent description is

obtained which yields results in agreement with experiments.

8
and the aberration
 
γn · v/c
n0 = (ν/ν0 ) (n − γv/c) 1 − , (45)
γ+1
 
γn0 · v/c
n = (ν0 /ν) (n0 + γv/c) 1 + . (46)
γ+1

Here, we have adopted the nomenclature of Mihalas & Mihalas (1984) and denoted quantities eval-
uated in the CMF with a subscribed ’0’. In one dimensional geometries, the aberration simplifies
to
µ − v/c
µ0 = , (47)
1 − v/cµ
µ0 + v/c
µ= . (48)
1 + v/cµ0

The consequences of this aberration effect are illustrated in Figure 1.

Figure 1: Effect of angle aberration. The polar plot shows the probability density of finding a
photon propagating in the direction µ in the LF under the assumption that the radiation field was
isotropic in the CMF, i.e. ρµ0 (µ0 ) = 1/2. Discuss in context of GRBs.

Using a series of simple “Gedankenexperiments”, Thomas (1930) derived the transformation

laws for the specific intensity
 3
ν
I(ν, µ) = I0 (ν0 , µ0 ) (49)
ν0
and the material functions emissivity
 2
ν
η(µ, ν) = η0 (ν) (50)
ν0

and opacity
 −1
ν
χ(ν, µ) = χ0 (ν0 ) . (51)
ν0

5.1 Moments
TBD..

6 Transfer equation : moving medium

When the fluid is accelerating (as during the SN explosion) or there is velocity gradient (as in
the coasting phase), it is difficult to solve the transfer equation with all quantities in the observer
frame, because ηn u and χν are anisotropic.

6.1 Observer frame equation with comoving η and χ from first-order

expansions
The attraction of retaining the observer frame is that the transfer equation is simpler than in the
comoving frame. Described in section 93 in Mihalas & Mihalas (1984).
The limitation of this method is when lines are important; the first-order expansions are then
insufficient.

9
6.2 Comoving frame equation
Described in section 95 in Mihalas & Mihalas (1984). The attraction of the comoving frame is the
isotropy of χ (and η) and easier matter-radiation interaction calculations. The drawback is a more
complex transfer equation. The equation is now not only a PDE in space and time, but also over
frequency and angle.
Discuss isotropy of η when scattering.
The fully relativistic comoving transfer equation even in spherical symmetry is a monster ex-
pression (Eq 95.9 in Mihalas & Mihalas 1984). The moment equations are also long (95.11 and
95.12).
If we retain only v/c terms, and also ignore the fluid acceleration term a, we get (Eq 95.17)
1 DI0
(52)
c Dt
µ ∂  2 
+ 2 r I0 (53)
   r ∂r 
∂ 1 µ0 v ∂v
+ (1 − µ20 ) + − I0
∂µ0 r c r ∂r
   
∂ 2 v µ20 ∂v
− ν0 (1 − µ0 ) + I0 (54)
∂ν0 cr c ∂r
 
2 v (1 + µ20 ) ∂v
+ (3 − µ0 ) + I0 (55)
cr c ∂r
= η0 − χ0 I0 (56)
Compare CMFGENs equation (Dessart Hiller 2012), which uses v/c ≪ 1 and homology (dv/dr =
v/r:
1 ∂I0 µ0 c + v ∂I0 (1 − µ20 ) ∂I0 vν0 ∂I0 3v
+ + − + I0 = η0 − χ0 I0 (57)
c ∂t c ∂r r ∂µ0 rc ∂ν0 rc
where D/Dt = ∂/∂t + v∂/∂r

7 Analytic solutions
7.1 The Time-Independent Formal Solution
Described in section 79 of Mihalas & Mihalas (1984). We restrict our discussion to plane-parallel
geometries. The time-independent radiative transfer equation
∂Iν
µ = Iν − Sν (58)
∂τν
can be rewritten as follows
∂[Iν exp (−τν /µ)] Sν exp (−τν /µ)
=− . (59)
∂τν µ
If the source function Sν is known, the radiative transfer problem can be solved by integration of
Eq. (59) Z τ2
I(τ1 , µ, ν) = I(τ2 , µ, ν)e−(τ2 −τ1 )/µ + µ−1 Sν (t)e−(t−τ1 )/µ dt. (60)
τ1

Eq. (60) is the so-called formal solution of the equation of transfer. For the outgoing intensity
(µ ≥ 0) in a semi-infinite medium τ1 = τν and τ2 = ∞. The formal solution is then given by:
Z ∞
I(τν ; µ, ν) = Sν (t)e−(t−τν )/µ dt/µ, (0 ≤ µ ≤ 1) (61)
τν

If we assume that no radiation is entering through the outer boundary (I − (0) = 0), we obtain
Z τν
I(τν ; µ, ν) = Sν (t)e−(τν −t)/(−µ) dt/(−µ), (−1 ≤ µ ≤ 0). (62)
0

for the incoming radiation field (µ ≥ 0).

10
Schwarzschild-Milne Equations
Using Eqs. (61) and (62) we can derive the following expression for the mean intensity
Z 1 Z ∞ Z 0 Z τν 
1
Jν (τν ) = dµ Sν (t)e−(t−τν )/µ dt/µ + dµSν (t) e−(τν −t)/(−µ) dt/(−µ) . (63)
2 0 τν −1 0

To simplify the problem, we substitute w = ±1/µ and change the order of integration:
Z ∞ Z ∞ Z τν Z ∞ 
1 e−w(t−τν ) e−w(τν −t)
Jν (τν ) = dtSν (t) dw + dtSν (t) dw . (64)
2 τν 1 w 0 1 w

By identifying the integrals that appear in Eq. (64) as the first exponential integral, we arrive at
a concise expression for the mean intensity
Z
1 ∞
Jν (τν ) = Sν (tν )E1 |tν − τν | dtν . (65)
2 0

The exponential integral is defined as follows

Z ∞ Z ∞
−n −xt n−1
En (x) ≡ t e dt = x t−n e−t dt. (66)
1 x

Similarly, expressions for the next two moments can be derived:

Z ∞ Z τν
Fν (τν ) = 2π Sν (tν )E2 (tν − τν ) dtν − 2π Sν E2 (τν − tν ) (67)
τν 0
Z ∞
1
Kν (τν ) = Sν (tν )E3 |tν − τν | dtν (68)
2 0

Due to the central importance of Eq. (65) in radiative transfer theory an abbreviated operator
notation has been introduced. Following this notation the mean intensity

Jν (τν ) = Λτν [S(τν )] (69)

is obtained by applying the so-called lambda operator

Z
1 ∞
Λτ [f (t)] ≡ f (t)E1 |t − τ | dt (70)
2 0

to the source function. Over time, however, the term ’lambda operator’ has taken on a much
broader meaning, being used to describe any procedure (including non-analytic) to obtain J from
S.

7.2 Wave Limit

In a vacuum (χν = ην = 0) the transfer equation (27) reduces to
 
1 ∂ ∂
+ Iν (x, n, t) = 0. (71)
c ∂t ∂s

Introducing
I + ≡ Iν (x, n, t) (72)
and
I − ≡ Iν (x, −n, t) (73)
we obtain the transfer equations
∂I + ∂I +
+c =0 (74)
∂t ∂s
and
∂I − ∂I −
−c =0 (75)
∂t ∂s

11
 from Eq. (71). Defining the mean-intensity-like quantity
1 +
j≡ (I + I − ) (76)
2
and the flux-like quantity
1 +
h≡(I − I − ) (77)
2
addition and subtraction of Eqs. (74) and (75) yields

∂j ∂h
+c =0 (78)
∂t ∂s
and
∂h ∂j
+c = 0. (79)
∂t ∂s
The quantities h and j will appear again later in the lecture when we will derive the Feautrier
Equations. Taking the partial derivative of Eq. (78) with respect to time and inserting Eq. (79) in
the resulting equation yields
∂2j 2
2∂ j
= c . (80)
∂t2 ∂s2
Similarly
∂2h ∂2h
2
= c2 2 (81)
∂t ∂s
can be obtained. Eqs. (80) and (81) are wave equations with the solutions

j(s, t) = A1 f1 (s − ct) + A2 f2 (s + ct) (82)

and
h(s, t) = B1 f1 (s − ct) + B2 f2 (s + ct). (83)
Here A1 , A2 , B1 , B2 are constants that are determined by the initial and boundary conditions.
One possible solution of Eqs. (82) and (83) is a monochromatic, plane wave:

I(x, t; n′ , ν ′ ) = I0 δ(s − ct)δ(n′ − n)δ(ν ′ − ν) (84)

In this case Jν = Hν = Kν and the Eddington factor fν = 1.

7.3 Diffusion limit

Advanced Reading
The discussion below assumes a static medium and that the radiation field and the medium are in
thermal equilibrium. Radiative diffusion in moving media and in nonequilibrium are discussed in
section 97 of Mihalas & Mihalas (1984).

Static, LTE Atmospheres

Described in section 80 of Mihalas & Mihalas (1984). For τν ≫ 1 the source function Sν (see
Eq. (23)) approaches the Planck function Bν . It is thus possible to write the source function at
any optical depth tν as a Taylor expansion around some reference optical depth τν :
X∞
∂ n Bν
Sν (tν ) = n
(tν − τν )n /n! (85)
n=0
∂τ ν

Recall Eq. (61) for the intensity of outgoing radiation in a planar, static medium - Thomas Jankas
effective absorption
Z ∞
′
I(τν ; µ, ν) = Sν (τν′ )e−(τν −τν )/µ dτν′ /µ, (0 ≤ µ ≤ 1). (86)
τν

12
The assumption of plane-parallel geometry is justified due to the small photon mean free paths.
Inserting the Taylor expansion of the source function into Eq. (86) yields
∂Bν (τν ) ∂ 2 Bν (τν )
I(τν , ν) = Bν (τν ) + µ + µ2 + ... (87)
∂τν ∂τν2
For −1 ≤ µ ≤ 0 the intensity can be calculated in a similar fashion using Eq. (62). The result is
identical to Eq. (87) apart from terms of the order of exp(−τ /µ), which vanish for τ → ∞. From
the defining Eqs. (6), (9) and (16) we find
1 ∂ 2 Bν (τν )
Jν = Bν (τν ) + + ... (88)
3 ∂τν2
1 ∂Bν (τν ) 1 ∂ 3 Bν (τν )
Hν = + + ... (89)
3 ∂τν 5 ∂τν3
1 1 ∂ 2 Bν (τν )
Kν = Bν (τν ) + + ... (90)
3 5 ∂τν2
for the moments of the specific intensity. Replacing derivatives by difference quotients
(i.e. ∂Bνn (τν )/∂τνn → Bν /τνn ), we see that the ratio of successive terms is of O(1/τν2 ) = O(λ2ν /l2 ).
Here λν denotes the photon mean free path and l is a characteristic structural length (e.g. a
pressure scale height in a stellar envelope.) Since the ratio λν /l is small (typical values in the sun
range from 10−7 to 10−10 ) it is sufficient to retain only the first terms in Eqs. (88) to (90):
Jν (τν ) = 3Kν (τν ) = Bν (τν ) (91)
Thus both the mean intensity Jν and the radiation pressure Kν have their equilibrium values. In
contrast to equilibrium the flux
1 ∂Bν (τν ) 1 ∂Bν dT
Hν = =− (92)
3 ∂τν 3χν ∂T dr
is nonzero. By integrating this expression over frequency, we obtain the total flux
Z ∞ 
1 ∂Bν (τν ) dT
F = −(4π/3) dν . (93)
0 χ ν ∂T dr
This equation is formally identical to Fourier’s Law for heat conduction. Introducing the so called
Rosseland mean opacity
Z ∞ Z ∞
−1 1 ∂Bν (τν ) ∂Bν (τν )
χR = dν (94)
0 χν ∂T 0 ∂T
we can define a radiative conductivity
4π dB 4
KR = = cλR aR T 3 (95)
3χR dT 3
in analogy to the thermal conductivity.
In conclusion, at large optical depths the transfer problem can be described by the single equation
Eq. (92), which behaves like a diffusion equation. The dimensionality of the problem has been
reduced from six to one dimensions.

7.4 The Grey Atmosphere

For a grey material the opacity is independent of frequency i.e. χν = χ. As a consequence the
radiation field becomes independent of the state of the material.

Applications
• Starting point in the calculation of more complex models. For stellar atmospheres it is typical
to proceed through a series of intermediate models with increasing physical complexity e.g.
LTE-gray → LTE → NLTE
• Test problem for numerical methods
• Provides boundary conditions for stellar structure calculations
• Neutron transport in heavy-water nuclear reactors

13
Basic Results
Frequency integration of the time independent, planar transfer equation yields:
∂I
µ =I −S (96)
∂τ
R∞ R∞
Here I ≡ 0 Iν dν and S ≡ 0 Sν dν denote the total intensity and source function respectively.
For a grey material the radiative equilibrium condition
Z ∞ Z ∞
χν Jν dν = χν Sν dν (97)
0 0

reduces to the simple requirement J = S. Using this result Eq. (96) can be simplified as follows:
∂I
µ =I −J (98)
∂τ
This is the transfer equation for a plane-parallel, grey atmosphere in radiative equilibrium. From
the first moment of the transfer equation
dH
=J −S =0 (99)
dτ
we see that the flux is constant throughout the atmosphere. This is a general result for static,
planar atmospheres in radiative equilibrium. The second moment equation is given by
dK
= H, (100)
dτ
which has the solution
1
K(τ ) = Hτ + c = F τ + c. (101)
4π
Recall that for large optical depths the specific intensity can be approximated as I(µ) = I0 + I1 µ,
which implies that K(τ ) = 1/3J(τ ). From this and Eq. (101) we infer that
3
J(τ ) → Fτ (τ ≫ 1). (102)
4π
The general solution can then be written as
3
J(τ ) = F (τ + q(τ )), (103)
4π
where q(τ ) denotes the as of yet undetermined Hopf function. It is possible to connect the constant
c in Eq. (105) for the second moment of the radiation field to the newly introduced Hopf function
by taking the limit of large optical depths:
 
1 1
lim J(τ ) − K(τ ) = F lim [τ + q(τ ) − τ − c] = 0 (104)
τ →∞ 3 4π τ →∞
Thus c = q(∞) and the second moment is given by
1
K(τ ) = Hτ + c = F τ + q(∞). (105)
4π

LTE
If LTE is assumed (i.e. Sν = Bν ) it is possible to assign a temperature T to the radiation field via
the radiative equilibrium equation:

J(τ ) = S(τ ) = B[T (τ )] = σR T 4 /π (106)

Defining the effective temperature Teff as the temperature a black body would have to reproduce
the emergent flux (i.e F = σR Teff
4
), we can rewrite Eq. (103) in terms of T
3 4
T4 = T [τ + q(τ )] (107)
4 eff

14
Mean Opacities
See section 3-2 in MH78. The goal is to define mean opacities in such a way that the general
equations
µ(∂Iν /∂z) = χν (Sν − Iν ) (108)
∂Hν /∂z = χν (Sν − Jν ) (109)
∂Kν /∂z = −χν Hν (110)
can be reduced to their grey counterparts

µ(∂I/∂z) = χ(S − I) (111)

dH/ dz = 0 (112)
dK/ dz = −χH. (113)
It is impossible to achieve a complete correspondence between the grey and nongrey problem.
However, suitable choices of the mean opacity can establish one to one correspondences for selected
quantities.

Rosseland Means
If the goal is to reproduce the correct integrated energy flux H, the mean opacity must be defined
as follows: Z ∞ Z ∞
− χ−1
ν (∂K ν /∂z) dν = Hν dν = H = −χ̄−1 (dK/ dz). (114)
0 0
Z ∞ Z ∞
χ̄ −1
= χ−1
ν (∂Kν /∂z) dν (∂Kν /∂z) dν (115)
0 0

Since Kν is not known a priori, it is necessary to find approximations in order to evaluate

the opacity. At high optical depths Kν → 1/3 Jν , Jν → Bν and we can write ∂Kν /∂z =
1/3 (∂Bν /∂T )(dT / dz). With these simplifications Eq. (115) can be written as

R ∞ 1 ∂Bν
1 0 χ ∂T
dν
= R ∞ ν∂Bν . (116)
χ̄R 0 ∂T dν

The assumptions used to derive the Rosseland mean are the same as those used in the diffusion
approximation. It is thus appropriate to use Rosseland means to describe radiative diffusion at
high optical depths. This allows the determination of the thermal structure of the atmosphere at
great depths via
3 4
T 4 = Teff (τ̄R + q(τ̄R )). (117)
4

Flux-Weighted Mean
To transform the nongrey equation for the second moment Kν [Eq. (110)] into the grey equation
(113) the mean opacity must be defined as follows
Z ∞
χ̄F ≡ H −1 χν Hν dν. (118)
0

This opacity is the flux-weighted mean of the frequency dependent opacity χν . We can verify that
the definition in Eq. (118) has the desired properties by integrating Eq. (110) over frequency:
Z ∞
− (dK/ dz) = χν Hν dν = χ̄F H (119)
0

Thus K(τ̄ ) = H τ̄ + c applies as in the grey case. This guarantees that the correct values for the
radiation pressure and radiation force are recovered. This is of relevance for the calculation of the
density structure of the atmospheres of early-type stars.

15
Planck and Absorption Means
Whereas the Rosseland mean is the appropriate average for optically thick systems the Planck
mean is suitable for optically thin systems. For details see section 3-2 in MH78.

Eddington approximation
We know that at high optical depths J = 3K holds. In the Eddington approximation this condition
is applied throughout the entire atmosphere. In combination with K = 1/(4π) · F τ + c this
simplifying assumption leads to the following expression for the mean intensity
3
JE = F τ + c′ . (120)
4π
Recall the formal solution for the flux at optical depth τν
Z ∞ Z τν
Fν (τν ) = 2π Sν (tν )E2 (tν − τν ) dtν − 2π Sν E2 (τν − tν ). (121)
τν 0

The flux at the outer boundary is then given by

Z ∞   
3 3 4
F (0) = 2π ′ ′
F τ + c E2 (τ ) dτ = 2πc E3 (0) + F − 2E4 (0) . (122)
0 4π 4 3
From F (0) = F follows that
3 E4 (0)
c′ = F. (123)
4π E3 (0)
Using the relation En (0) = 1/(n−1) we find that c′ = F/(2π). This implies that the mean intensity
is given by  
3 2
Je = F τ+ . (124)
4π 3
Thus the Hopf function in the Eddington approximation is q(τ ) = 2/3. In LTE the temperature
structure is determined by  
4 3 4 2
T = Teff τ + . (125)
4 3
Since T = Teff for τ = 2/3, this optical depth is commonly identified as the effective depth
of continuum formation. Despite the simplyfing Eddington approximation Eq. (125) provides a
quite accurate description of the thermal structure of a grey atmosphere. We expect the greatest
departures from the analytic solution to occur close to the boundary. However, the ratio of the
boundary temperature T0 to the effective temperature in the Eddington approximation T0 /Teff =
0.841 still agrees fairly well with the analytic solution T0 /Teff = 0.8114.

Solution with discrete ordinates

For an in-depth discussion see section 3-2 of MH78. Both approximate and exact solution can be
obtained by replacing the integrals in the transfer equation
Z
1 1
µ[∂I(τ, µ)/∂µ] = I(τ, µ) − I(τ, µ) dµ (126)
2 −1
with quadrature sums i.e.
Z n
1 1
1 X
I(τ, µ) dµ ≈ aj Ij (µj ) (127)
2 −1 2 j=−n
This reduces the integro-differential equation to a system of 2n coupled ordinary differential equa-
tions. In the limit of n → ∞ this approximation becomes exact and allows the derivation of the
analytic solution.

Spherical Geometry
A discussion of grey spherical atmospheres in LTE and radiative equilibrium can be found in section
7-6 of Mihalas (1978).

16
8 The Elementary Supernova Model – Formal Solution for
Line Transport in Moving Media
Having presented the formal solution to the radiative transfer equation, this concept can be used
to obtain a basic model for line formation in supernova ejecta. For this, a number of important
physical concepts have to be introduced. We first refresh the definition and implications of the two
important references frames in radiative transfer (see Section 5). Afterwards, we introduce atomic
line interactions using the two level atom and introduce with the Sobolev theory an important
approximation for treating such interactions in expanding flows. Finally, we set up the elementary
supernova model and use it to qualitatively understand line formation in supernova ejecta during
the photospheric phase.

8.1 Atomic Line Interactions: The Two-Level Atom

To determine the opacity and the source function appropriate to describe atomic line interactions,
we consider a very simple and idealised situation, namely that of a medium composed of “two-
level atoms”. Such a two-level atom has two isolated energy levels, lower l and upper u, with the
statistical weights gl and gu which are energy separated by
Eu − El = hνlu . (128)
Electrons can transition between these two levels by either absorbing or emitting a photon, be it
spontaneous or stimulated (see illustration in Figure 2). Einstein introduced transition probabilities

nu , g u Eu

hνlu = Eu − El

nl , g l El

Figure 2: Illustration of the two-level atom and the three fundamental processes which can occur.
Clockwise from the top right: absorption, spontaneous emission and stimulated emission.

for these processes and derived fundamental relations between them, which thus bear his name
(“Einstein relations”). The rate of energy absorbed via absorption processes in the line transition
per unit volume is given by
Ȧν = (Blu hνlu /4π)nl φν Iν , (129)
introducing the first Einstein coefficient Blu . Analogously, the second Einstein coefficient can be
defined via the rate of energy emitted per unit volume due to spontaneous emissions
Ėνspont = (Aul hνlu /4π)nu φν . (130)
The last Einstein coefficient describes stimulated emissions. In this process, an incoming photon
induces the emission of another photon with the same frequency and propagating in the same
direction as the incident quantum. Consequently, the specific intensity of the radiation field appears
in the expression of the rate of energy emitted by stimulated emission processes per unit volume
Ėνstim = (Bul hνlu /4π)nu φν Iν (131)
and ensures that the emitted radiation field has the same angular distribution as the incident one.
Stimulated emission thus has a functional dependence similar to Equation (129) and can be treated
as “negative absorption”.2
2 Assuming complete redistribution, i.e. that the emission line profile is identical to the absorption one.

17
 In thermodynamic equilibrium, one can easily relate these energy flow rates since detailed
balance must hold
Ȧν − E˙νstim = Ėνspont . (132)
In thermodynamic equilibrium, Iν = Bν and
nu gu
= exp(−hνlu /kT ) (133)
nl gl

holds. After rearranging Equation (132) into

Aul 1
Bν = gl Blu
(134)
Bul gu Bul exp(hνlu /kT ) − 1

one can easily identify the following relationships between the Einstein coefficients

gl Blu = gu Bul , (135)

3
Aul = (2hνlu /c2 )Bul , (136)

since
2hν 3 1
Bν = 2
. (137)
c exp(hν/kT ) − 1
With these, the line opacity, corrected for stimulated emission can be derived
 
g l nu
χν = nl (Blu hνlu /4π) 1 − φν . (138)
g u nl

The transition probabilities have to be determined experimentally (or by detailed quantum me-
chanical calculations). Often the line opacity is expressed in terms of the oscillator strength fl (see
e.g. Rybicki & Lightman 1979)
 
πe2 g l nu
χν = f l nl 1 − φν , (139)
me c g u nl

which is related to the Einstein coefficient via

4π 2 e2
Blu = fl . (140)
hνul me c
At the end, we note that the Einstein relations do not rest on the assumptions of thermodynamic
equilibrium or complete redistribution, but hold in general.

8.2 Line Transfer in expanding media - Sobolev Theory

As a photon moves along its trajectory through its environment, its frequency in the local rest
frame, i.e. in the CMF will continuously change, depending on the instantaneous velocity of the
material relative to the propagation direction. Thus, the photon can continuously shift into reso-
nance with an atomic line transition, whenever it overlaps with the non-negligible part of the line
profile function φ(ν). This behaviour makes a treatment of line interactions in general very chal-
lenging. Nevertheless, in rapidly expanding media, a number of simplifications can often be safely
applied, which significantly reduce the complexity of the problem. The formal theory underlying
these approximations has been first developed by Sobolev (1960)3 and is thus referred to as Sobolev
theory or Sobolev approximation.
Here, we refrain from a formal discussion of this field, which may be found e.g. in Castor
(1970); Rybicki & Hummer (1978); Hummer & Rybicki (1985), but instead provide a heuristic
but illustrative introduction to the field, following mostly Lamers & Cassinelli (1999). In this
illustration we assume a spherically symmetric environment and only consider one line. Naturally,
the Sobolev theory does not rest on these simplifications but can be derived for the fully three
dimensional case and be easily generalized to account for multiple lines. When working with photon
propagation in a spherical environment, it is often convenient to work with the so-called impact
geometry. To do so, we introduce the coordinate z, which is measured along the axis from the

18
 p

θ
r r θ

Figure 3: Illustration of the impact geometry coordinate system.

center of the spherically symmetric environment to the observer. Perpendicular to this axis, the
impact parameter p is measured. This construction is illustrated in Figure 3. Obviously,
z
µ=
, (141)
p r
r = p2 + z 2 , (142)

relate µ, r, z and p.
Consider a photon starting its propagation in the z-direction at the point (z1 , p) with the LF
frequency ν. The line optical depth it will encounter along this path is obtained by integrating
over the line opacity Z ∞
τ= κν (z)ρ(z)dz. (143)
z1

Here, we have introduced the specific line interaction cross section κν , which has dimensions of
cm2 g−1 . As discussed previously (see Section 8.1), the frequency dependent line opacity after
accounting for stimulated emission is
 
πe2 nu g l
χν = κ ν ρ = f l nl 1 − φ(ν). (144)
me c nl g u

In general, performing the integral is non-trivial, since the photon changes its CMF frequency
along its trajectory and can shift into and out of resonance with the line:
Z ∆ν(∞)  
πe2 nu g l dz
τ= fl nl 1 − φ(∆ν)d∆ν. (145)
me c ∆ν(z1 ) nl gu d∆ν

Here, we have used the variable ∆ν which describes the distance of the photon’s current CMF
frequency from the line center, ν0 :
 
v(r) z
∆ν(z) = ν 1 − . (146)
r c

We now take a closer look at the line profile function and investigate in more detail under which
conditions photons can actually interact with atomic lines. In general, the line profile function in
astrophysical applications is typically best described by a Voigt profile, consisting of Lorentzian
wings due to natural and collisional broadening and a Gaussian core, which accounts for Doppler
broadening due to thermal and microturbulent motions. For simplicity, however, we will only
consider the Gaussian core in the following discussion4
 
1 1 ∆ν 2
φ(∆ν) = √ exp − 2 . (147)
π ∆νG ∆νG
3 The Russian original was published already in√1947.
4 The characteristic width of a Voigt profile is 2 times the standard deviation of the Gaussian core.

19
Dicuss possibilities of microturbulence. As already mentioned, the mean velocities of the thermal
and microturbulent motion dictate the width of the Gaussian profile
r
ν0 2
∆νG = 2 i + hv 2
(hvth turb i). (148)
c 3
At ∆ν = ±1.5∆νG , the line interaction probability has already dropped to about 10 per cent of the
peak value at the line centre. Thus, we can define an effective line interaction region by requiring
that the local CMF frequency of the photon has to lie within
 
3 v(r) z 3
ν0 − ∆νG ≤ ν 1 − ≤ ν0 + ∆νG . (149)
2 r c 2
The key insight of the Sobolev approximation is to relate the extent of the line interaction region
to the typical length scales on which the plasma state of the ambient material changes significantly.
In order to benefit from the simplifications of the Sobolev approximation it is essential that the
radial velocity gradient is monotonous throughout the environment
sgn(dv/dr) = const. (150)
The concept of interaction regions can be formalised by introducing the so-called Sobolev length
vG
L(r, µ) = (151)
d(µv(r))/dz
At any point, these vectors span a surface which encompasses the region in which a photon emitted
at r by the line into the direction µ can be reabsorbed in the same line transition. The line
interaction region defined by (149) is then equal to three Sobolev lengths. If the plasma state
does not vary throughout this region, the line profile in the optical depth integration can be safely
replaced by a δ-distribution. This is the core consequence of the Sobolev approximation.5
With the Sobolev approximation, the line optical integration simplifies to
 
dz
τ = (κl ρ)rS (152)
d∆ν rS
and only depends on the local conditions at the Sobolev point, i.e. where
 v
ν0 = ν 1 − µ (153)
c rs
With Equation (146), this can be expanded into the final expression for the Sobolev optical depth
!
c/ν0
τ = κl ρ . (154)
(1 − µ2 ) v(r)
r +µ
2 dv(r)
r rS

In supernova ejecta, v(r) = r/t typically holds. In this case, the Sobolev optical depth further
simplifies to  
c
τ = κl ρ . (155)
tν0 rS

8.3 The Elementary Supernova Model

8.3.1 Overview
The elementary supernova model (ES) was developed by Jeffery and Branch and is extensively
described in Jeffery & Branch (1990). It provides a very simple but effective description of line
formation in supernova (SN) ejecta. In its most basic form, the elementary supernova model rests
on the following simplifications and assumptions:

Spherical symmetry This assumption reduces the dimensionality of the problem and according
to polarimetry studies should be globally a reasonable approximation. Naturally, any chemical
inhomogeneities are neglected in this description.
5 Replacing the line profile function with the δ-distribution is mathematically and physically correct when the

plasma state is constant in the interaction region.

20
Homologous expansion Shortly after explosion, the ejecta expand in a force-free ballistic fash-
ion. In this situation, the velocity of a mass element is given by distance to the explosion center,
r, and by the time since explosion

v = r/t. (156)

Photosphere A key simplification of the ES is the division of the SN ejecta into an inner
region, in which continuum processes are important and in which any energy generation processes
occur (radioactivity) and into an outer atmosphere regime in which line formation happens and
continuum processes are neglected. The two regimes are separated by a photosphere, which serves
as a computational boundary at which radiation streams as a black-body field at a certain effective
temperature (often determined from a black-body fit to the observed spectrum) into the ejecta

Iνph = Bν (Teff ). (157)

This is certainly an oversimplification since the opacity in supernova ejecta is typically scattering
dominated. Thus, the thermalization region does not coincide with the photosphere but lies much
deeper inside the ejecta. However, this photosphere description provides still a useful tool for line
identification and for determining relative line strengths.

Pure Scattering Within the ES, line interactions are treated as pure resonance scatterings.
This is equivalent to describing each line transition by a two-level atom. Again, this is an oversim-
plification and only adequate for line transitions in which the involved upper and lower levels do
not strongly interact with other atomic energy levels. This description typically fails for Hα and
also neglects fluorescence and multi-line effects which are crucial for the overall spectrum formation
in SNe Ia.

Simplified line optical depth Finally, the ES assumes that the line optical depth only varies
with the ejecta density and thus neglects any ionization, excitation (and composition) changes
within the ejecta. The line optical depth at any radius is then given by
ρ(r)
τ (r) = τ (rref ) (158)
ρ(rref )
which requires its value at a certain reference location, e.g. at the photosphere. Jeffery & Branch
(1990) discuss in detail the implications and justification of this simplification and conclude that
for ground and metastable levels it is appropriate and for non-metastable levels less so.

8.3.2 Details
The main purpose of the ES is to understand line formation in supernova ejecta. To simply this
discussion, we consider an isolated line transition (the concept is easily generalized to having an
entire series of possible line interactions in the ejecta). Since we are interested in the observable
line profile, we consider only interactions which re-emit photons parallel to the z-axis (we again
work in the impact geometry coordinate system). Now consider an arbitrary location in ejecta with
(z, p) (see illustration in Figure 4). A scattering in a line with rest frame wavelength λ0 occurring
in this region will contribute to the observed (i.e. LF) spectrum at wavelength6
 
1z
λ = λ0 1 − (159)
tc
This is a simple consequence of the requirement that the photon is scattered into the line of sight,
i.e. that
z
µ= p (160)
z + p2
2

and of homology p
r z 2 + p2
v= = . (161)
t t
6 In the ES discussion we will only consider the first order Doppler effect, which is the dominant relativistic

process.

21
 p

scattered into LOS

absorbed

Photosphere z
Rph observer

scattered out of LOS R

Figure 4: Illustration of the Elementary Supernova model geometry.

Equation (159) tells us that contributions to the line profile at a certain (LF) wavelength are
confined to a constant z slice through the ejecta (i.e. parallel to the p-axis). Contributions from the
z > 0 hemisphere will be blue-shifted with respect to the rest-frame wavelength of the transitions.
By contrast, the red-shifted contributions will originate from the z < 0 hemisphere.
For a quantitative prediction of the line profile, we need for every wavelength to integrate over
the respective line forming region Z ∞
Fν = 2π dp pIνem. (162)
0
Here, the emergent specific intensity is given by

 ph
S(r){1 − exp[−τ (r)]} + Iν exp[−τ (r)] if Rph ≤ z ≤ R and p ≤ Rph
em
Iν = S(r){1 − exp[−τ (r)]} if p > Rph and − R ≤ z ≤ R (163)

 ph
Iν else

The optical depth is given by Equation (158) and the source function follows
S(r) = Iνph W (r) (164)
which is simply a consequence of geometric dilution (see Mihalas 1978, p. 120)
 q 
1
W (r) = 1 − 1 − (Rph /r)2 (165)
2
Equations (158) and (162) to (165) build a simple method of calculating line profiles emerging
from SN ejecta. One would typically use a numerical quadrature procedure to solve Equation
(162). After choosing a LF wavelength/frequency for which the line flux should be determined,
Equation (159) is inverted and returns the z-coordinate of the line forming region. If |z| > R, the
integration is trivial, since the resonance region lies outside of the ejecta and only the photosphere
continuum contributes to the emergent spectrum. Otherwise, one would determine the extent of
the resonance region in p space via p
p = R2 − z 2 (166)
and then numerically integrate Equation (162) from 0 to pmax .

22
 Parameter Value
t 3000 s
vmax 0.01c
λ0 1216.7
ve 108 cm s−1
vref 108 cm s−1

Table 2: Parameters used in the ES line profile calculations shown in Figure 5.

8.4 Applications
Single-line illustration Figure 5 shows a few example line profiles determined with this tech-
nique. They only differ in the reference line optical depth, which increases from 0.1 to 1 to finally
10. The remaining parameters are listed in Table 2. Here, we assumed the following density profile
 
vref − v
ρ(v) ∝ exp . (167)
ve

Figure 5: Some example P-Cygni line profiles determined with the ES method. The parameters
are listed in Table 2. Different reference line optical depths were used.

Multi-line application: SYN++ and SYNAPPS The principles underlying single-line for-
mation in the elementary supernova model can be easily generalised to situations with a multitude
of line transitions. In this case, the calculation of the emergent line flux (i.e. Equation 163) is
replaced by a recursive procedure (e.g. Lucy 1999)

Ikr = Ikb exp(τk ) + Sk (1 exp(τk )). (168)

Here, the subscript k marks the contributions of the different line transitions and the superscripts
r and b denote the “red” and the “blue” wing of the line. When calculating the emergent flux at
the LF frequency ν, for each ray, p, one determines the maximum and minimum CMF frequencies
which can scatter into the LOS (using Equation 159). All lines, k, with natural frequency falling
into this interval contribute to the emergent flux and thus determine the extent of the above
recursion process. The procedure then starts at −z end of the ray p and proceeds from one
line resonance point to the next until the ejecta surface facing the observer is reached. Hereby,
Ikb = Ik−1
r
, assuming that the line list is ordered such that natural frequencies decrease with k 7 .
This simple line formation approach underlies the widely used radiative transfer program
Syn++ (and its extension Synapps, see Thomas et al. 2011). The programs are freely avail-
able at [Link] In Figure 6, a simple example calculation is shown.8

Figure 6: Example calculation with Syn++.

Improving spectrum formation in Monte Carlo calculations Lucy (1999) used the exact
same principles just discussed to improve the calculation of emergent spectra from SNe ejecta in
Monte Carlo simulations. The key feature of this approach is to use the Monte Carlo calculation
to determine the line source functions in the ejecta and then determine the final spectrum by
performing a recursive formal integration along the line of sight as outlined above. Following this
approach a virtually noise-free spectrum can be calculated9 , as illustrated in Figure 7.

7 This simple recursion relies on the fact that no line overlaps occur in the Sobolev approximation.
8 Based on the setup specified in the example configuration file shipped with Syn++, syn++.yaml.
9 Note however, that the determination of the source functions on which the integration is based, is still subject

to stochastic fluctuations.

23
Figure 7: Taken from Lucy (1999), Fig. 4. The formal integral method yields the thick solid line
spectrum, which is virtually noise free.

9 Numeric solutions : explicit scattering coupling

If emission and absorption coefficients are numerically known, numeric solution to the transfer
equation in any of its formats is straightforward - not necessarily easy but in principle not prob-
lematic. Such solutions are called formal solutions.
Some very simple problems fall in this category, and the challenge lies mainly in designing as
fast as possible algorithm.
Normally, however, η and χ depend on the radiation field somehow. When this is so, one
solution approach is (lambda) iteration : alternate formal solutions with recomputation of η and
χ. In many physical problems such brute-force split-up of the equation system and solution with
iteration works well.
For some transfer problems this works perfectly well, particularly when optical depths are low
or moderate. Indeed, with modern computing power this approach should probably be the first one
considered. This approach may not be particularly fast. But many alternative ’clever’ approaches
invented in the RT literature since 1960s were sought mainly because computers were too slow
to do this iteration in a reasonable time. This may not at all be the case today, and little is
computationally gained from studying and implementing certain alternative techniques (although
they may still provide physical insight).
However, this is not the whole story. When optical depths are high, lambda iteration can
converge so slowly, that for practical purposes it becomes a non-convergent method. When this is
so, this equation system splitting has to be avoided/reduced. We will first look at a simplified case
of pure scattering in the computational domain.

Figure 8: Illustration of convergence problems in Lambda iteration, from Auer 1991. The true
solition of S is the dotted line. The solid lines show successive iterations. At i = 20, changes are
so slow that an apparent false convergence is obtained. The true solution is eventually obtained :
but only after about 1000 iterations.

One may solve explicit scattering problems by brute force : the equation system is of size
Nd Nµ Nν , and we may attempt matrix inversion, which costs Nd3 Nµ3 Nν3 . But if we have say 100
points in each of these 3 dimensions, we need to store a matrix of size 106 × 106 , which is not
possible. Even if it was, inversion would cost 1018 flops.
Instead, better solution approaches are one of three kinds; the Feautrier method, the Rybicki
method, and the Variable Eddington Factor method. Some confusion then easily arises,
because ’flavours’ of the Feautrier method is applied also as part of the Variable Eddington Factor
method. Different articles and books can appear contradictory because authors differ in their exact

24
definition of terms like ’Feautrier method’.
One should be wary of generic statements in RT theory. Here are nevertheless some comments
on the benefits of differential approaches compared to integral ones from the literature (Mihalas
1986)
• (+) Easy to incorporated new physics
• (+) More naturally posed for dynamic problems and link to radiation hydro codes
• (+) For difference schemes using symmetric and anti-symmetric variabes (such as Feautrier
scheme) : flux can be recovered more robustly as numeric cancellation is avoided. This
improves energy conservation etc.
• (-) Ill-suited for material interfaces and unresolved shock fronts, here generally only integral
approaches survive

9.1 Feautrier method

For illustration of this method, we first consider a problem in planar geometry:
dIν
µ = Iν (µ) − Sν (µ) (169)
dτν
For rays going inward (-) and outward (+):

dIν+
µ = Iν+ (µ) − Sν+ (µ) (170)
dτν
dI −
−µ ν = Iν− (µ) − Sν− (µ) (171)
dτν
Adding:  
dIν+ dI −    
µ − ν = Iν+ (µ) + Iν− (µ) − Sν+ (µ) + Sν− (µ) (172)
dτν dτν
Subtracting:  
dIν+ dI −    
µ + ν = Iν+ (µ) − Iν− (µ) − Sν+ (µ) − Sν− (µ) (173)
dτν dτν
The last equation can also be written
 
  dIν+ dIν−  
Iν+ (µ) − Iν− (µ) =µ + + Sν+ (µ) − Sν− (µ) (174)
dτν dτν

which when put back into Eq. 172, combined with an assumption that S is isotropic (so S + (µ) +
S − (µ) = Sν and S + (µ) − S − (µ) = 0) gives

∂ I¯ν (z, µ)
µ2 = I¯ν (z, µ) − Sν (z) (175)
∂τν2
where
1
I¯ν (µ) ≡ [Iν (µ) + Iν (−µ)] . (0 ≤ µ ≤ 1) (176)
2
Writing out the source function (isotropy not enforced):
R∞R1
¯
2 ∂ Iν (z, µ) σ(z, µ′ − µ, ν ′ − ν)I¯ν ′ (z, µ′ )dµ′ dν ′
µ 2
= I¯ν (z, µ) − Sν (z) + 0 0
t
(0 ≤ µ ≤ 1) (177)
∂τν χν (z, µ)
Note that the demand of scattering isotropy already removed lab-frame differential velocity
field problems, although some tricks salvage some such situations (Rutten, page 118).
At first glance it looks like each equation constains I¯ν ′ (µ′ ) for all angles and frequencies, which
we cannot know unless solving for all angles and frequencies at once, e.g. a full coupling, and full
inversion is needed.
But Feautrier (1964) came up with a method that is significantly cheaper. The essence lies in
that the matrix system is block tridiagonal : for each angle and frequency combination, there is

25
only coupling in space between three points (any second order spatial derivative has this property).
Then, one can come up with a backsubstitution algorithm that is cheaper than brute force matrix
inversion. Thus, Feautrier’s method is in essence a pure ’math trick’. But the method (or very
similar methods) are used also in alternative approaches such as VEF : it is therefore quite broad
in its application and important to understand.
Solution method. For each depth point d, create vectors Īd+1/2 than contain I¯ for each
angle-frequency combination at that depth (so the Īd+1/2 vectors are of length Nµ Nν ). We have,
at depth d, for angle-frequency point number i (Eq 83.17 in M84)
   
1 1 1 1 1 ¯
µ2i I¯d+3/2,i − + I¯d+1/2,i + Id−1/2,i (178)
∆τd+1/2,i ∆τd+1,i ∆τd,i ∆τd+1,i ∆τd,i
P
t
ηd+1/2 + i′ σd+1/2,i,i′ I¯d+1/2,i′
¯
= Id+1/2,i − (179)
χd+1/2,i

Here quantities centred at cell centres have indeces of halv (e.g. 1.5 for the first cell) and
quantities centred at cell edges have integer indeces (e.g. 1 for the first cell).
For each depth, there is spatial derivative coupling between Īd+1/2 and Īd−1/2 and Īd+3/2 for the
same i. Then, we can create (diagonal) matrices at each depth, Ad+1/2 and Cd+1/2 , containing
the discretized derivative factors for I¯d−1/2,i (third term on LHS) and I¯d+3/2,i (first term on LHS),
respectively:
1 1
Aii 2
d+1/2 = µi (180)
∆τd+1/2,i ∆τd,i
ii 1 1
Cd+1/2 = µ2i (181)
∆τd+1/2,i ∆τd+1,i

Create also a (full) matrix Bd+1/2 , which on the diagonal has derivative operators for Īd+1/2
(2nd term on LHS), a -1 (first term on RHS), and a coherent scattering term (for no deflection
coherent scattering i = i′ component of last summation term on RHS). The rest of the matrix is
filled with terms from discretization of the scattering term where angle and/or frequency changes
(so the matrix is full).
 
ii 2 1 1 1 σd+1/2,i,i
Bd+1/2 = −µi + −1+ (182)
∆τd+1/2,i ∆τd,i ∆τd+1,i χd+1/2
ij σd+1/2,i,j
Bd+1/2 = (183)
χd+1/2

On the RHS we obtain a vector Ld+1/2 with the fixed source term (numerically known):
t
ηd+1/2,i
Lid+1/2 = −Sd+1/2,i
t
=− (184)
χd+1/2,i

Note that for the same frequency, each entry in L is the same for isotropy in χ. The transfer
equations for each angle and frequency, for each depth d, is now a matrix system:

− Ad+1/2 Īd−1/2 + Bd+1/2 Īd+1/2 − Cd+1/2 Īd+3/2 = Ld+1/2 (185)

Lower boundary condition. Specify an outgoing intensity Iν+ (τmax , µ) from the lower bound-
ary d = D. One may show (page 370 in Mihalas)

∂ I¯ν
µ = Iν+ (τmax , µ) − I¯ν (τmax , µ) (186)
∂τν taumax
Discretize:
I¯D−1/2,i − I¯D+1/2,i
µi +
= ID+1/2,i − I¯D+1/2,i (187)
∆τD+1/2,i

26
Giving
µi
Aii
D+1/2 = − (188)
∆τD+1/2,i
ii
CD+1/2 = 0 (189)
ii µi
BD+1/2 = −1 (190)
∆τD+1/2,i
ij
BD+1/2 = 0 (191)
+
Lii
D+1/2 = ID+1/2,i (192)

Upper boundary condition. The upper (outer) boundary condition is normally I − = 0. Then
one can show (Eq. 83.20 in M84)
∂ I¯ν (µ)
µ = I¯ν (µ) (193)
∂τν
or (Eq. 83.36 in M84)

I¯5/2,i − I¯3/2,i  −1   

µi = I¯3/2,i 1 + 1/2∆τ3/2.i /µi + ∆τ3/2,i /µi I¯3/2,i − S3/2,i (194)
∆τ2,i

so (to be completed)

Aii
3/2 = 0 (195)
ii µi
C3/2 = (196)
∆τ3/2,i
ii µi
B3/2 = −1 (197)
∆τ2,i
Lii
3/2 = ... (198)

Solution technique. For the first depth point (d = 3/2), we get from Eq. 185

Ī3/2 = B−1 −1
3/2 C3/2 Ī5/2 + B3/2 L3/2 (199)

Define D3/2 = B−1 −1

3/2 C3/2 and X3/2 = B3/2 L3/2 . By substituting Eq. 199 into Eq. 185 for d = 5/2,
we get
Ī5/2 = D5/2 Ī7/2 + X5/2 (200)
where D5/2 can be calculated from the known A5/2 , B5/2 , C5/2 and D3/2 .
The method is to perform a forward-backward sweep to calculate all the D matrices, starting
at the upper boundary d = 1:
1. Compute all the Ad+1/2 , Bd+1/2 , Cd+1/2 matrices at each depth d.

2. Compute D3/2 = B−1

3/2 C3/2

3. Compute X3/2 = B−1

3/2 L3/2
 −1
4. Compute next Dd+1/2 = Bd+1/2 − Ad+1/2 Dd−1/2 Cd+1/2
 −1  
5. Compute next Xd+1/2 = Bd+1/2 − Ad+1/2 Dd−1/2 Ld+1/2 + Ad+1/2 Xd−1/2
6. When you arrive at last depth point D, solve Īd+1/2 = XD+1/2 .
7. Determine all other Īd+1/2 = Dd+1/2 Īd+3/2 + Xd+1/2

Computational cost. Each matrix inversion is now limited to matrixes of size Nµ Nν . We have
4 inversions of them at each depth point, so in total a computational inversion cost of 4Nd Nµ3 Nν3 in-
stead of Nd3 Nµ3 Nν3 in the fully coupled brute force method, a gain by factor Nd2 /4, or 2500 for Nd = 100.
With samplings of order 100 points in each dimension, CPU time (1 GHz) is about 24 hours, but
drops to seconds in the gray or coherent scattering case.

27
9.1.1 Comments.
When is the Feautrier method used? According to Castor 2004, it has been used in the vast
majority of slab geometry work since 1964. A natural application area is when we want to avoid
Lambda iteration (high optical depth) and explicitly include scattering directly in the solution.
But if η and χ still depend on the radiation field, lambda iteration is still needed. Rutten (page)
122) describes that also when S is numerically known (as in Lambda iteration), the Feautrier
method is a very efficient Lambda operator, faster than integral solutions (e.g. using exponential
functions). Thus, the method finds use also in lambda iterations.
The method scales with Nµ3 Nν3 ; thus when very high resolution is needed in these dimensions,
the method may become to expensive. The Rybicki method is an alternative similar method which
is preferred in certain problems where Nν is large (see next section).
Because of the need of isotropy, the method does not work for moving media in the lab frame,
or any other situation with angle-dependent scattering (Castor 2004). However, looking at the
basic equations, it is not clear that this should be a limitation. In the addition equation we can
with no problem retain S(µ) and S(−µ). We will also get first order derivatives in S(µ) and S(−µ),
but this will just add off-diagonal entries in the A and C matrices. Thus - the basic algorithm
still works, but inversion becomes more expensive when A and C are non-diagonal. Discuss this
statement. Feuatrier’s original (2-page) paper uses a known isotropic S, and this may have led to
subsequent confusion that the method is limited to this (?).
One advantage of Feautrier method is that it works with the second order form of the transfer
equation, which is second order accurate and numerically more benign than the first order form
(Castor 2004).

9.1.2 Issue of isotropy

Mihalas 1986 derives an equation (his 3.12) in the observer frame for plane-parallell moving media.
Because of the motions, the source function is not isotropic

∂ 2 I¯∗
µ2 = I¯∗ (µ) − Sν (µ) (201)
∂τν2

Note here a special definition of I¯ is used. Clearly, this can be solved exactly as for when S is
isotropic.
Later on, M86 discusses the more generic case (Eqs 172 and 173), with a reference to Milkey
1975 for discretization.

9.1.3 Spherical symmetry

A similar second-order equation can be derived in spherical symmetry (see more below), so the
method generalizes to spherical symmetry.
This statement appears to hold strictly for isotropic scattering, Mihalas1986 tells us ’Observer-
frame partial redistribution (non-isotropic) calculations in spherical symmetry have never been
carried out using Feautrier variables.’ (but one goes to comoving frame).

9.1.4 2D/3D
A generalization of Feautrier method to 2D (and sketching for 3D) is presented in Cannon 1970.
A statement from Mihalas 1986 is that ’the approach is quite costly in 2 dimensions’. One of the
spatial dimensions can be ’recursed’ away, but matrices of size Nd1 Nµ Nν now need to be inverted,
becoming Nd1 Nd2 Nµ Nν in 3D.

9.1.5 Comoving frame

The Feautrier method (or similar derivatives) can be applied for comoving frame formulations, e.g.
Noerdlinger and Rybicki 1974 (plane parallell), Mihalas and Kunasz 1976 (spherical symmetry).
Noerdlinger and Rybicki 1974 uses the two first order moment equations, rather than combining
them into a second-order one.

28
9.2 Rybicki method
The Rybicki method switches the inner and outer ordering of the Feautrier method, defining
vectors for each angle-frequency point instead of at each depth point. This replaces the Nd Nν3 Nµ3
cost of Feautrier with Nd2 Nν Nµ + Nd3 (see M84 page 377 for details) and is thus preferential for
problems with many angle-frequency points.
Details...
However, it appears there are some subtle caveats with this method for certain problems when
coupling to other physical constraints. (M84 page 380, M86).

9.2.1 2D/3D
A generalization to 2D is presented in Mihalas 1978.

9.3 Variable Eddington Factor method

Assume we have isotropic scattering. In the Feautrier method, this means only that the off-diagonal
elements (with the same frequency) in the B matrix become identical. This does not translate
to any computational speedup, and cost remains at 4Nd Nµ3 Nν3 . Using a different approach where
the moment equations are iterated with formal solutions to determine closure then becomes much
faster. R1 R1
Integration over angle of Eq 177 gives (1/2 0 I¯ν µ2 dµ = 1/2 −1 Iν µ2 dµ = Kν )
R
∂ 2 Kν t R(ν, ν ′ )Jν ′ dν ′
2
= J ν − S ν − (202)
∂τν χν
We have defined earlier the Eddington factor fν = Kν /Jν , so we can write this as
R
∂ 2 f ν Jν t R(ν, ν ′ )Jν ′ dν ′
= J ν − S ν − (203)
∂τν2 χν
This equation has the same form as Eq. 177, except that the variable is now Jν instead of
(angle-dependent) I¯ν . There may still be off-diagonal entries in B due to frequency redistribution.
Thus, Eq. 203 can be solved with a similar backsubstitution scheme as the Feautrier method.
We have got rid of the angles, bringing each iteration down from 4Nd Nµ3 Nν3 to 4Nd Nν3 . However,
iteration is now needed. The total gain is then Nµ3 /Niter , assuming neglegible computation time
for the formal solutions. In many applications a few iterations is sufficient. This benign property
comes from the fact that the Eddington factor measures a degree of asymmetry of the radiation
field; an integral quantity that is robustly recovered in a few iterations.

9.3.1 Spherical symmetry

An equation can be derived by combination of the zeroth and first moment equations (Eq. 83.74
in M84)  
1 ∂ r2 ∂(fν qν Jν ) r2
= (Jν − Sν ) (204)
qν ∂τν qν ∂τν qν
where qν is the sphericality factor, and is fully specified by fν :
Z r
ln qν = [(3fν − 1)/r′ f )ν] dr′ (205)
rc

Eq 204 has again a second order spatial derivative and can be discretized and solved by the
Feautrier scheme.

9.3.2 2D/3D
2D/3D cases? ...

9.3.3 Comoving frame

In the comoving frame in spherical symmetry, two Eddington factors are needed instead of one
(Mihalas1986-II).

29
9.3.4 An example in spherical symmetry : CMFGEN
CMFGEN (Hillier and Dessart 2012) solves time and frequency-dependent RT in spherical sym-
metry (comoving frame), with a velocity gradient, and assumption of isotropic electron scattering.
Bound-bound, bound-free processes are treated implicitly by (Lambda) iterating solutions to en-
ergy and NLTE equations, but electron scattering is treated explicitly in the source function.
It is not clear from Hillier and Dessart 2012 how the zeroth and first moment equations are
solved, given a f factor. The standard method is a a Feautrier-like backsubstitution scheme (see
above), and presumably this is what is used. However, it is also possible that the moment equations
are not actually solved; There code uses a linearization technique, where level populations are solved
by eliminating J using the transfer equations.
The formal solutions to obtain the f factors are used with the tangent ray method.

10 Numeric solutions : Formal solutions

Formal solutions, for numerically known emission and absorption coefficients (or equivalently source
function), need to be carried out either in simple problems where these are truly known, or in iter-
ative methods such as Lambda iteration or Variable Eddington Factor methods. This is generally
quite straightforward and costs of order Nd Nν Nµ in 1D as matrix inversion is generally not needed.

10.1 Planar geometry

...

10.2 Spherical symmetry

Page 380 in M84..

10.2.1 Tangent ray

Cost is ∼ Nd2 Nν in 1D (M86). Nd3 Nν in 2D, Nd4 Nν 3D?
A drawback of the tangent ray method is that it cannot be made fully consistent with the
moment equations (M86); this is possible instead in the discrete space method (although normally
more expensive).

10.2.2 Discrete space

Cost is Nd Nµ3 Nν in 1D (M86), making it preferable when a small number of angle points is needed.
...

10.2.3 An example : CMFGEN

section 3.3 in Hillier2012.

10.3 2D/3D
...

11 Numeric solutions : Lambda iteration

Consider the transport equation
dI
µ =I −S (206)
dτ
For pure scattering S = χJ.
We cannot integrate along rays, because

11.1 Approximate Lambda Iteration

Rutten sec 5.3 good

30
12 Moment methods for neutrino-hydrodynamics
In this chapter we make quite an abrupt jump from radiative transfer to radiation-hydrodynamics
and from photons to neutrinos. In practice, doing radiation-hydrodynamics means employing
significant approximations – often more so than in radiative transfer – in order for the computa-
tional costs to remain tractable. A particularly cumbersome aspect of radiation-hydrodynamics is
the rather complex entanglement between physical assumptions, suitable discretization schemes,
parallelization strategies, and the overall computational efficiency. Likewise to the radiative trans-
fer methods discussed previously, the question how efficient and accurate a given method is can
hardly be answered generally, but instead should be addressed on the basis of a given model
setup and suitably posed diagnostics. In this chapter, we outline the probably most popular
radiation-hydrodynamics schemes, namely so-called moment closure methods, which are based on
the evolution of the angular moments of the radiation intensity. Large parts of this chapter are
adopted from (Just et al. 2015), in which the ALCAR code is described. Several other approximate
methods to incorporate neutrinos in hydrodynamical simulations exist, e.g. neutrino leakage and
light-bulb schemes. In this chapter, however, we will not discuss those.
Please note the following changes of nomenclature in this chapter compared to other chapters:
• We replace the frequency ν by the energy ε = hν
• The symbol ’ν’ is now used to denote (certain species of) neutrinos
• If not specified otherwise, all radiation moments E, F , P , etc. are taken to be in the comoving
frame while the corresponding lab-frame moments are given the subscript “lab”
• Energy-integrated (i.e. gray) quantities are denoted using an overbar, e.g.,
Z
Ē = E dε (207)

• The letter f is now used to denote the flux factor F/(cE).

12.1 Introduction: Neutrino transport in core-collapse supernovae and

neutron-star mergers
Compared to the typical cross section of photons, σT ≈ 6.65 × 10−25 cm2 , the cross section of
neutrinos is extremely tiny, σw ≈ (ε/(me c2 ))2 × 10−44 cm2 . The associated mean free path in an
environment with n nucleons per unit volume is roughly
 
1 ρ
λν = ∼ 1020 cm . (208)
nσw 1 g cm−3
Hence, under “normal” conditions in a stellar environment neutrinos interact only once, namely
when they are produced. From the technical point of view, radiation-hydrodynamics in this
optically-thin regime is most straightforward because the fluid is only subject to local losses of
energy and lepton number, the rates of which can be computed solely from the thermodynamic
properties of the fluid.
The situation changes dramatically when matter is in more extreme conditions, such as in
core-collapse supernovae (CCSNe) or neutron-star (NS) mergers (see Figs. 9, 10). The following
characteristic properties can be identified in both environments:
• Once densities of ρ ∼ 1011−12 g cm−3 are reached, the neutrino mean free path becomes
comparable to characteristic length scales of the radiating object (such as the pressure scale
height) and the neutrino transport reaches diffusion conditions (small diffusive flux, diagonal
Eddington tensor).
• The surface (radius) at which neutrinos effectively decouple from the medium is called the
neutrino surface (neutrinosphere radius). Here, the strongest cooling of matter by neutrinos
takes place.
• The net-cooling region is enclosed by a net-heating region, at which a small fraction of
neutrinos (typically a few percent) are reabsorbed. The surface (radius) between both regions
is called the gain surface (gain radius).

31
Figure 9: Six phases of neutrino production and its dynamical consequences in a core-collapse
supernova (from top left to bottom right). In the lower halves of the plots the composition of the
stellar medium and the neutrino effects are sketched, while in the upper halves the flow of the
stellar matter is shown by arrows. Inward pointing arrows denote contraction or collapse, outward
pointing arrows expansion or mass ejection. Radial distances R are indicated on the vertical axes,
the corresponding enclosed masses M (r) are given on the horizontal axes. RFe , Rs , Rν , Rg , and
Rns denote the iron-core radius, shock radius, neutrinospheric radius, gain radius (which separates
neutrino cooling and heating layers), and proto-neutron star (PNS) radius, respectively. MCh
defines the effective Chandrasekhar mass, Mhc the mass of the homologously collapsing inner core
(where velocity u ∝ r), ρc the central density, and ρ0 ≈ 2.7 × 1014 g cm−3 the nuclear saturation
density. (Figure taken from Janka et al. 2007) (Figure caption taken from Janka 2017)

Figure 10: Sketch of the physical conditions in a black-hole torus systems as a possible remnant of
a neutron-star neutron-star or a neutron-star black-hole merger.

• Scientifically, neutrino heating is of paramount importance: In the “neutrino-driven mecha-

nism” of CCSNe it is responsible for the revival of shock expansion and therefore the entire
explosion. Moreover, after explosion, or during the evolution of the remnant of a NS merger
neutrino heating leads to the acceleration of a thermal wind, called neutrino-driven wind,
which is hypothesized as an important site of heavy element nucleosynthesis.
• The main contribution of neutrino heating takes place right above the gain surface in the
semi-transparent regime, in which neutrinos can neither be described well by diffusion nor by
free streaming. Thus, an accurate transport scheme is highly desirable to obtain conclusive
results concerning the impact of neutrino heating.

12.2 Specifics of neutrino transport

We summarize some distinct properties and technicalities encountered when dealing with neutrino
transport:

• Neutrinos come in six different species νe , ν̄e , νµ , ν̄µ , ντ , ν̄τ . In principle, one needs to solve
transport equations for each species. In practice, however, it often suffices to evolve the four
non-electron type neutrinos using a single representative species (“X-neutrinos”), because
absorption/emission processes with muons and tauons are subdominant compared to those
with electrons.
• As opposed to photons, each lepton family (electron-, muon-, tau-family) satisfies a number
conservation law, which should be fulfilled also by the numerical transport scheme. In par-
ticular, the sum of the numbers of electrons plus electron neutrinos minus positrons minus
electron antineutrinos is conserved.
• Neutrinos satisfy Fermi statistics. Hence, the equilibrium distribution differs from that of
photons and Fermi blocking effects need to be taken into account. Also, (electron) neutrinos
can be ascribed a non-vanishing chemical potential such that in chemical equilibrium (also
called β-equilibrium)

µνe + µn = µp + µe− (209)

holds.
• In photon transport one often deals with line opacities (in stellar atmospheres, stellar ejecta)
or energy-independent Thomson scattering (accretion disks). In neutrino transport, in con-
trast, energy-dependent continuum opacities dominate, most importantly that of electron-
/positron-capture and their inverse, which are ∝ ε2 . See Table 3 for a summary of the
most important interaction channels. Due to the strong energy-dependence, gray neutrino-
transport schemes typically underperform in capturing the mean energies and therefore the
net-neutrino heating rates.

32
Table 3: Most important neutrino processes in supernova and proto-neutron star matter. (Table
taken from (Janka 2017))
Process Reactiona
Beta-processes (direct URCA processes)
electron and νe absorption by nuclei e− + (A, Z) ←→ (A, Z − 1) + νe
electron and νe captures by nucleons e− + p ←→ n + νe
positron and ν̄e captures by nucleons e+ + n ←→ p + ν̄e
“Thermal” pair production and annihilation pro-
cesses
Nucleon-nucleon bremsstrahlung N + N ←→ N + N + ν + ν̄
Electron-position pair process e− + e+ ←→ ν + ν̄
Plasmon pair-neutrino process γ
e ←→ ν + ν̄
Reactions between neutrinos
Neutrino-pair annihilation νe + ν̄e ←→ νx + ν̄x
Neutrino scattering νx + {νe , ν̄e } ←→ νx + {νe , ν̄e }
Scattering processes with medium particles
Neutrino scattering with nuclei ν + (A, Z) ←→ ν + (A, Z)
Neutrino scattering with nucleons ν + N ←→ ν + N
Neutrino scattering with electrons and positrons ν + e± ←→ ν + e±
a
N means nucleons, i.e., either n or p, ν ∈ {νe , ν̄e , νµ , ν̄µ , ντ , ν̄τ }, νx ∈ {νµ , ν̄µ , ντ , ν̄τ }

12.3 Governing equations

The equations for the radiation transport part depend on the choice of the frame in which certain
quantities are defined in (see Sec. 5). The comoving-frame moment equations of order O(v/c) can
be derived from the comoving-frame transfer equation (i.e. Eq. 56 generalized to the multi-D case)
by performing the same angular integration operations as done when computing E and F from I:

∂t E + ∇j F j + ∇j (v j E) + (∇j vk )P jk − (∇j vk )∂ε (εP jk ) = C (0) , (210a)

i 2 ij j i j i ijk (1),i
∂t F + c ∇j P + ∇j (v F ) + F ∇j v − (∇j vk )∂ε (εQ )=C . (210b)

Note that t and x (with respect to which the spatial derivatives are taken) are Eulerian (lab-frame)
coordinates, while ε, E, F i as well as all quantities entering the collision integrals C (0) , C (1),i are
measured in the comoving frame. Note also that in deriving Equations (210b) time derivatives of
the velocity (i.e. acceleration terms) as well as products of the velocity and time derivatives of
radiation moments are ignored, because they are considered to be O(v 2 /c2 ) based on dimensional
arguments (e.g. Mihalas & Mihalas 1984; Rampp & Janka 2002).
The above equations couple to the hydrodynamics equations,

∂t ρ + ∇j (ρv j ) = 0 , (211a)
j
∂t (ρYe ) + ∇j (ρYe v ) = QN , (211b)
i i j
∂t (ρv ) + ∇j (ρv v + Pg ) = , QiM (211c)


∂t et + ∇j v (et + Pg ) = QE + vj QjM .
j
(211d)

(with mass density, ρ, electron-fraction Ye , gas pressure Pg , and total fluid energy et ≡ ρv 2 /2 + ei ,
and fluid internal energy ei ) by means of the source terms
X
QE = − C̄ (0) , (212a)
species
1 X (1),i
QiM = − C̄ , (212b)
c2
species
Z ∞ " (0)   (0)  #
C C
QN = −mB − dε , (212c)
0 ε νe ε ν̄e

(with mB being the atomic mass unit) describing exchange of energy, momentum and lepton
number, respectively, between the fluid and neutrinos.

33
 Another widely used formulation of the above moment equations is the “mixed-frame” formula-
tion (e.g. Hubeny & Burrows 2007; Krumholz et al. 2007; Cardall et al. 2013), in which the radi-
ation moments are defined in the lab-frame, while all interaction coefficients entering the collision
integral are in the comoving frame. It is a priori not clear which formulation is more convenient.
However, in the mixed-frame formalism the treatment of dynamic diffusion (advection of radiation
under optically thick conditions) is somewhat more difficult because the lab-frame flux

i 4
F̄lab = F̄ i + v i Ē + vj P̄ ij → F̄ i + v i Ē , (213)
3
is dominated by the advection flux ∝ v i Ē, which may lead to large truncation errors for the
subdominant comoving-frame radiation flux F̄ i . On the other hand, the mixed-frame equations
provide the numerically convenient property that the radiation moments truly decouple from the
fluid in the optically thin limit (of vanishing radiation–matter interaction), while when using the
comoving-frame equations the advection/aberration terms couple the moments to the local fluid
velocity in all regimes.

12.4 Moment closure methods

The full information contained in the Boltzmann equation can be captured equally well by an
infinite series of conservation equations for the angular moments, in which the evolution equation
for a moment of rank m contains the moment of rank m+1 within the divergence operator. Instead
of solving the infinite series of moment equations, the series can be truncated at the level of the
(m + 1)-th moment, provided the (m + 1)-th moment is available to close the set of m equations.
In order to determine the (m + 1)-th moment in a manner fully consistent with the Boltzmann
equation, one needs to solve the latter additionally. However, since its solution is only used for
the closure, in practice it is sufficient to solve a somewhat simpler, downgraded version of the
Boltzmann equation, also called “model-Boltzmann equation”. This approach is pursued in the
VERTEX code (Rampp & Janka 2002), which employs a tangent-ray Boltzmann-solver to compute
the closure. The accuracy of this algorithm is comparable to one solving directly the full Boltzmann
equation; it is hence called a “Boltzmann-solver”. A limitation of VERTEX is, however, that it is
not fully multidimensional, but instead it solves several 1D transfer problems quasi-independently
along radial rays (called “Ray-by-Ray-plus” approach).
A computationally cheaper and fully multidimensional, though much more approximate option
is to assume that an algebraic closure relation holds between the evolved moments and the (m+ 1)-
th moment. This is what defines the algebraic closure methods, such as flux-limited diffusion
(FLD) and the algebraic-Eddington-factor method (AEF or often called M1). Essentially, this
corresponds to imposing additional conditions or symmetries on the local radiation field. The
consequence is that two (out of seven in the general case) independent variables describing the
angular dependence of the radiation field disappear from the treatment. Evidently, the tradeoff for
this computational simplification is that an algebraic closure method may strongly vary in quality
between different physical setups. For example, since in an algebraic closure method the consistent
evolution of higher-order angular moments is ignored, it appears likely that the quality of the
scheme appreciably depends on the geometric complexity of the radiation field, or equivalently on
the shape and number of the individual radiation sources. Moreover, connected to this issue is the
circumstance that in the optically thin limit of vanishing source terms an algebraic moment scheme
is generally not able to accurately describe the unperturbed superposition of multiple radiation
fronts, simply on account of the closure being a local, non-linear function of the evolved quantities.
Despite these conceptual deficiencies algebraic closure methods can in many cases offer an excellent
compromise between efficiency and accuracy when performing energy-dependent, multidimensional
radiation-hydrodynamics simulations.
Independent of the rank at which the scheme is truncated, any closure prescription should agree
with certain consistency requirements that directly follow from the definition of the moments or
from the Boltzmann equation. Using the normalized moments f ≡ F/(cE), where f ≡ |f| is the
flux-factor, Dij ≡ P ij /E, the Eddington tensor, and q ijk ≡ Qijk /(cE), the normalized 3rd-moment

34
tensor, it follows from the definition of the angular moments that
|f| ≤ 1 , (214)
ij
D ≤ 1, (215)
X
Djj = 1 , (216)
j

|q ijk | ≤ 1 , (217)
X X X
q ijj = q jij = q jji = f i , (218)
j j j

must hold at any time. In the free-streaming limit all of the radiation propagates into a single
direction away from its source and it must hold
f = 1 , Dij = niF njF , , q ijk = niF njF nkF , (219)
with niF ≡ F i /|F| denoting the direction of the flux density. In the opposite limit of very frequent
interactions, i.e. in the diffusion limit, the specific intensity is approximately isotropic. Ignoring
velocity terms, the radiation moment equations degenerate in this limit to the diffusion equation
 
c
∂t E + ∇i − ∇ E = C (0)
i
(220)
3κtra
and the relations
1 ∇E 1 ij
f = − , Dij = δ . (221)
3κtra E 3

12.4.1 Flux-limited diffusion method

The approach of FLD is to truncate the set of moment equations at the level of the 1st-moment
equation and to derive an expression for the flux density based on the diffusion limit described by
Eqs. (220) and (221). Introducing the “Knudsen number” R = |R|, with
1 ∇E
R≡ , (222)
ω κtra E
where κtra is the transport opacity and
ω ≡ (κs E + κa E eq )/(κtra E) (223)
is called “effective albedo”, the flux density FFLD is prescribed as
FFLD = − Λ(R) RcE , (224)
in which Λ(R) is called the “flux-limiter”. The latter is constructed such that the flux density
correctly preserves the constraints of Eqs. (219), (221). To this end the limits
lim Λ(R) R = 1, (225)
R→∞

and
1
lim Λ(R) = , (226)
R→0 3
respectively, have to be fulfilled (note that ω → 1 in the diffusion limit). Three prescriptions are
widely used in the literature (Wilson et al. 1975; Liebendörfer et al. 2004; Levermore & Pomraning
1981):
1
ΛWilson (R) = , (227a)
3+R
ΛBruenn (R) =
  √ 
min Λ 1+ 1−(rν /r)2
Wilson (R), , r > rν
2R
(227b)

ΛWilson (R) , else ,
 
1 1
ΛLevermore (R) = coth R − . (227c)
R R

35
The limiter in Eq. (227b) is only designed for the spherically symmetric case in which r is the
radius coordinate and rν is the radius of a (properly defined) neutrinosphere. The limitation for
r > rν is intended to ensure that the neutrino flux cannot be higher than if the neutrinos were
distributed isotropically into a finite cone subtending the sphere of radius rν .
In order to evolve Eq. (210a), one also needs the 2nd-moment tensor P ij . This can be obtained
from additional assumptions relating the flux-limiter with the Eddington factor (Levermore 1984).
Essentially, one uses the energy- and flux-density to compute the Eddington factor using the closure
relations given in Sec. 12.4.2.
The main drawbacks of FLD are: First, the prescription of the flux density is in general not
consistent with the 1st-moment equation. As a direct consequence, the full RHD system suffers
from momentum and therefore energy non-conservation whenever momentum transfer between
matter and radiation takes place. Second, in more than one dimension a complication arises from
the fact that the flux density vector is always directed opposite to the gradient of the energy density
since the pressure is effectively isotropic (cf. Eq. (221)): Radiation in the free-streaming limit will
not keep its original flux direction after closely passing opaque objects. Instead it behaves like a
gas and fills up space in every direction, unable to form persistent shadows.
The third issue is a purely computational aspect: The energy equation evolved in FLD is – at
least whenever f 6= 1 – of parabolic mathematical nature. As such, it comes with the property that
the operator ∇ · F needs to be treated time-implicitly in most practical cases. This is because the
characteristic timescale τFLD associated with ∇·F can become extremely short in the optically thin
limit κtra → 0. One can roughly estimate the local timescale τFLD by thinking of the operator ∇·F
as being locally a linear convex combination of the advection operator αfc∇E and the diffusion
operator (1 − α)(−Λc/κtra )∇2 E, with some weighting factor 0 < α < 1. A heuristic dimensional
analysis then gives
E c c
∼ α f E + (1 − α) ΛE (228a)
τFLD ∆x κtra ∆x2
 −1
αf (1 − α)3Λ
⇒ τFLD ∼ + , (228b)
τadv τdiff
where ∆x is the local grid size and
τadv ≡ ∆x/c , (229a)
2
τdiff ≡ 3κtra ∆x /c (229b)
are the characteristic timescales of advection and diffusion, respectively. Hence, owing to the fact
that τdiff → 0 in optically thin regions, τFLD can drop significantly below τadv .

12.4.2 Algebraic Eddington factor method

In the AEF method the truncation of the moment equations takes place at the level of the 2nd-
moment equation, i.e. the Eddington tensor is expressed as a function of evolved quantities. Besides
resolving by construction the first two drawbacks of the FLD scheme mentioned in Sec. 12.4.1,
another important computational difference to the FLD scheme is that the equations solved in the
AEF scheme are intrinsically hyperbolic, which allows to use explicit time integration methods (at
least for all but the source terms) with time steps not lower than the minimum of τadv over the
computational domain.
Below we list some closure relations expressing the Eddington factors as functions of the flux
factors, f , and normalized energy densities, e = (hc)3 /(4πε)3 E (Minerbo 1978; Dubroca & Feugeas
1999; Janka 1991; Cernohorsky & Bludman 1994):
1 1

χMinerbo(f ) = + 6f 2 − 2f 3 + 6f 4 , (230a)
3 15
3 + 4f 2
χM1 (f ) = p , (230b)
5 + 2 4 − 3f 2
 
1 1 3
χJanka (f ) = 1 + f 1.31 + f 3.56 , (230c)
3 2 2
 
2 f 1
χ[Link]. (e, f ) = (1 − e)(1 − 2e)σ + , (230d)
3 1−e 3

36
where σ(x) ≡ x2 (3 − x + 3x2 )/5 in Eq. (230d).
To extend the one-dimensional Eddington factor χ to the multidimensional Eddington tensor
Dij and the 3rd-moment tensor q ijk , we make use of the underlying assumption that these tensors
are only a function of the scalar E and the vector F. From symmetry arguments it follows that
Dij must be symmetric with respect to rotation around the flux direction, nF ≡ F/|F|, which
is only fulfilled if Dij is a linear combination of the two tensors δ ij and niF njF . After using the
trace condition of Eq. (216) the two remaining coefficients can be expressed as functions of a single
parameter, χ, which is the multidimensional generalization of the Eddington factor and is defined
as R 2
dΩ (n · nF ) F
χ≡ R . (231)
dΩ F
The Eddington tensor then reads (e.g. Levermore 1984):
1 − χ ij 3χ − 1 i j
Dij = δ + nF nF . (232)
2 2
The energy-dependent comoving-frame moment equations, Eqs. (210), also contain the 3rd
moments, q ijk . In analogy to the derivation of Dij above, the condition that this tensor only
depends on E and F must result in q ijk being invariant under rotation around nF , which is
only fulfilled if q ijk is a linear combination of niF δ jk , njF δ ik and nkF δ ij as well as niF njF nkF . The
corresponding coefficients can be eliminated using the trace conditions of Eq. (218) in favor of a
single parameter, q, defined as R 3
dΩ (n · nF ) F
q≡ R , (233)
dΩ F
yielding for the 3rd-moment tensor:
f − q i jk 5q − 3f i j k
q ijk = (nF δ + njF δ ik + nkF δ ij ) + nF nF nF . (234)
2 2
The 3rd-moment factor q explicitly depends on the distribution function. Therefore, only closures
that dictate an explicit functional form of the distribution function are suited for the computation
of the 3rd moment, unless additional assumptions are made in the construction of the closure. For
the Minerbo closure, the factor q can be calculated in a straightforward manner in analogy to the
derivation of χ (Minerbo 1978) and reads
f
qMinerbo (f ) = (45 + 10f − 12f 2 − 12f 3 + 38f 4 − 12f 5 + 18f 6 ) . (235)
75

12.4.3 Numerical treatment of the AEF/M1 scheme

Owing to the fact that the advection-type operators on the left-hand side of the two-moment sys-
tem, Eqs. (210), are of hyperbolic mathematical nature, we can employ Godunov-type finite-volume
methods commonly used in numerical hydrodynamics to discretize these operators. However, in
regions of strong coupling with matter the source terms become stiff and the moment equations
approach the parabolic diffusion limit. Hence, the time integration is performed in a mixed explicit–
implicit manner, in which all terms on the left-hand side of Eqs. (210) are treated explicitly in
time while the source terms on the right-hand side of Eqs. (210) are handled implicitly (at least
whenever being in the stiff regime). In that way the overall time step used for integration of the
whole scheme is constrained by the Courant condition to be on the order of the advection timescale
τadv ≡ ∆x/c, i.e. the light-crossing time of grid cells with width ∆x. The alternative would be
to integrate the full two-moment system implicitly in time. In that case the computational cost
per time step would be significantly higher (particularly in the multidimensional case) but on the
other hand this would allow to employ a larger time step which is closer to the fluid-dynamical
time step ∆x/v.

Hyperbolic part The notion is to exploit a Godunov method as the basis for a high-resolution
shock capturing scheme that solves the local Riemann problems between discontinuous states at
the interfaces of grid cells. We start the presentation of its working method by considering the
one-dimensional system    
E F
∂t + ∂x = 0, (236)
F c2 χE

37
where the algebraic closure χ = χ(e, f ) is a function of e and f . This system is hyperbolic if the
Jacobian matrix J of the vector of fluxes (F, c2 χE)T ,
 
0 1
J = , (237)
c2 (χ + e ∂χ ∂χ ∂χ
∂e − f ∂f ) c ∂f

has real eigenvalues λ1D

± , given by
s
c ∂χ c ∂χ 2 ∂χ ∂χ
λ1D
± = ± + 4(χ + e −f ). (238)
2 ∂f 2 ∂f ∂e ∂f

All of the closures listed in Eqs. (230) fulfill the condition of hyperbolicity and lead to the following
properties: In the free-streaming limit, f → 1, we have
∂χ
χ = 1 , λ1D 1D
+ = +c , λ− = ( − 1)c , (239)
∂f
while in the diffusion limit, f → 0, one obtains
1 1
χ= , λ1D
± = ±√ c . (240)
3 3
That is, the limiting cases for the Eddington factor and the wave speeds are consistent with what
is dictated by the Boltzmann equation.
In the multidimensional generalization of Eq. (236) the matrix eigenvalues contain an additional
dependence on the direction cosine µ ≡ cos αF , where αF is the angle between the direction of the
radiation flux vector F and the coordinate direction with respect to which the derivative is taken.
The wave speeds are now obtained as roots of a cubic polynomial leading, at least in terms of a
general closure, to rather large expressions. For practical purposes we do not take into account
the exact angular dependence of the eigenvalues λexact
± (µ) but we instead approximate the latter
using the following 1st-order expansion in µ:


λ± (µ) = λexact
± (0) + |µ| λexact
± (1) − λexact
± (0) , (241)

where

λexact
± (1) = λ1D
± , (242a)
s
exact c ∂χ 1 ∂χ
λ± (0) = ± 2(1 − χ − e ) + (1 + 2f 2 − 3χ) . (242b)
2 ∂e f ∂f

In a fashion that is commonly employed in numerical hydrodynamics, we use the above velocities
as signal speeds for an approximate Riemann solver in order to compute the numerical fluxes
through each cell interface. We use the two-wave solver by Harten, Lax and van Leer (“HLL-
solver”), which approximates the final numerical interface fluxes as functions of the left-/right-hand
side fluxes FL/R (with F ∈ {F i , c2 P ij }) and states UL/R (with U ∈ {E, F i }) as


HLL λHLL L
+ F − λ− F
HLL R λHLL
+ λ−
HLL
UR − UL
F ≡ + , (243)
λHLL
+ − λHLL
− λHLL
+ − λHLL
−

with the signal speeds λHLL

+ + ) and λ−
= max(0, λL+ , λR HLL
− ). All quantities labeled by
= min(0, λL− , λR
L/R in this flux formula are computed using the cell-interface reconstructed moments Ê L/R , F̂ i,L/R .
Yet, there is a caveat we have to face when approaching the parabolic diffusion limit (cf.
Eq. (221)) with the scheme described above, since the latter is originally designed only for hyper-
bolic systems. In contrast to the hyperbolic system, the parabolic diffusion equation is not asso-
ciated with characteristic waves propagating information between cells with finite speeds. Hence,
the ansatz of using a Riemann solver that tracks characteristics via upwinding and captures shocks
by adding diffusivity is no longer justified in the parabolic diffusion regime. Instead, the fluxes in
the diffusion regime should be of central type (i.e. symmetric with respect to the cell interface)
and they should lead to as little as possible numerical diffusivity in order not to spoil the effects

38
Figure 11: Neutrino radiation field of a static proto-neutron star. We compare for two models
(see top title) between different transport schemes (see line labels in the second panel from the
top). In the panels are displayed from top to bottom the properties (density ρ, temperature T
and electron fraction Ye ) of the hydrodynamic q background, the mean flux-factor f¯ ≡ F̄ /(cĒ),
R R
luminosity L ≡ 4πr2 F̄ and rms-energy hǫirms ≡ εE(ε)dε/ N (ε)dε of electron neutrinos, and
the source terms QE /nB , QN /ρ for the gas energy density (cf. Eq. (212a)) and electron-number
density (cf. Eq. (212c)), respectively. Note that in cases when the dotted green line is invisible it
is superimposed by the solid green line.

of the physical diffusivity. To handle this issue, we employ a simple switch between the two types
of fluxes according to: 
FHLL 1 if Pi+ 12 < 1 ,
i+
2 
FHLL,∗
i+ 12
= 1
(244)
 F 1 +F 1
2
L
i+
R
i+
if Pi+ 21 > 1 ,
2 2

where the index “i” denotes a representative grid index for any coordinate direction and the “stiffness
parameter”
∆x
P ≡ κtra ∆x = , (245)
λν
is a measure of the degree of neutrino–matter coupling relative to numerically resolved scales of
length and time: For P > ∼ 1 neutrino interactions proceed on spatial and temporal scales smaller
than the grid scale ∆x and shorter than the numerical time step ∆x/c, respectively. Hence, for
P exceeding unity the source terms become stiff and thereby cause the two-moment system to
undergo the transition from a hyperbolic to a parabolic system.

Source terms The numerical integration of the interaction source terms deserves special care
because the characteristic neutrino-interaction timescale,

τint ≡ (cκtra )−1 = λν /c , (246)

can become shorter than τadv by up to many orders of magnitude. In this case, i.e. for stiffness
parameters P > 1 (cf. Eq. (245)), the moment equations are stiff and a fully explicit time inte-
gration would lead to numerical instability. Hence, the source terms make an implicit treatment
indispensable.

12.4.4 Selected test problems

Neutrino radiation field of a static proto-neutron star
See Fig. 11.

Neutrino field around a torus

See Fig. 12.

References
Cardall, C. Y., Endeve, E., & Mezzacappa, A. 2013, Phys. Rev. D, 87, 103004
Castor, J. I. 1970, MNRAS, 149, 111
Cernohorsky, J. & Bludman, S. A. 1994, ApJ, 433, 250
Dubroca, B. & Feugeas, J. 1999, Academie des Sciences Paris Comptes Rendus Serie Sciences
Mathematiques, 329, 915
Hubeny, I. & Burrows, A. 2007, ApJ, 659, 1458
Hummer, D. G. & Rybicki, G. B. 1985, ApJ, 293, 258
Janka, H.-T. 1991, PhD thesis, , Technische Universität München, (1991)

39
Figure 12: Comparison of the neutrino radiation field emerging from the configuration of model
M3A8m3a5 at t = 50 ms as computed with the ray-tracing scheme (labeled RT) and with the
algebraic-Eddington-factor method (labeled AEF). Panels (a), (d) show color coded the energy-
integrated energy densities, E, and Panels (b), (e) the (absolute values of) flux densities, F , of
electron neutrinos for both schemes, while in Panels (c), (f) the respective relative differences
of the latter quantities between both schemes are depicted. The arrows in Panels (b) and (e)
indicate the flux-density vectors multiplied by 4πr2 scaled such that the maximum arrow length
corresponds to ∼ 4 × 1053 erg s−1. The mean energies, ε, listed in Panels (a), (b) are computed
as ratios of luminosity to total number flux, both calculated as integrals of the corresponding
flux densities over a sphere at radius r = 300 km. The luminosities, L, are given on top of
Panels (b) and (e). Panels (g), (i) show the net heating and cooling rates due to β-processes with
both electron-neutrino species, Qν . Finally, Panels (h), (j) display the electron fractions for local
neutrino-capture equilibrium, Yeν . The yellow line in each panel demarcates the net cooling from
the net heating region, i.e. it coincides with Qν = 0.

Janka, H.-T. 2017, ArXiv e-prints

Janka, H.-T., Langanke, K., Marek, A., Martinez-Pinedo, G., & Mueller, B. 2007, Phys. Rept.,
442, 38
Jeffery, D. J. & Branch, D. 1990, in Supernovae, Jerusalem Winter School for Theoretical Physics,
ed. J. C. Wheeler, T. Piran, & S. Weinberg, 149
Just, O., Obergaulinger, M., & Janka, H.-T. 2015, MNRAS, 453, 3386
Krumholz, M. R., Klein, R. I., McKee, C. F., & Bolstad, J. 2007, ApJ, 667, 626
Lamers, H. J. G. L. M. & Cassinelli, J. P. 1999, Introduction to Stellar Winds (Cambridge, UK:
Cambridge University Press)
Levermore, C. D. 1984, J. Quant. Spec. Radiat. Transf., 31, 149
Levermore, C. D. & Pomraning, G. C. 1981, ApJ, 248, 321
Liebendörfer, M., Messer, O. E. B., Mezzacappa, A., et al. 2004, ApJS, 150, 263
Lucy, L. B. 1999, A&A, 345, 211
Mihalas, D. 1978, Stellar Atmospheres, 2nd edn. (San Francisco: W. H. Freeman and Co)
Mihalas, D. & Mihalas, B. W. 1984, Foundations of Radiation Hydrodynamics (New York: Oxford
University Press)
Minerbo, G. N. 1978, J. Quant. Spec. Radiat. Transf., 20, 541
Rampp, M. & Janka, H. 2002, A&A, 396, 361

Rybicki, G. B. & Hummer, D. G. 1978, ApJ, 219, 654

Rybicki, G. B. & Lightman, A. P. 1979, Radiative Processes in Astrophysics (New York: Wiley-
Interscience)
Sobolev, V. V. 1960, Moving Envelopes of Stars (Cambridge, MA: Harvard University Press)
Thomas, L. H. 1930, The Quarterly Journal of Mathematics, os-1, 239
Thomas, R. C., Nugent, P. E., & Meza, J. C. 2011, PASP, 123, 237
Wilson, J. R., Couch, R., Cochran, S., Le Blanc, J., & Barkat, Z. 1975, in Annals of the New York
Academy of Sciences, Vol. 262, Seventh Texas Symposium on Relativistic Astrophysics, ed. P. G.
Bergman, E. J. Fenyves, & L. Motz, 54–64

40