- Email: [email protected]
A numerical solution of the line transfer equation in a spherically symmetric atmosphere
Journal of Quantitative Spectroscopy & Radiative Transfer 62 (1999) 599}608
Note
A numerical solution of the line transfer equation in a spherically symmetric atmosphere D. Mohan Rao Indian Institute of Astrophysics, Koramangala, Bangalore 560 034, India Received 14 September 1998
Abstract
A numerical solution based on "nite di!erences is proposed for solving the line transfer equation in a spherically symmetric atmosphere. We discretize the curvature term in the equation by the backwarddi!erence scheme for the outward-going ray from the atmosphere and by the forward-di!erence scheme for the backward-going ray. As the radiation "eld is highly anisotropic in the far-wings of the line, the accuracy of the method is improved by computing the source function "rst using the numerical scheme for a limited set of rays. The intensity in the line is then computed with the known source function using the formal solution of the transfer equation with a large grid of rays. The solution is compared with the published results using di!erent representation of the transfer equation and is found to be quite accurate. Another salient feature of the method is that one can enforce the non-negativity of the solution quite easily. 1999 Elsevier Science Ltd. All rights reserved.
1. Introduction For modeling the supergiants and giants one approximates the atmospheres of these stars as spherically symmetric. Several workers in this "eld proposed diverse techniques for solving the radiative transfer equation in spherical geometry. Here we summarize in brief some of the computational techniques which are based on "nite di!erences. We will not comment on the techniques like operator perturbation techniques [2] and integral equation methods [3]. Hummer and Rybicki [4] developed a method based on iterating the space-dependent Eddington factor for the conservative grey case. It can be [5, 6] extended to the case of linearly polarized radiation and for the problem of line formation. This technique takes into account the outward peaking of the radiation "eld that occurs in extended spherical systems. In the presence of radial velocity "elds [6] the transfer equation is transformed to (p, z) coordinate system where p"r sin h, z"r cos h, r is the radial coordinate and h is the angle made by the ray relative to the radial direction. The equation is di!erenced along the rays and the resulting set of linear equations 0022-4073/99/$ - see front matter 1999 Elsevier Science Ltd. All rights reserved PII: S 0 0 2 2 - 4 0 7 3 ( 9 8 ) 0 0 1 2 3 - X
600
D. Mohan Rao / Journal of Quantitative Spectroscopy & Radiative Transfer 62 (1999) 599}608
are solved using Rybicki's [7] method. All these techniques are highly accurate and found to be stable. However, it may be di$cult to generalize them to problems with complicated phase functions and anisotropic velocity "elds. In contrast to the above approach, Peraiah and Grant [8] integrate the monochromatic transfer equation over a radial and angular mesh and the resulting equations are solved using the discrete space theory of radiative transfer. The solution of the transfer equation is non-negative provided one satis"es the conditions on radial and angular size of the mesh. The same technique is applied [9] for the problem of spectral line formation. However, in the line wings one "nds a stringent condition for the solution to be non-negative. The negative intensities are avoided in some situations by calculating the source functions from the solution and using the formal solution of the transfer equation [10]. An integral operator technique was developed [11] within the framework of discrete theory of radiative transfer theory. The resulting algorithm is quite complicated and involves a large number of matrices. In this paper we develop a method which follows the same approach [8] but di!ers in the treatment of the discretization of the curvature term. We derive a stability condition which can be enforced easily. The method can be easily extended to the cases with arbitrary phase functions. In Section 2., we brie#y describe the method and the results are given in Section. 3.
2. A description of the numerical method 2.1. The transfer equation and the boundary conditions The transfer equation in spherical symmetry under the assumption of complete redistribution is given by $k
*I!(r, k, x) 1!k *I!(r, k, x) $ #p(r) (x)I!(r, k, x)"p(r) (x)S(r, $k) *k r *r
(1)
where I is the speci"c intensity of the radiation at a spatial position r, with frequency x in the direction given by h"cos\ k. The absorption coe$cient is denoted by p(r) and (x) is the pro"le function. The sign$refers to the radiation #owing outward and inward through the atmosphere, respectively. The S(r, k) is the line source function given by 1 P($k, k ) (x ) I(r, k , x ) dx dk #e(r)B(r), (2) 2 \ \ where e(r) is the probability that the photon is destroyed by collisional de-excitation and the P(k, k ) is the phase function. We assume the medium to have an outer radius R surrounding a core of radius R ("1) and consider two cases: (1) the core is hollow; (2) the core is perfectly absorbing and emits radiation in a speci"ed direction. In the "rst case S(r, $k)"(1!e(r))
I(R , #k, x)"I(R , !k, x) (re#ection at the lower boundary). In the second case I(R , #k, x)"1.0.
D. Mohan Rao / Journal of Quantitative Spectroscopy & Radiative Transfer 62 (1999) 599}608
601
In both the cases I(R, !k, x)"0, (i.e. incident radiation is zero at the top of the atmosphere). 2.2. Discretization of the scattering integral We discretize the angular integral in the scattering term by a Gaussian quadrature formula of order J and the frequency integral by a Trapezoidal formula of order I. The scattering integral can be written as
\ \ I(r, k , x )P($k, k ) (x ) dx dk " \ I(r, k , x )P($k, k ) (x ) dx dk # \ I(r, !k , x )P($k, k ) (x ) dx dk
( ' ( ' (3) " P($k, k )I(r, k , x )c a # P($k, !k ) I(r, !k , x ) c a , H H G H G H H G H G H G H G where k and c are the Gaussian roots and weights. x and a are the trapezoidal points.We use the H H G G notation I>(r, k , x)"I(r, #k , x) and I\(r, k , x)"I(r, !k , x). H H H H
(4)
2.3. Finite-diwerence scheme Here we shall describe a "nite-di!erence representation of Eq. (1). The medium is divided into n shells r "R'r '2'r "1. If we integrate Eq. (1) over a shell [r , r ] and apply ,> L L> a backward-di!erence scheme for the curvature term in the case of the outward-going ray k (I>(r , k , x )!I>(r , k , x )) H L H G L> H G (1!k) (I>r , k , x )!I>(r ,k ,x) H L> H G L> H\ G *r # L r (k !k ) L> H H\ #p(r ) (x )I>(r , k , x )*r "p(r ) (x )(S(r , k )*r L> G L> H G L L> G L>/ H L for j"1, 2 , J, i"1, 2 , I.
(5)
To compute I>(r, k , x) we need to know I(r, k "0, x). The intensity values I>(r, k , x), H j"2, 2 , J depend sensitively on this quantity and is di$cult to determine. Hence we adopt the following approach. We set (1!k) *I!(r, k, x) "0 for q (x)510 *k r since the radiation "eld is nearly isotropic if the line optical depth is large.
(6)
602
D. Mohan Rao / Journal of Quantitative Spectroscopy & Radiative Transfer 62 (1999) 599}608
By applying the linear interpolation formulas I!(r , k , x) between 0 and k , and eliminating L>/ I!(r , 0, x), we get L> k , k , x)!I\(r , k , x)) I>(r , k , x)!I\(r , k , x)" (I>(r L> L> L> L>/ k for 10\4q (x)(10. (7) In the very outer layers of the medium, for q (x)(10\, we apply the above formula at r instead L of r . L> If we apply a forward-di!erence scheme for the inward-going ray we get !k (I\(r , k , x )!I\(r , k , x )) H L H G L> H G , k , x )!I\(r ,k,x) (1!k) (I\(r L> H> G L> H G *r H ! L (k !k ) r H> H L> #p(r ) (x ) I\(r , k , x )*r "p(r ) (x )(S(r , !k )*r L> G L> H G L L> G L> H L for j"1, 2 , J, i"1, 2, 2 , I, where *r "r !r . We will approximate I!(r ) by the &&diamond rule'' L L L> L> I!(r )"0.5(I!(r )#I!(r )). L> L> L Now we shall write Eqs. (5) and (8) in a matrix form. Using Eq. (2), We have
(8)
(9)
M[;>!;> ]#0.5o">[;> #;>]#0.5*q [;> #;>]# L L> L> L L L> L "0.25(1!e )'>>(;> #;>)*q #0.25(1!e )'>\(;\ #;\)*q L> L> L L L> L> L L #e B *q (10) L> L> L !M[;\!;\ ]#0.5o"\[;\ #;\]#0.5*q [;\ #;\] L L> L> L L L> L "0.25(1!e )'\\(;\ #;\)*q #0.25(1!e )'\>(;> #;>)*q L> L> L L L> L> L L #e B *q (11) L> L> L where o"*r /r and *q "p *r . The quantities e ,p and B are replaced by L L L L> L L> L> L> their cell-averages. At j"J, we set the curvature term to zero as k K1. The matrices M, '!!, ;! ( L and the curvature matrices "! are de"ned in the Appendix A. We can arrange these equations in a canonical form
;> t(n, n#1) r(n#1, n) ;> &> L " L> # L> ;\ r(n, n#1) t(n#1, n) ;\ &\ L> L L> for 1(2(2(n"N,
(12)
where r and t are the re#ection and transmission matrices. They are given in the Appendix B. These equations can be solved using the algorithm given in [8] for obtaining the emergent and internal radiation "eld.
D. Mohan Rao / Journal of Quantitative Spectroscopy & Radiative Transfer 62 (1999) 599}608
603
2.4. Non-negativity of the numerical solution In this section we shall derive a condition for the re#ection and transmission matrices r and t to be non-negative. For this we must have *>, *\, S>>, S\\5O, where O is a null matrix (see Appendix B). From Appendix A, we note that the o!-diagonal elements of (*>) \ are negative and the diagonal elements of the matrix are positive. The matrix (*>)\"(a ) is diagonally dominant I I (see [1]), i.e., I'( a 5 "a " I I I I I I $I
(13)
The matrix is also strictly diagonally dominant as the strict inequality is valid for all the elements except at j"1. It is also irreducible as all the o!-diagonal entries of any row or column of the matrix are not zero. Hence from Varga [1] (Theorem 3.11, Corollary 1), *>'O. For S>>'O, we must have
(1!k) *r H *r !0.5 (x )*q '0 N0.5 k !0.5 (1!k)4k . H L G L H r r L> L>
Since we use logarithmic spacing in *q, this condition is naturally satis"ed in the outer layers (*rK10\, *qK10\). As the radiation "eld is nearly isotropic in the deeper layers, we can violate the condition there without much trouble. Similar arguments can be used for *\, S\\ and the stability condition remains the same. One should note that the approach for discretization of the curvature term in [8] gives diagonal and o!-diagonal elements for ">, "\ with large positive and negative signs. Hence, it is di$cult to obtain a non-negative condition in the far wings of the line at the limb (kK0).
3. Results and discussion For checking the accuracy of the method, we have chosen two cases. In both the cases, only isotropic scattering is assumed (i.e., p(k, k )"1). Case 1: We adopted the model used by Avertt and Loeser [3]. We chose a spherical medium with the outer radius r "30 and the inner absorbing core with r "1. The absorption coe$cient a is varied as r\. The radial optical thickness from r to r is 4. The , source function for monochromatic scattering is given by S "aJ #(1!a)S?, G G G
(14)
where a is the scattering coe$cient. In Table 1 the mean intensity J values are tabulated along with those given in [3] for a scattering coe$cient a"1/2 and S?"1.0. The di!erences in the values G between ours and [3] are of the order 7% in the outer layers of the atmosphere. Case 2: We adopted the same model as that of Hummer and Kunasz [6]. We have chosen in all the cases p(r)J1/r, B(r)"1 and the total optical thickness ¹"1000. The Gaussian pro"le
604
D. Mohan Rao / Journal of Quantitative Spectroscopy & Radiative Transfer 62 (1999) 599}608
Table 1 Mean intensity J at selected values of r, from di!erent calculations r
J (Mihalas et al.)
J (Rogers)
J (Avrett et al.)
J (our results)
30 29.6 29 28 26 24 20 16 12 9 6 4 2 1
0.0636 0.0671 0.0716 0.0784 0.092 0.106 0.141 0.182 0.26 0.345 0.493 0.638 0.864 0.964
0.0638 0.0674 0.0716 0.0785 0.0921 0.107 0.14 0.187 0.255 0.336 0.468 0.615 0.849 0.961
0.0642 0.0676 0.0721 0.0791 0.0928 0.108 0.141 0.19 0.261 0.34 0.476 0.627 0.859 0.986
0.0686 0.0717 0.0761 0.0855 0.098 0.113 0.147 0.196 0.262 0.341 0.472 0.619 0.853 0.988
Fig. 1. The emergent #uxes are plotted against the frequency. The curves are labeled by the radius of extension of the medium.
D. Mohan Rao / Journal of Quantitative Spectroscopy & Radiative Transfer 62 (1999) 599}608
605
Table 2 Comparison of source functions for ¹"10, e"10\ R
30 100 300
S(q"1)
S(q"10)
S(q"100)
Our work
Kunasz et al.
Our work
Kunasz et al.
our work
Kunasz et al.
4.14E-04 2.8E-04 2.32E-04
3.9E-04 2.6E-04 2.2E-04
1.84E-03 1.56E-03 1.47E-03
1.7E-03 1.5E-03 1.5E-03
1.23E-02 1.2E-02 1.19E-02
1.2E-02 1.2E-02 1.2E-02
Table 3 Comparison of source functions for ¹"10, e"10\ R
30 100 300
S(q"1)
S(q"10)
S(q"100)
Our work
Kunasz et al.
Our work
Kunasz et al.
our work
Kunasz et al.
1.4E-03 7.14E-04 4.17E-04
1.4E-03 6.7E-04 3.9E-04
4.4E-03 2.61E-03 1.9E-03
4.4E-03 2.5E-03 1.8E-03
2.67E-02 2.30E-02 2.18E-02
2.5E-02 2.1E-02 2.0E-02
function is assumed, i.e., 1
(x)" e\V. (n
(15)
The source function for this case is given by Eq. (2). We choose 15 Gaussian angles and 50 radial shells. When the medium is highly extended (R53), we are unable to match the line pro"les in the wings (x52) with that of [6]. The reason is that the radiation "eld is highly anisotropic in the wings throughout the medium whereas in the core it is anisotropic only in the outer layers. Hence, we computed the radiation "eld using 15 angles and determined the source function. With the known source function , we can use the same algorithm [Eqs. (5) and (8)] for the formal solution of the transfer equation for computing the radiation "eld ray by ray. This does not involve the necessity of using the large matrices as only the intensity information at the previous angle is needed. We use 100 Gaussian angles to compute the line pro"les and plot the #ux pro"les in Fig. 1. We "nd that our "gure matches with that of Fig. 8 in [6]. The values of source functions are tabulated in Tables 2 and 3 along with those of [6]. The maximum di!erence in the results are of the order of 10% at S(q"100) for the cases in which R"100 and 300. We also performed the calculations using 4, 8, 15, 20 angles in the main algorithm and used the formal solution of radiative transfer equation. We found that at least 10}15 angles have to be used in the main algorithm to obtain the correct source function.
606
D. Mohan Rao / Journal of Quantitative Spectroscopy & Radiative Transfer 62 (1999) 599}608
4. Conclusions We have demonstrated the accuracy of the numerical method based on a simple "nite-di!erence representation of the curvature term in the transfer equation. We compared the values of the source function obtained applying this method and those of the other methods using a di!erent form of the transfer equation. The di!erences are not signi"cant. Though we use only isotropic scattering in this paper for simplicity, the extension of this technique to the problem of anisotropic phase function is straightforward.
Appendix A By writing M "[a ]"k d for j"1, 2 , J, j "1, 2 , J H HH ( HH
(A.1)
">"[l> ] ( HH (1!k) H l "0 for j"j "1" for j "j, j"2, 2 , J HH (k !k ) H H\ (1!k) H for j "j!1, j"2, 2 , J. "! (k !k ) H H\
(A.2)
"\"[m ], ( HH (1!k) H for j "j, j"1, 2 , J!1 m " HH (k !k ) H> H (1!k) H for j "j#1, j"1, 2 , J!1, "! (k !k ) H> H
(A.3)
"[ ]"( )d , I II II
" (x , k ) and k"j#(i!1)J, 14k4K"IJ, I G H
(A.4)
;!"[;> , ;> 2 , ;> 2 , ;> ], L L L G L ' L ;! "[I!(q , k , x ), I!(q , k , x ), 2 , I!(q , k , x )], G L L G L G L ( G
(A.5)
D. Mohan Rao / Journal of Quantitative Spectroscopy & Radiative Transfer 62 (1999) 599}608
607
where t denotes the transpose.
M> O ( O M> ( M>" O
"\ (
"\"
\
O
,
M> ; ( ' '
O "\ \ ( \
"> (
">"
O "> (
\
,
"> ; ( ' '
,
"\ ; ( ' '
'!!"[a ], k"j#(i!1)J, k "j #(i !1)J, II
(A.6)
where a "p($k , $k ) (x ) (x )a c , for i, i "1, 2, 2 , I, for j, j "1, 2 , J. H H G G G H II Appendix B *>"[M>#0.5o">!0.25(1!e)*q '>>#0.5*q ]\, L L S>>"[M>!0.5o">#0.25(1!e)*q '>>!0.5*q ], L L S>\"0.25(1!e)*q '>\, L *\"[M\#0.5o"\!0.25(1!e)*q '\\#0.5*q ]\, L L S\\"[M\!0.5o"\#0.25(1!e)*q '\\!0.5*q ], L L S\>"0.25(1!e)*q '\>, L r>\"*>S>\,
r\>"*\S\>,
t>"[I!r>\r\>]\,
t\"[I!r\>r>\]\,
t(n, n#1)"t>[*>S>>#r>\r\>], t(n#1, n)"t\[*\S\\#r\>r>\], r(n#1, n)"t>[r>\#r>\*\S\\], r(n, n#1)"t\[r\>#r\>*>S>>], >
"*q et>[*>r>\#*\B] L L> \ "*q et\[*>r\>#*>B]. L> L
(B.1)
608
D. Mohan Rao / Journal of Quantitative Spectroscopy & Radiative Transfer 62 (1999) 599}608
Acknowledgements The author likes to thank B.A.Varghese for his help in the preparation of the manuscript and Dr. K.E. Rangarajan for his invaluable suggestions.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11]
Varga RS. Matrix iterative analysis. Englewood Cli!s, NJ: Prentice-Hall, Inc., 1962;23}96. Kubat J. Astron Astrophys 1996;305:255}264. Avrett EH, Loeser R. Methods in radiative transfer. Cambridge: Cambridge University Press, 1984, 341. Hummer DG, Rybicki GB. Mon Not R Astron Soc 1971;152:1. Cassinelli FP, Hummer DG. Mon Not R Astron Soc 1971;153:9. Kunasz PB, Hummer DG. Mon Not R Astron Soc 1974;166:19. Rybicki GB. JQSRT 1971;11:589. Peraiah A, Grant IP. J Inst Maths Appl 1973;12:75}90. Grant IP, Peraiah A. Mon Not R Astron Soc 1972;160:239. Peraiah A. Kodaikanal Obs Bull Ser A bf 1978;2:180}183. Peraiah A, Varghese BA. Astrophys J 1985;290:411.
Note
A numerical solution of the line transfer equation in a spherically symmetric atmosphere D. Mohan Rao Indian Institute of Astrophysics, Koramangala, Bangalore 560 034, India Received 14 September 1998
Abstract
A numerical solution based on "nite di!erences is proposed for solving the line transfer equation in a spherically symmetric atmosphere. We discretize the curvature term in the equation by the backwarddi!erence scheme for the outward-going ray from the atmosphere and by the forward-di!erence scheme for the backward-going ray. As the radiation "eld is highly anisotropic in the far-wings of the line, the accuracy of the method is improved by computing the source function "rst using the numerical scheme for a limited set of rays. The intensity in the line is then computed with the known source function using the formal solution of the transfer equation with a large grid of rays. The solution is compared with the published results using di!erent representation of the transfer equation and is found to be quite accurate. Another salient feature of the method is that one can enforce the non-negativity of the solution quite easily. 1999 Elsevier Science Ltd. All rights reserved.
1. Introduction For modeling the supergiants and giants one approximates the atmospheres of these stars as spherically symmetric. Several workers in this "eld proposed diverse techniques for solving the radiative transfer equation in spherical geometry. Here we summarize in brief some of the computational techniques which are based on "nite di!erences. We will not comment on the techniques like operator perturbation techniques [2] and integral equation methods [3]. Hummer and Rybicki [4] developed a method based on iterating the space-dependent Eddington factor for the conservative grey case. It can be [5, 6] extended to the case of linearly polarized radiation and for the problem of line formation. This technique takes into account the outward peaking of the radiation "eld that occurs in extended spherical systems. In the presence of radial velocity "elds [6] the transfer equation is transformed to (p, z) coordinate system where p"r sin h, z"r cos h, r is the radial coordinate and h is the angle made by the ray relative to the radial direction. The equation is di!erenced along the rays and the resulting set of linear equations 0022-4073/99/$ - see front matter 1999 Elsevier Science Ltd. All rights reserved PII: S 0 0 2 2 - 4 0 7 3 ( 9 8 ) 0 0 1 2 3 - X
600
D. Mohan Rao / Journal of Quantitative Spectroscopy & Radiative Transfer 62 (1999) 599}608
are solved using Rybicki's [7] method. All these techniques are highly accurate and found to be stable. However, it may be di$cult to generalize them to problems with complicated phase functions and anisotropic velocity "elds. In contrast to the above approach, Peraiah and Grant [8] integrate the monochromatic transfer equation over a radial and angular mesh and the resulting equations are solved using the discrete space theory of radiative transfer. The solution of the transfer equation is non-negative provided one satis"es the conditions on radial and angular size of the mesh. The same technique is applied [9] for the problem of spectral line formation. However, in the line wings one "nds a stringent condition for the solution to be non-negative. The negative intensities are avoided in some situations by calculating the source functions from the solution and using the formal solution of the transfer equation [10]. An integral operator technique was developed [11] within the framework of discrete theory of radiative transfer theory. The resulting algorithm is quite complicated and involves a large number of matrices. In this paper we develop a method which follows the same approach [8] but di!ers in the treatment of the discretization of the curvature term. We derive a stability condition which can be enforced easily. The method can be easily extended to the cases with arbitrary phase functions. In Section 2., we brie#y describe the method and the results are given in Section. 3.
2. A description of the numerical method 2.1. The transfer equation and the boundary conditions The transfer equation in spherical symmetry under the assumption of complete redistribution is given by $k
*I!(r, k, x) 1!k *I!(r, k, x) $ #p(r) (x)I!(r, k, x)"p(r) (x)S(r, $k) *k r *r
(1)
where I is the speci"c intensity of the radiation at a spatial position r, with frequency x in the direction given by h"cos\ k. The absorption coe$cient is denoted by p(r) and (x) is the pro"le function. The sign$refers to the radiation #owing outward and inward through the atmosphere, respectively. The S(r, k) is the line source function given by 1 P($k, k ) (x ) I(r, k , x ) dx dk #e(r)B(r), (2) 2 \ \ where e(r) is the probability that the photon is destroyed by collisional de-excitation and the P(k, k ) is the phase function. We assume the medium to have an outer radius R surrounding a core of radius R ("1) and consider two cases: (1) the core is hollow; (2) the core is perfectly absorbing and emits radiation in a speci"ed direction. In the "rst case S(r, $k)"(1!e(r))
I(R , #k, x)"I(R , !k, x) (re#ection at the lower boundary). In the second case I(R , #k, x)"1.0.
D. Mohan Rao / Journal of Quantitative Spectroscopy & Radiative Transfer 62 (1999) 599}608
601
In both the cases I(R, !k, x)"0, (i.e. incident radiation is zero at the top of the atmosphere). 2.2. Discretization of the scattering integral We discretize the angular integral in the scattering term by a Gaussian quadrature formula of order J and the frequency integral by a Trapezoidal formula of order I. The scattering integral can be written as
\ \ I(r, k , x )P($k, k ) (x ) dx dk " \ I(r, k , x )P($k, k ) (x ) dx dk # \ I(r, !k , x )P($k, k ) (x ) dx dk
( ' ( ' (3) " P($k, k )I(r, k , x )c a # P($k, !k ) I(r, !k , x ) c a , H H G H G H H G H G H G H G where k and c are the Gaussian roots and weights. x and a are the trapezoidal points.We use the H H G G notation I>(r, k , x)"I(r, #k , x) and I\(r, k , x)"I(r, !k , x). H H H H
(4)
2.3. Finite-diwerence scheme Here we shall describe a "nite-di!erence representation of Eq. (1). The medium is divided into n shells r "R'r '2'r "1. If we integrate Eq. (1) over a shell [r , r ] and apply ,> L L> a backward-di!erence scheme for the curvature term in the case of the outward-going ray k (I>(r , k , x )!I>(r , k , x )) H L H G L> H G (1!k) (I>r , k , x )!I>(r ,k ,x) H L> H G L> H\ G *r # L r (k !k ) L> H H\ #p(r ) (x )I>(r , k , x )*r "p(r ) (x )(S(r , k )*r L> G L> H G L L> G L>/ H L for j"1, 2 , J, i"1, 2 , I.
(5)
To compute I>(r, k , x) we need to know I(r, k "0, x). The intensity values I>(r, k , x), H j"2, 2 , J depend sensitively on this quantity and is di$cult to determine. Hence we adopt the following approach. We set (1!k) *I!(r, k, x) "0 for q (x)510 *k r since the radiation "eld is nearly isotropic if the line optical depth is large.
(6)
602
D. Mohan Rao / Journal of Quantitative Spectroscopy & Radiative Transfer 62 (1999) 599}608
By applying the linear interpolation formulas I!(r , k , x) between 0 and k , and eliminating L>/ I!(r , 0, x), we get L> k , k , x)!I\(r , k , x)) I>(r , k , x)!I\(r , k , x)" (I>(r L> L> L> L>/ k for 10\4q (x)(10. (7) In the very outer layers of the medium, for q (x)(10\, we apply the above formula at r instead L of r . L> If we apply a forward-di!erence scheme for the inward-going ray we get !k (I\(r , k , x )!I\(r , k , x )) H L H G L> H G , k , x )!I\(r ,k,x) (1!k) (I\(r L> H> G L> H G *r H ! L (k !k ) r H> H L> #p(r ) (x ) I\(r , k , x )*r "p(r ) (x )(S(r , !k )*r L> G L> H G L L> G L> H L for j"1, 2 , J, i"1, 2, 2 , I, where *r "r !r . We will approximate I!(r ) by the &&diamond rule'' L L L> L> I!(r )"0.5(I!(r )#I!(r )). L> L> L Now we shall write Eqs. (5) and (8) in a matrix form. Using Eq. (2), We have
(8)
(9)
M[;>!;> ]#0.5o">[;> #;>]#0.5*q [;> #;>]# L L> L> L L L> L "0.25(1!e )'>>(;> #;>)*q #0.25(1!e )'>\(;\ #;\)*q L> L> L L L> L> L L #e B *q (10) L> L> L !M[;\!;\ ]#0.5o"\[;\ #;\]#0.5*q [;\ #;\] L L> L> L L L> L "0.25(1!e )'\\(;\ #;\)*q #0.25(1!e )'\>(;> #;>)*q L> L> L L L> L> L L #e B *q (11) L> L> L where o"*r /r and *q "p *r . The quantities e ,p and B are replaced by L L L L> L L> L> L> their cell-averages. At j"J, we set the curvature term to zero as k K1. The matrices M, '!!, ;! ( L and the curvature matrices "! are de"ned in the Appendix A. We can arrange these equations in a canonical form
;> t(n, n#1) r(n#1, n) ;> &> L " L> # L> ;\ r(n, n#1) t(n#1, n) ;\ &\ L> L L> for 1(2(2(n"N,
(12)
where r and t are the re#ection and transmission matrices. They are given in the Appendix B. These equations can be solved using the algorithm given in [8] for obtaining the emergent and internal radiation "eld.
D. Mohan Rao / Journal of Quantitative Spectroscopy & Radiative Transfer 62 (1999) 599}608
603
2.4. Non-negativity of the numerical solution In this section we shall derive a condition for the re#ection and transmission matrices r and t to be non-negative. For this we must have *>, *\, S>>, S\\5O, where O is a null matrix (see Appendix B). From Appendix A, we note that the o!-diagonal elements of (*>) \ are negative and the diagonal elements of the matrix are positive. The matrix (*>)\"(a ) is diagonally dominant I I (see [1]), i.e., I'( a 5 "a " I I I I I I $I
(13)
The matrix is also strictly diagonally dominant as the strict inequality is valid for all the elements except at j"1. It is also irreducible as all the o!-diagonal entries of any row or column of the matrix are not zero. Hence from Varga [1] (Theorem 3.11, Corollary 1), *>'O. For S>>'O, we must have
(1!k) *r H *r !0.5 (x )*q '0 N0.5 k !0.5 (1!k)4k . H L G L H r r L> L>
Since we use logarithmic spacing in *q, this condition is naturally satis"ed in the outer layers (*rK10\, *qK10\). As the radiation "eld is nearly isotropic in the deeper layers, we can violate the condition there without much trouble. Similar arguments can be used for *\, S\\ and the stability condition remains the same. One should note that the approach for discretization of the curvature term in [8] gives diagonal and o!-diagonal elements for ">, "\ with large positive and negative signs. Hence, it is di$cult to obtain a non-negative condition in the far wings of the line at the limb (kK0).
3. Results and discussion For checking the accuracy of the method, we have chosen two cases. In both the cases, only isotropic scattering is assumed (i.e., p(k, k )"1). Case 1: We adopted the model used by Avertt and Loeser [3]. We chose a spherical medium with the outer radius r "30 and the inner absorbing core with r "1. The absorption coe$cient a is varied as r\. The radial optical thickness from r to r is 4. The , source function for monochromatic scattering is given by S "aJ #(1!a)S?, G G G
(14)
where a is the scattering coe$cient. In Table 1 the mean intensity J values are tabulated along with those given in [3] for a scattering coe$cient a"1/2 and S?"1.0. The di!erences in the values G between ours and [3] are of the order 7% in the outer layers of the atmosphere. Case 2: We adopted the same model as that of Hummer and Kunasz [6]. We have chosen in all the cases p(r)J1/r, B(r)"1 and the total optical thickness ¹"1000. The Gaussian pro"le
604
D. Mohan Rao / Journal of Quantitative Spectroscopy & Radiative Transfer 62 (1999) 599}608
Table 1 Mean intensity J at selected values of r, from di!erent calculations r
J (Mihalas et al.)
J (Rogers)
J (Avrett et al.)
J (our results)
30 29.6 29 28 26 24 20 16 12 9 6 4 2 1
0.0636 0.0671 0.0716 0.0784 0.092 0.106 0.141 0.182 0.26 0.345 0.493 0.638 0.864 0.964
0.0638 0.0674 0.0716 0.0785 0.0921 0.107 0.14 0.187 0.255 0.336 0.468 0.615 0.849 0.961
0.0642 0.0676 0.0721 0.0791 0.0928 0.108 0.141 0.19 0.261 0.34 0.476 0.627 0.859 0.986
0.0686 0.0717 0.0761 0.0855 0.098 0.113 0.147 0.196 0.262 0.341 0.472 0.619 0.853 0.988
Fig. 1. The emergent #uxes are plotted against the frequency. The curves are labeled by the radius of extension of the medium.
D. Mohan Rao / Journal of Quantitative Spectroscopy & Radiative Transfer 62 (1999) 599}608
605
Table 2 Comparison of source functions for ¹"10, e"10\ R
30 100 300
S(q"1)
S(q"10)
S(q"100)
Our work
Kunasz et al.
Our work
Kunasz et al.
our work
Kunasz et al.
4.14E-04 2.8E-04 2.32E-04
3.9E-04 2.6E-04 2.2E-04
1.84E-03 1.56E-03 1.47E-03
1.7E-03 1.5E-03 1.5E-03
1.23E-02 1.2E-02 1.19E-02
1.2E-02 1.2E-02 1.2E-02
Table 3 Comparison of source functions for ¹"10, e"10\ R
30 100 300
S(q"1)
S(q"10)
S(q"100)
Our work
Kunasz et al.
Our work
Kunasz et al.
our work
Kunasz et al.
1.4E-03 7.14E-04 4.17E-04
1.4E-03 6.7E-04 3.9E-04
4.4E-03 2.61E-03 1.9E-03
4.4E-03 2.5E-03 1.8E-03
2.67E-02 2.30E-02 2.18E-02
2.5E-02 2.1E-02 2.0E-02
function is assumed, i.e., 1
(x)" e\V. (n
(15)
The source function for this case is given by Eq. (2). We choose 15 Gaussian angles and 50 radial shells. When the medium is highly extended (R53), we are unable to match the line pro"les in the wings (x52) with that of [6]. The reason is that the radiation "eld is highly anisotropic in the wings throughout the medium whereas in the core it is anisotropic only in the outer layers. Hence, we computed the radiation "eld using 15 angles and determined the source function. With the known source function , we can use the same algorithm [Eqs. (5) and (8)] for the formal solution of the transfer equation for computing the radiation "eld ray by ray. This does not involve the necessity of using the large matrices as only the intensity information at the previous angle is needed. We use 100 Gaussian angles to compute the line pro"les and plot the #ux pro"les in Fig. 1. We "nd that our "gure matches with that of Fig. 8 in [6]. The values of source functions are tabulated in Tables 2 and 3 along with those of [6]. The maximum di!erence in the results are of the order of 10% at S(q"100) for the cases in which R"100 and 300. We also performed the calculations using 4, 8, 15, 20 angles in the main algorithm and used the formal solution of radiative transfer equation. We found that at least 10}15 angles have to be used in the main algorithm to obtain the correct source function.
606
D. Mohan Rao / Journal of Quantitative Spectroscopy & Radiative Transfer 62 (1999) 599}608
4. Conclusions We have demonstrated the accuracy of the numerical method based on a simple "nite-di!erence representation of the curvature term in the transfer equation. We compared the values of the source function obtained applying this method and those of the other methods using a di!erent form of the transfer equation. The di!erences are not signi"cant. Though we use only isotropic scattering in this paper for simplicity, the extension of this technique to the problem of anisotropic phase function is straightforward.
Appendix A By writing M "[a ]"k d for j"1, 2 , J, j "1, 2 , J H HH ( HH
(A.1)
">"[l> ] ( HH (1!k) H l "0 for j"j "1" for j "j, j"2, 2 , J HH (k !k ) H H\ (1!k) H for j "j!1, j"2, 2 , J. "! (k !k ) H H\
(A.2)
"\"[m ], ( HH (1!k) H for j "j, j"1, 2 , J!1 m " HH (k !k ) H> H (1!k) H for j "j#1, j"1, 2 , J!1, "! (k !k ) H> H
(A.3)
"[ ]"( )d , I II II
" (x , k ) and k"j#(i!1)J, 14k4K"IJ, I G H
(A.4)
;!"[;> , ;> 2 , ;> 2 , ;> ], L L L G L ' L ;! "[I!(q , k , x ), I!(q , k , x ), 2 , I!(q , k , x )], G L L G L G L ( G
(A.5)
D. Mohan Rao / Journal of Quantitative Spectroscopy & Radiative Transfer 62 (1999) 599}608
607
where t denotes the transpose.
M> O ( O M> ( M>" O
"\ (
"\"
\
O
,
M> ; ( ' '
O "\ \ ( \
"> (
">"
O "> (
\
,
"> ; ( ' '
,
"\ ; ( ' '
'!!"[a ], k"j#(i!1)J, k "j #(i !1)J, II
(A.6)
where a "p($k , $k ) (x ) (x )a c , for i, i "1, 2, 2 , I, for j, j "1, 2 , J. H H G G G H II Appendix B *>"[M>#0.5o">!0.25(1!e)*q '>>#0.5*q ]\, L L S>>"[M>!0.5o">#0.25(1!e)*q '>>!0.5*q ], L L S>\"0.25(1!e)*q '>\, L *\"[M\#0.5o"\!0.25(1!e)*q '\\#0.5*q ]\, L L S\\"[M\!0.5o"\#0.25(1!e)*q '\\!0.5*q ], L L S\>"0.25(1!e)*q '\>, L r>\"*>S>\,
r\>"*\S\>,
t>"[I!r>\r\>]\,
t\"[I!r\>r>\]\,
t(n, n#1)"t>[*>S>>#r>\r\>], t(n#1, n)"t\[*\S\\#r\>r>\], r(n#1, n)"t>[r>\#r>\*\S\\], r(n, n#1)"t\[r\>#r\>*>S>>], >
"*q et>[*>r>\#*\B] L L> \ "*q et\[*>r\>#*>B]. L> L
(B.1)
608
D. Mohan Rao / Journal of Quantitative Spectroscopy & Radiative Transfer 62 (1999) 599}608
Acknowledgements The author likes to thank B.A.Varghese for his help in the preparation of the manuscript and Dr. K.E. Rangarajan for his invaluable suggestions.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11]
Varga RS. Matrix iterative analysis. Englewood Cli!s, NJ: Prentice-Hall, Inc., 1962;23}96. Kubat J. Astron Astrophys 1996;305:255}264. Avrett EH, Loeser R. Methods in radiative transfer. Cambridge: Cambridge University Press, 1984, 341. Hummer DG, Rybicki GB. Mon Not R Astron Soc 1971;152:1. Cassinelli FP, Hummer DG. Mon Not R Astron Soc 1971;153:9. Kunasz PB, Hummer DG. Mon Not R Astron Soc 1974;166:19. Rybicki GB. JQSRT 1971;11:589. Peraiah A, Grant IP. J Inst Maths Appl 1973;12:75}90. Grant IP, Peraiah A. Mon Not R Astron Soc 1972;160:239. Peraiah A. Kodaikanal Obs Bull Ser A bf 1978;2:180}183. Peraiah A, Varghese BA. Astrophys J 1985;290:411.