Radiative Transfer for the FIRST ERA, Astrofizyka i kosmologia
[ Pobierz całość w formacie PDF ]
261
RADIATIVE TRANSFER FOR THE FIRST ERA
J. Trujillo Bueno
1
,
2
1
Instituto de Astrofısica de Canarias, E-38200, La Laguna, Tenerife, Spain
2
Consejo Superior de Investigaciones Cientıficas, Spain.
Abstract
This paperpresentsabriefoverviewofsomerecentad-
vancesin numericalradiativetransfer,whichmayhelpthe
molecular astrophysics community to achieve new break-
throughs in the interpretation of spectro-(polarimetric)
observations.
to the transfer of polarized radiation in magnetized plas-
mas including anisotropic pumping (Trujillo Bueno, 1999;
2001), which may be of interest for modelling polarization
phenomenainMASERS.ThecaseofRTinmolecularlines
is presented in an extra contribution at this conference by
Asensio Ramos, Trujillo Bueno and Cernicharo (2001).
Keywords:Methods:numerical–radiativetransfer–Stars:
atmospheres – Missions: FIRST
2. RT methods based on Gauss-Seidel iteration
The essential ideas behind the iterative schemes on which
our NLTE multilevel transfer codes are based on can be
easily understood by considering the “simplest” NLTE
problem: the coherent scattering case with a source func-
tion given by
1. Introduction
The development of novel Radiative Transfer (RT) meth-
ods often leads to important breakthroughs in astrophysi-
cal plasma spectroscopy because they allow the investiga-
tion of problems that could not be properly tackled using
the methods previously available. This RT topic is also of
great interest for the “Promise of FIRST”, since a rigor-
ous interpretation of the observations will require to carry
out detailed confrontations with the results from NLTE
RT simulations in one-, two-, and three-dimensional ge-
ometries.
The e/cient solution of NLTE multilevel RT problems
requires the combination of a highly convergent iterative
scheme with a very fast formal solver of the RT equation.
This applies to the caseof unpolarizedradiationin atomic
lines, to the promisingtopic ofthe generationand transfer
of polarized radiation in magnetized plasmas and to RT
in molecular lines.
The “dream” of numerical RT is to develop iterative
methodswhereeverythinggoesassimplyaswith thewell-
known Λ-iteration scheme, but for which the convergence
rate is extremely high. In this contribution we present
an overview of some iterative methods and formal solvers
we have developed for RT applications. Our RT methods
are based on Gauss-Seidel iteration and on the
non-linear
multigrid method. These new RT developments are of in-
terest because they allow the solution of a given RT prob-
lem with an order-of-magnitudeof improvementin the to-
tal computational work with respect to the popular ALI
method (on which most present NLTE codes are based
on). Our RT methods have been succesfully applied to as-
trophysical problems of unpolarized radiation in atomic
lines (in 1D, 2D and 3D with multilevel atoms) and also
S=(1
−
)J +
B
,
(1)
with
the NLTE parameter, J the mean intensity and
B the Planck function. The mean intensity at point “
i
”
is the angular average of incoming (“
in
”) and outgoing
(“
out
”) contributions. For instance, for a one-point angu-
lar quadrature
J
i
=J
i
in
+J
i
out
≈
1
2
(I
i
i
+I
out
i
)
.
(2)
iteration scheme is to do the fol-
lowing in orderto obtain the “new” estimate of the source
function at each spatial grid-point “
i
”:
S
new
i
−
=(1
−
)J
old
i
+
B
i
,
(3)
i
is the mean intensity at the grid-point “
i
”cal-
culated using the previous known values of the source
function (i.e. using S
ol
i
).
For a given spatial grid of NP points the formal so-
lution of the transfer equation can be symbolically repre-
sented as
I
Ω
=
Λ
Ω
[
S
]+
T
Ω
,
(4)
T
Ω
gives the transmitted specific intensity due to
theincidentradiationattheboundaryand
NP
operator whose elements depend on the optical distances
between the grid-points. Thus, the mean intensity at the
grid-point ”
i
”wouldbe:
J
i
=Λ
i
,
1
S
1
+
...
+Λ
i
,
i
−
1
S
i
−
1
+Λ
i
,
i
S
i
+
Λ
Ω
isaNP
×
(5)
+Λ
i
,
i+1
S
i+1
+
...
+Λ
i
,
NP
S
c
NP
+T
i
.
Proc. Symposium
‘The Promise of the Herschel Space Observatory’
12–15 December 2000, Toledo, Spain
ESA SP-460, July 2001, eds. G.L.Pilbratt, J.Cernicharo, A.M.Heras, T.Prusti, & R. Harris
The well-known Λ
where J
old
where
262
J. Trujillo Bueno
)
/
2(withthe
sum applied to all the directions of the numerical angular
quadrature)and
a
,
b
and
c
aresimplysymbolsthat weuse
as a notational trick to indicate below whether we choose
the“old”orthe“new”valuesofthesourcefunction.Thus,
forinstance,theΛ-iterationmethodconsistsincalculating
J
i
choosing
a
=
b
=
c
=old,whichgivesJ
i
=J
old
InthislastexpressionΛ
ij
=
(Λ
i
ij
+Λ
out
ij
a formal solution and S
new
k
as dictated by Eqs. (10) and
(7).
4) Go to the next point
k
+ 1 and repeat what has
just been indicated in the previous point until arriving
to the other “boundary point”. Having reached this stage
iniciate againthe same process,but choosingnow as “first
point
i
= 1” the above-mentioned “boundary point”.
The result of what we have just indicated is a pure GS
iteration. Actually, after an
incoming
and
outgoing
pass
we get two GS iterations! A SOR iteration is obtained by
doing the corrections as follows:
δ
S
SO
i
=
ωδ
S
G
i
,
(11)
where
ω
is a parameter with an optimal value between 1
and 2 that can be found easily (see Trujillo Bueno and
Fabiani Bendicho, 1995).
Figure 1 shows an example of the convergence rate of
alltheseiterativemethodsappliedto aNLTEpolarization
transferproblem in a stellarmodel atmospherepermeated
byaconstantmagneticfieldthatproducesaZeemansplit-
tingof3Dopplerwidths.Wepointoutthatthecomputing
time
per iteration
is similar for all these methods and that
matrix inversions are not performed. Thus, the important
point to keep in mind is that our implementation of the
GS method is 4 times faster than Jacobi(i.e. than the ALI
method on which most NLTE codes are based on), while
our SOR method for radiative transfer applications is 10
times faster.
i
as
indicated in Eq. (3).
The Jacobi method, known in the RT literature as the
OAB method (from Olson, Auer and Buchler, 1986), and
on which most NLTE codes are based on, is found by
choosing
a
=
c
= old, but
b
= new, which yields
J
i
=J
old
i
+Λ
ii
δ
S
i
(6)
In fact, using this expression instead of J
ol
i
in Eq. (3)
one finds that the resulting Jacobi iterative scheme is
S
new
i
i
+Λ
ii
(S
new
i
−
S
ol
i
)=J
old
=S
old
i
+
δ
S
i
,
(7)
with the correction
δ
S
i
=
[(1
−
)J
old
i
+
B
i
−
S
ol
i
]
,
(8)
[1
−
(1
−
)Λ
ii
]
where Λ
ii
is the diagonal element of the Λ-operator as-
sociated to the spatial grid-point “
i
”. Note that the cor-
rectioncorrespondingtothe slowlyconvergentΛ-iteration
method is given by Eq. (8), but with Λ
ii
=0.
A superior type of iterative schemes are the Gauss-
Seidel (GS) based methods of Trujillo Bueno and Fabi-
ani Bendicho (1995), which can also be suitably general-
ized to the polarization transfer case (cf. Trujillo Bueno
& Landi Degl’Innocenti, 1996; Trujillo Bueno and Manso
Sainz, 1999). This type of iterative schemes are obtained
by choosing
c
=oldand
a
=
b
= new. This yields
J
i
=J
old and new
i
+Λ
ii
δ
S
i
,
(9)
is the mean intensity calculated using
the “new” values of the source function at grid-points
1,2,...,
i−
1and the “old” values at points
i, i
+1
,i
+2,....,
NP. The iterative correction is given by
δ
S
GS
i
=
[(1
−
)J
old and new
i
+
B
i
−
S
ol
i
]
(10)
[1
−
(1
−
)Λ
ii
]
It is important to clarify the meaning of this last equa-
tion:
1) First, at point
i
= 1 (which we can freely be as-
signed to any of the two boundaries of the medium under
consideration) use “old” source function values to calcu-
late J
1
via a formal solution. Apply Eqs. (10) and (7) to
calculate S
ne
1
.
2) Go to the next point
i
=2anduseS
new
1
and the
Figure 1. The variation of the maximum relative change ver-
sus the iteration number for several types of iterative meth-
ods applied to the NLTE Zeeman line transfer problem dis-
cussed by Trujillo Bueno and Landi Degl’Innocenti (1996).
The NLTE parameter
=10
−
4
. Dotted line: the
Λ
-iteration
method. Dashed line: the Jacobi-based ALI method. Solid line:
the GS-based method. Dashed-dotted line: the SOR method.
at points
j
=2
,
3
, ...,
NP
to get J
2
via a formal solution. Apply Eqs. (10) and (7)
to calculate S
ne
2
.
3)Gotothenextspatialpoint
k
andusethepreviously
obtained“new”sourcefunctionvaluesat
j
=1
,
2
, ..., k−
For pedagogical reasons we have chosen here a NLTE
linear
problem in order to explain in simple terms our
GS-based iterative schemes. The generalization to the full
non-linear
multilevel problem can be carried out as in-
dicated in the Appendix of the paper by Trujillo Bueno
1,
but the still “old” ones at
j
=
k, k
+1
, ...,
NP to get J
k
via
where J
old and new
i
“old” source-function values S
old
j
Radiative Transfer for the FIRST ERA
263
(1999). The critical point is always to remember that the
approximations one introduces for achieving the required
linearityof the statisticalequilibriumequationsateachit-
erative step should treat adequately the
coupling
between
transitions and the
non-locality
of the problem (see Socas-
Navarro & Trujillo Bueno, 1997).
3. The non-linear multigrid RT method
Our GS and SOR radiative transfer methods are based,
like the Jacobi-like ALI method, on the idea of operator
splitting. Therefore, all of them are characterized by a
convergent rate which
decreases
as the resolution of the
spatial grid is increased. As a result, if NP is the number
of grid-points in a computational box of
fixed
dimensions,
the computing time or computational work (
MG
)required
by the three previous iterative methods to yield the self-
consistent atomic (or molecular) level populations scales
withNPasfollows(cf.TrujilloBueno&FabianiBendicho,
1995):
–
Jacobi-based ALI method
→W∼
NP
2
–
Our GS-based RT method
→W∼
NP
2
/
4
W
–
Our SOR RT method
→W∼
NP
√
NP
Is there any suitable RT multilevel method character-
izedby
Figure 2. Variation with the grid-spacing
∆
z of the maximum
eigenvalue of the iteration operator corresponding to several
multilevel iterative schemes. The MG symbol refers to our non-
linear multigrid code, while MUGA to our multilevel GS-based
code. MALI refers to the Jacobi-based multilevel ALI method
of Rybicki & Hummer, (1991).
NP?Thiswouldbe ofgreatinterestfor3DRT
with multilevel atoms where NP
4. Formal solvers for RT applications
10
6
. The answer is af-
firmative. This has been worked out by Fabiani Bendicho,
Trujillo Bueno and Auer (1997) who considered the appli-
cation of the
non-linear
multigrid method (see Hackbush,
1985) to multilevel RT.
Theiterativeschemeofthenon-linearmultigridmethod
is composed of two parts: a
smoothing
one where a small
number of GS-based iterations on the desired
finest
grid
are used to get rid of the high-frequency spatial compo-
nents of the errorin the currentestimate, and a correction
obtainedfromthesolutionofanerrorequationina
coarser
grid. With our non-linear multigrid code the total com-
putational work scales simply as NP, although it must be
saidthatthecomputingtimeperiterationisabout4times
larger than that required by the Λ
∼
TheformalsolutionroutinesofourNLTEcodes(forunpo-
larized or polarized radiation and for atomic or molecular
species) are based on improvements and generalizations
of the
short-characteristics
(SC) technique introduced by
Kunasz & Auer (1988). Let us recall it briefly indicating
also our generalizationto 3D radiative transfer and to the
case of polarized radiation.
The scalar RT equation for the specific intensity is
dI
ν
d
s
=
χ
ν
(S
ν
−
I
ν
)
,
(12)
iteration method.
In order to compare the convergence rate of all these
iterative methods we present in Fig. 2 an estimate of
the maximum eigenvalues (
ρ
) of the corresponding iter-
ation operator, which controls the convergence properties
of such iterative schemes. The knowledge of this maxi-
mum eigenvalue (
ρ
) is useful because errors decrease as
ρ
itr
,where“
itr
” isthe iterativestep.As it canbe notedin
Fig.2theconvergencerateofboth, theMALI andMUGA
schemesdecreaseswhenthespatialresolutionofthegridis
improved,whilethemaximumeigenvalueofour
non-linear
multigrid method is always very small (
ρ ∼
−
where
s
is the geometric distance along the ray propagat-
ing in a certain direction in a 3D medium,
χ
ν
is the total
opacity and S
ν
the source function.
PointOisthegrid-pointofinterestatwhichonewishes
to calculate the specific intensity I
O
, for a given frequency
(
ν
) and a direction (
). Point M is the the intersection
point with the grid-plane that one finds when moving
along
Ω
.Atthis
upwind
point the specific intensity I
M
(for the same frequency and angle) is known from pre-
vious steps. In a similar way, point P is the intersection
point with the grid-plane that one encounters when mov-
ing along
−
Ω
. We also introduce the optical depths along
the ray between points M and O (∆
τ
M
) and between
points O and P (∆
τ
P
). From the formal solution of the
previous transfer equation one finds that
I
O
=I
M
e
−
∆
τ
M
+
∆
τ
M
0
Ω
0
.
1) and
in-
sensitive
to the grid-size. A maximum eigenvalue
ρ
=0
.
1
means that the error decreases by one order of magnitude
eachtimeweperformaniteration!Thisexplainsthat,typ-
ically, two multigrid iterations are su/cient to reach the
self-consistent solution for the atomic level populations.
S(
t
)e
−
(∆
τ
M
−
t
)
dt,
(13)
with the optical depth variable measured from M to O.
W∼
264
The integral of this equationcanbe solved
analytically
by integrating along the
short-characteristics
MO assum-
ing that the source function S(
t
)varies
parabolical ly
along
M,O and P. The result reads:
I
O
=I
M
e
−
∆
τ
M
+Ψ
M
S
M
+Ψ
O
S
O
+Ψ
P
S
P
,
(14)
The main changes when going to 3D imposing hori-
zontal periodic boundary conditions lie in the interpola-
tion. We have assumed that I
M
is known but, in most
cases, the M-point (like the point P) will not be a grid-
point of the chosen 3D spatial grid. The intensity at this
M-point has to be calculated by interpolating from the
available information at the nine surrounding grid-points,
as we must also do for obtaining the opacities and source
functions at M and P. Parabolic interpolation can how-
ever generate spurious negative intensities if the spatial
variation of the physical quantities is not well resolved
by the spatial grid. This happens, for instance, if one
tries to simulate the propagation of a beam in vacuum
using a three dimensional grid. To avoid these problems
we have improved the 1D monotonic interpolation strat-
egy of Auer and Paletou (1994), and generalized it to the
two-dimensional parabolic interpolation that is required
for 3D RT calculations with multilevel atomic models (see
Fabiani Bendicho & Trujillo Bueno, 1999).
where Ψ
X
(with X either M, O or P) are given in terms of
the quantities ∆
τ
M
and ∆
τ
P
that we evaluate numerically
byassumingthatln(
χ
)varieslinearlywiththegeometrical
depth,
χ
being the opacity.
If the interest lies in the generation and transfer of
po-
larized
radiation in magnetized astrophysical plasmas (cf.
Trujillo Bueno and Landi Degl’Innocenti, 1996; Trujillo
Bueno,1999;2001)the situationisabitmorecomplicated
because, instead of having to solve the previous RT equa-
tionforthespecificintensity,onehastosolve,ingeneral,a
vectorialtransfer equationfor the
four
Stokesparameters.
Forexample, forthe standardcaseofpolarizationinduced
by the Zeeman effect, the Stokes-vector at the grid-point
Ois
I
M
+
∆
τ
M
I
O
=
O
(0
,
∆
τ
M
)
0
O
(
t,
∆
τ
M
)
S
(
t
)
dt,
(15)
5. Conclusions
The RT methods presented here are especially attractive
because of their direct applicability to a variety of compli-
cated RT problems of astrophysical interest. We empha-
size that their convergence rate are extremely high, that
they do not require the construction and the inversion of
any large matrix and that the computing time per itera-
tion is very small.
4
Mueller matrix of the atmospheric slab between
t
and
∆
τ
M
). In general, this evolution operator does not have
aneasyanalyticalexpressionandthe integralofthe previ-
ousequationcannotbesolvedanalytically.However,ifthe
4
O
(
t,
∆
τ
M
) is the evolution operator (i.e. the 4
×
4 absorption matrix conmutes between depth points M
and O (e.g. because one assumes the absorption matrix to
be constant between M and O and equal to its true value
at the middle point) the evolution operator reduces then
to an expression given by the exponential of the absorp-
tion matrix. The integral of Eq. (15) can then be solved
analytically provided that the source function vector
References
is
assumed to vary parabolically along M, O and P. Our for-
mal solution method for Zeeman line transfer is precisely
based on this idea and the results reads:
S
Auer, L., Fabiani Bendicho, P., & Trujillo Bueno, J., 1994,
A&A, 292, 599
Auer, L., & Paletou, F.., 1994, A&A, 285, 675
Fabiani Bendicho, P., & Trujillo Bueno, J., 1999, in
Solar Po-
larization
, K. N. Nagendra & J. O. Stenflo (eds.), Kluwer
Academic Publishers, p. 219-230.
Fabiani Bendicho, P., Trujillo Bueno, J., & Auer, L.,1994,
A&A, 324, 161
Hackbush, W., 1985,
Multi-Grid Methods and Applications
,
Springer Verlag, Berlin
Kunasz, P., & Auer, L. H., 1988, JQRST, 39, 67
Landi Degl’Innocenti, E., & Landi Degl’Innocenti, M., 1995,
Solar Physics, 97, 239
Olson, G. L., Auer, L., & Buchler, J. R., 1986, JQRST, 35, 431
Rybicki, G. B., & Hummer, D. G., A&A, 245, 171
Socas Navarro, H., & Trujillo Bueno, J., 1999, ApJ, 516, 436
Trujillo Bueno, J., 1999, in
Solar Polarization
,K.N.Nagendra
& J. O. Stenflo (eds.), Kluwer Academic Publishers, p. 73-
96
Trujillo Bueno, J., 2001, in
Advanced Solar Polarimetry
,M.
Sigwarth (ed.) ASP Conf. Series, in press
Trujillo Bueno, J., & Fabiani Bendicho, P., 1995, ApJ, 455, 646
Trujillo Bueno, J., & Landi Degl’Innocenti, E., 1996, Solar
Physics, 164, 135
Trujillo Bueno, J., & Manso Sainz, R., 1999, ApJ, 516, 436
I
O
=
O
(0
,
∆
τ
M
)
I
M
+
Ψ
M
S
M
+
Ψ
O
S
O
+
Ψ
P
S
P
,
(16)
4 matrices
which are given in terms of ∆
τ
M
and ∆
τ
P
,intermsof
the inverse of the absorption matrix and in terms of the
analytical expression given by Landi Degl’Innocenti and
Landi Degl’Innocenti (1985)for the evolutionoperatorfor
the case in which the absorption matrix is assumed to be
constant between the spatial grid points. An alternative
formal solver of the Stokes-vector transfer equation suit-
able for NLTE applications is the one given by Eqs. (1-4)
of Socas-Navarro, Trujillo Bueno & Ruiz Cobo (2000).
The application of these formal solution methods in
1D is straightforward. The generalization to 2D geome-
tries with horizontal periodic boundary conditions is sub-
stantially more complicated. A suitable strategy has been
described by Auer, Fabiani Bendicho and Trujillo Bueno
(1994).
Ψ
X
(with X either M, O or P) are 4
×
where
×
where
[ Pobierz całość w formacie PDF ]