Numerical solution of the radiation transport equation in disk geometry

J. Quant. Spectrosc. R&at. Transfer Vol. 31, No. 6, pp. 565-580, Printed in Great Britain. All rights reserved

NUMERICAL TRANSPORT

1987 Copyright

SOLUTION EQUATION

0

0022-4073/87 1987 Pergamon

$3.00 + 0.00 Journals Ltd

OF THE RADIATION IN DISK GEOMETRY

GEORGEF. SPAGNAJR and CHUN MING LEUNG Department of Physics, Rensselaer Polytechnic Institute, Troy, NY 12180-3590, U.S.A. (Received

8 August 1986)

Abstract-We describe an efficient numerical method for solving the problem of radiation transport in a dusty medium with two dimensional (2-D) disk geometry. It is a generalization of the one-dimensional quasi-diffusion method in which the transport equation is cast in diffusion form and then solved as a boundary value problem. The method should Lx applicable to a variety of astronomical sources, the dynamics of which are angular-momentum dominated and hence not accurately treated by spherical geometry, e.g. protoplanetary nebulae, circumstellar disks, interstellar molecular clouds, accretion disks, and disk galaxies. The computational procedure and practical considerations for implementing the method are described in detail. To illustrate the effects of 2-D radiation transport, we compare some model results (dust temperature distributions and i.r. flux spectra) for externally heated, interstellar dust clouds with spherically symmetric and disk geometry.

1.

INTRODUCTION

The presence of dust grains plays an important role in the thermodynamics and evolution of many astronomical objects. The association of this dust with stars, nebulae, interstellar molecular clouds, and nuclei of some galaxies is well established from observations in the infrared. Since thermal emission by dust is believed to be the dominant mechanism for producing the i.r. radiation, any interpretation of i.r. observations based on modeling of the sources will invariably involve studies of radiation transport through dust. A key theoretical problem not addressed by existing models of radiation transport in dust is the geometry of the i.r. source. Although there is growing observational and theoretical evidence for a large number of disk-shaped or toroidal objects of astrophysical interest (e.g. circumstellar disks, protoplanetary disks, protostellar accretion disks, bipolar molecular flows, and disk galaxies), the majority of radiative transfer models currently in use invoke the assumption of spherical geometry. This assumption is made because spherically symmetric geometry is the only one-dimensional (1-D) geometry which accounts for the finite dimensions of a system in all directions. Establishing the source geometry can provide severe constraints on the origin, dynamics, and properties of the i.r. sources. Hence it is crucial to solve the problem of radiation transport in a dusty medium with two-dimensional (2-D) disk geometry and apply the results to interpreting astronomical observations in the infrared. Progress toward the solution of multi-dimensional transport problems, though substantial, has been largely concerned with 2-D planar geometry,’ with emphasis on the interpretation of fine spatial structure on the solar surface. Literature on multi-dimensional transfer calculations has been reviewed by Jones and SkumanichZ and a survey of work on heat transfer problems in twoand three-dimensional geometries with engineering applications is given by Shih.3 Very few attempts have been made to solve the problem of scattering, absorption and emission by dust grains in a medium with disk geometry. Recently, Lefevre et aL4 have performed Monte-Carlo simulations of ellipsoidal dust shells around cool stars. Their approach may lack flexibility since they report that the angular and spatial grids used are not well suited for modeling highly flattened ellipsoids. Here we describe a self-consistent, numerical solution to the problem of radiation transport in a dusty medium with 2-D disk geometry. The method is based on the use of the so-called variable Eddington factors in modeling stellar atmospheres. 54 This technique, better described as a quasi-diffusion method (QDM), has been applied9 to a dusty medium with 1-D spherical geometry. The QDM proved to be very flexible and efficient, requiring only very modest computing effort while treating the transport problem (e.g. anisotropic scattering, radiation field anisotropy) accurately. The basic ideas of the QDM as applied to a dusty medium with 1-D geometries have 565

566

GEORGEF. SPAGNAJR and Cm

MING LEUNG

been described in detail by Leung9,” and a FORTRAN computer code applying the method to spherically symmetric geometry has been documented by Spagna and Leung.” Recently, this computer code has been updated and generalized to include the other two 1-D geometries (planar and cylindrical).‘2 As a result, the method/computer program has been adopted by many authors and applied successfully to a wide variety of astronomical problems directly related to infrared observations. In this paper, we generalize the QDM fqr a dusty medium with 2-D disk geometry. In this method the transport equation is first cast as a set of moment equations of the radiation field. The system of equations is then placed in quasi-diffusion form by introducing a set of auxilliary functions which allow the moment equations to be closed. The first set of functions, called anisotropy factors, forms the elements of the anisotropy tensor (also called the Eddington tensor in the literature of stellar atmospheres13) and measures the anisotropy of the radiation field. The second is the configuration function which depends on an integral involving the anisotropy factors. This function plays a similar role as the sphericality factor, first introduced by Auer14 in modeling extended stellar atmospheres. By first assuming the radiation field to be isotropic (diffusion approximation), we initially fix the anisotropy factors and the configuration function. We then linearize and solve the coupled nonlinear system of moment and energy balance equations using the Newton-Raphson iteration method to estimate the mean intensity and the dust temperature distribution, i.e. the source function. From the source function, we solve for the angular distribution of the radiation field by a ray-tracing technique and use the solution to update the anisotropy factors and the configuration function. This iteration procedure is repeated until the solution converges. Physically, the success of the QDM rests upon the fact that the temperature distribution depends strongly on the radiation energy density and only weakly on the anisotropy of the radiation field. This weak coupling between the thermal structure and radiation field anisotropy allows the problem to be separated into two distinct steps: solution of the moment and energy-balance equations for a given radiation field anisotropy, followed by solution of the ray equations to determine the angular distribution of the radiation field for a given source function. The outline of this paper is as follows. In Section 2, we derive the moment equations and the associated boundary conditions. The derivation of the ray equations, the associated boundary equations, as well as the evaluation of the anisotropy factors and the configuration function are discussed in Section 3. In Sections 4-6, we describe respectively the numerical solutions of the moment equations, the ray equations, and the system of linearized moment and energy-balance equations. In Section 7 we outline the practical considerations for implementing the numerical method. In the last section, we discuss the application of the method to two specific problems: the first is a conservative gray case which admits an analytic solution and the second is an example to determine the temperature structure and emergent flux for an externally heated dust cloud. The physics of radiation transport in disk geometry and its observational implications will be discussed more fully in an upcoming paper.i5 2. MOMENT

EQUATIONS

AND BOUNDARY

CONDITIONS

In our formulation below, we shall use cylindrical coordinates (r, z, 8,+) and omit the frequency subscript for all frequency-dependent quantities. Unless otherwise stated, we shall consider only one frequency at a time. Consider an element dr at some point (T, z) along a photon trajectory in the disk. The polar angle made with the figure (symmetry) axis is 8 while the azimuthal angle to the outward radial direction is 4, as shown in Fig. 1. The time-independent equation of radiation transport for the monochromatic specific intensity Z(r, z, 8,4) in the case of anisotropic scattering and isotropic source may be written as dZ/ds = cos 8 (aZ/az) + sin 8 cos #J (al/&) - (l/r) sin 8 sin 4 (al/a&) 2n

53

Z’p(cos a) sin 8’ de’ d4’], (1) s @=Os e=o where rr’ and es are respectively the volume absorption and scattering coefficients, and B is the Planck function corresponding to the temperature of dust grains. The scattering phase function, = -

aqz

-

B) -

d[Z

-

(l/411)

Disk geometry solution of the radiation transport equation

561

Fig. 1. The direction of photon propagation, denoted by s, and its relation to the spatial (r, z) and angular (0,4) variables in disk geometry are shown.

p(cos u), determines the probability of scattering between two directions separated by an angle a and satisfies the following normalization condition: p(cos a) da = 1,

(1/4x) s

(2)

where cosa=cos8cosO’+

sin0sin8’(cos4cos$‘+sind,sin$‘)

(3)

defines the scattering angle a and the integration is taken over all solid angles. The scattering phase function is usually expanded in a series of Legendre polynomials.‘6 For linear anisotropic scattering, p(cos a) is then approximated as (1 + 3g cos a), where g is the scattering asymmetry parameter. For isotropic scattering g = 0, for forward scattering g = 1, and for backward scattering g = - 1. In the following derivation and for all the models reported here we have assumed isotropic scattering for simplicity. We now define the a-component of the nth angular moment of the specific intensity Z as M”, = (1/47r)

ZQ dR. s

We may decompose the direction unit vector R into its components directions in the local coordinate system as Q = sin

8 cos 4,

(4) along the three orthogonal

(5)

R, = cos 8,

(6)

R, = sin 8 sin 4,

(7)

568

GEORGEF. SPAGNAJR and CHUN MING LEUNG

and the integration

over all solid angles is Sd*=~~~~~~=~sinBdBdm.

(8)

We identify the zeroth moment (MO) of Z as the mean intensity .Z and the three components of the first moment (Ma, a = r, z, 4) as the Eddington fluxes along the three orthogonal directions in the local coordinate system: H,, H,, and Hb. Of the nine components for the second moment which we call the K-integrals, we only need four (AZ;, /I = rr, zz, C#JC#, rz, using ther notation Qg, = f&Q) for our formulation. The K-integrals measure the anisotropy of the radiation field and are proportional to the elements of the radiation pressure tensor or the electromagnetic stress tensor.13 Specifically, we have J(r, z) = (1/4x)

P2n rn

2n fW,

z>

=

(1/4x)

H,(r, z) = (1/4x)

z) =

(1/4x)

Mr,

z)

=

(1/4x)

Ur,

z>

=

(1/4x)

J I$=0J e=o

z)

=

(1/4x)

z)

=

(1/4x)

Z(r, z, 8,4)

sin28 cos 4 de d$,

Z(r, z, 0,4) sin 0

cos 8

(10)

de d&

(11)

n

I$=0 e=o

(12)

Z(r, Z, 8,#) sin30 COS’~de d$,

(13)

Z(r, Z, 8, 4) cos2 &I sin 8 de d4,

(14)

Z(r, z, &#I) sin38 sin24 de d$,

(1%

II

f o=o s e=o 2R

Z(r, z, tl,t$) sin28 sin 4 de d4,

11

s ,#I=0s e=o 2n

II

s .$=of e=o 2n

KAr,

(9)

ss 2n

&&,

s

sin 0 de d#,

x

f #=O 0=0

2n

H&r,

Z(r, z, 44)

J #=OJ e=o

11

Z(r, Z,&$)COS

8

sin28

cos tp de

d#.

s +5=0s e=o

(16)

If we integrate equation (1) over all solid angles, i.e. take its zeroth angular moment, we have (aH,Pz)

+

W)[W%Wrl

= v - xJ,

(17)

where the volume emission coefficient r~is given by ‘1 = o”J + a”B and x = 8 + & is the volume extinction coefficient. Note that the scattering terms in equation (1) cancel out in the above zeroth moment equation. The first r- and z-moments of equation (1) are obtained by integrating over all solid angles with the respective weights R, and Q: (&JJz)

+ (KMr) G%Pz)

Next, we define the anisotropy radiation field as

+ (llr)(K, +

-&)

= xH,,

Wr)kWKzY~rl = - xK.

factors (also known as variable Eddington

(18) (19) factorssp6J3) of the

_L = WJ,

(20)

_fiz= KzlJ,

(21)

f& = K&J,

(22)

_fL= LIJ.

(23)

Disk geometry solution of the radiation transport equation Table 1. Behavior of the anisotropy Radiation field

factors

Anisotropy

Isotropic Radially streaming Axially streaming Anisotropic

569

factors

f”=L,=f&=f._L=O L = 1, Lz = & = L = 0 L,,=Lf,=&=.L=O f,+f,+joo=L_Cz+O

It can be shown” that, in general, the sum xr + f, + f44= 1. Furthermore, for a completely isotropic radiation field (diffusion approximation), fr, = f, = f+,+ = l/3, and f,= = 0. For radiation streaming along a given direction, the conjugate anisotropy factor approaches unity. When radiation is peaking sideward relative to a given direction, the conjugate anisotropy factor drops below l/3 and approaches zero. The behavior of the anisotropy factors under various conditions is summarized in Table 1. For the extreme cases of complete isotropy and of radiation streaming either radially or axially, frz = 0, but for an anisotropic radiation field, frzis nonzero in general. After introducing the following integrating factor (which we call the configuration function 5) into equation (18), (24) the first-order moment equations, equations (18)-(19), may be rewritten as

(25) P(_AzJY~zl+ w)[wLJ>/~~l=

-

xffz*

(26)

Substituting ZZ,and ZZ,from the above equations into equation (17) yields the combined moment equation (CME) of radiation transport for disk geometry:

+ (l/r)(a/ar){(r/~)[d(f,~)/azi

+ (ri~5)[a(f,,5Z)/ari)

= XJ - V. (27)

Although the CME has been derived assuming isotropic scattering, equation (27) is also valid for general anisotropic scattering provided the expressions for 1 and q are modified as discussed in Leung.9 If 2 and q are known (i.e. knowing the grain temperature distribution) and with the assumption of diffusion approximation, the CME becomes linear in J and can be solved directly without resorting to any linearization scheme. This step is used to start the Newton-Raphson iteration procedure to determine the radiation field consistent with an initial guess of the dust temperature distribution. Boundary equations are now developed from the first-order moment equations [equations (25) and (26)] by specifying the net fluxes. Because of the assumed axial symmetry we need to solve the CME only for a single quadrant of the (r-z)-plane, bounded by the figure axis (r = 0), the midplane (z = 0), the disk edge (r = R), and the disk surface (z = Z). At each boundary we decompose each of the two Eddington fluxes (ZZ,and ZZ,) into its receding (H-) and approaching (H+ ) components (with respect to an external observer): n/2 ’ Z+ (r, z, 8,4) sin*8 cos $J de d+, H,+ (r, z) = (1/4x) (28) s q5=-n/2 s e=o 3n/2 (r, z) =

(1/4x)

H,+ (r, z) =

(1/4x)

H;

(1/4x)

H;

2n

(29)

n/2

z+ (r, Z, O,C#J)sin 8 cos 8 de d& s +=o s L9=0 2n

(r, 2) =

n

~= ,2 B DZ-(I,z,e,~)sin28cosddBd~, II s = s

(30)

1

I- (r, Z, 8, 4) sin 8 cos 8 de d$. 5 +=o s e=n/2

(31)

GEORGE F. SPAGNA JR and CHUNMING

570

LEUNG

Physically, I+ corresponds to radiation approaching an external observer and I- to radiation receding from the observer. Usually ZZ- is known at the outer boundary while H+ is given at the inner boundary. The incident radiation field at either boundary is usually assumed to be isotropic and is given as an input parameter in the model. By introducing an appropriate boundary factor, defined as the ratio of the sum of the flux components to the mean intensity at the given boundary, we may express the net fluxes in equations (25) and (26) as (for c1= r or z) Ha = faBJ - 2H;

(at the outer boundary B),

(32)

H, = 2H,+ - AcJ

(at the inner boundary C),

(33)

H, = 0

(along figure axis r = 0 and along midplane z = 0),

(34)

where we have defined the following boundary factors:

z)l/JW, z),

(35)

Z)l/J(r, Z),

(36)

_L = [H,‘(R, z) + H; (R, frB = [HZ (r, 2) + H; (r,

f< = [Hz (rc, z) + H; (rc, z)llJ(r,, fi

= [HZ (r, 2,) + H; (r,

z),

(37)

z,)llJ(r, z,),

(38)

Here r, and z, denote the (r-z)-coordinates of the central core at the inner boundary. For externally heated models with no inner core, rc = z, = 0. We now express x and q in terms of the dust number density nd, the grain radius a, as well as the absorption and scattering efficiency factors (Qabsand Q,) from the Mie scattering theory, viz. X = %]Qabsna* +

Qsca~~*l,

(39)

v = %[Q,~a*J + Q,~,~~‘W)I, where at each point (r, z) in the medium, the equilibrium grain temperature by the condition of energy balance

w-9 T(r, z) is determined

s

[q - xJ] dv = (W/471).

(41)

Here, W represents the amount of energy lost by the grains through nonradiative processes, e.g. energy transfer to the gas component through gas-grain collisions or molecule formation. If we neglect such processes, the resulting equation of radiative equilibrium reduces to the usual form [q - xJ] dv = Qabsrra2[B(T) - J] dv = 0. (42) s f If we integrate the zeroth moment equation [equation (1711 over all frequencies and apply the constraint of radiative equilibrium, as given by equation (42), we obtain the condition for the conservation of integrated net flux in disk geometry by integrating over r and z I

I r’H,(r’,

S[f

0

z) dr’ + r

H,(r, z’) dz’ s0

1

dv = H*(r, z) = constant,

(43)

where H*(r, z) is the frequency integrated net flux through an imaginary cylindrical surface of radius r and half-thickness z, concentric with the disk shaped cloud. For externally heated models without a central source H* = 0, so that the integrated net flux vanishes for all such surfaces inside the cloud. On the other hand, for models with a central heat source (rc # 0, z, # 0), Ha = L */( 16rr*), where L * is the luminosity of the central source. A self-consistent solution to the problem of grain temperature distribution and radiation transport then amounts to solving the CME [equation (27)] in conjunction with the boundary equations [equations (25), (26), (32)-(34)], subject to the constraint of radiative equilibrium which is given by equation (42). In addition, if the diffusion approximation is not used, we must also determine the angular distribution of the radiation field to evaluate the anisotropy factors and the configuration function in the CME. The flux conservation condition in equation (43) is used as an important check on the accuracy and self-consistency of the numerical solution.

Disk geometry solution of the radiation transport equation

571

Fig. 2. Shown here are the definitions of variables y and z in an impact plane defined by an impact parameter [ in the solution of the ray equation. Also shown is their relation to a photon trajectory s lying in the plane and the angle 0.

3. THE

RAY

EQUATIONS

AND

BOUNDARY

CONDITIONS

We determine the angular distribution of the radiation field by a ray-tracing technique, i.e. the transport equation for each frequency is integrated along parallel rays parametrized by their distance of closest approach (impact parameter) to the figure axis and making an angle 8 with the z-direction. In particular, consider a plane (which we call an impact plane) parallel to the figure axis and defined by an impact parameter c, as shown in Fig. 2. On this impact plane we establish a (y-z)-coordinate system with z = 0 at the midplane and y = 0 at the symmetry line of the plane (defined by r = [). In terms of the disk coordinates (r, Cp)and the impact parameter c, y is given by y = r cos 4 = (r* - [2)1/2.

(W

Physically, y represents the projection of a ray onto a plane parallel to the midplane, with the ray constrained to lie in an impact plane. In the numerical solution we assign to c the same discrete set of values as r. For a ray confined to an impact plane and making an angle 8 with the direction along the figure axis, we may write the transport equation either as a 1-D equation along the ray or as a 2-D equation in the (y-z)-coordinates as dZ/ds = ,@Z/~Z) + &YZ/Q) = ‘1 - xZ,

(45)

where ZA= cos 8 and /I = sin 8. For a fixed angle 8, this is the first-order ray equation for transport in a disk. We can recast this equation into second-order form by separating Z into its positive and negative components of ~1,i.e. I+ (1 3 ~12 0) and I- (0 > p 2 - l), and using their symmetric (intensity-like) and antisymmetric (flux-like) averages along the ray u = (1/2)[z+ + I-], v = (1/2)[z+ -1-l.

(46)

512

GEORGEF. SPAGNAJR and CHUN MING LEUNG

Physically, I+ corresponds to radiation approaching the observer and Z- to radiation receding from the observer. By alternately adding the subtracting the transport equation (45) written respectively for I+ and I-, we obtain the following coupled equations for u and v: (dulds) = /@/a~)

+ /?(au/ay) = - XV,

(dv/ds) = /@/a~)

+ /?(&J/@) = n - xu.

(48) (49)

Eliminating v between the two equations we arrive at the second-order (CRE) either in 1-D form as (l/X)(d’uldS’)

+

combined ray equation

Ml/xYdsl(du/ds) = P - rl

(50)

or, in 2-D form, as (1/X)[CL2(~2u/dz2)+ pz(a2U/ay’) + 2~~(a%/ayaz)] +

wwlxY~4 + BwlxY~Y1)Mwaz) + mu/ay)l= xu - v.

(51)

Note that unlike the CME, the scattering terms on the right-hand side of the CRE do not cancel since q depends also on J. Given the source function and opacity (i.e. given x and n), the CRE becomes linear in u and can easily be solved for u from which we can evaluate the angular moments (the mean intensity ,Z and the K-integrals), and hence the anisotropy factors and the configuration function. Unlike the case of the CME, in general the ray equation must be solved for all four quadrants of the (y-z)-plane for each angle 8 and for each impact parameter [. For the special case of rays parallel to the midplane (0 = 90”) the 2-D form of the ray equation is reduced to 1-D form (in y only) and is solved only for the first quadrant. The boundary conditions for the CRE are developed in complete analogy with the CME by using the first-order ray equation in equation (45) and expressing the flux-like quantity u in terms of u, I+, and I- as follows: du/ds = /@/a~)

+ B(au/+)

= x(Z- - u),

(52)

du/ds = ~(~u/~z) + /~(Ju/$Y) = x(u -I+),

(53)

For p = 0, we impose the additional symmetry condition at y = 0 dujdr =

ss +=o e=o

u(r,z,&$)sin@dt)d$,

(55)

L(r,z)=(l/471)

2n n U(I, s S$=os .9=0

Z, 0, 4)

sin%

_L(r,z) = (1/4x)

2n R u(r, s ,$=os 8=0

z, 8, 4)

cos2q5

Z, 8, C#J)

sin38 sin24 d0 d$/J(r,

ss s2nsn 2n

f&(r,

z) =

(1/4x)

Mr,z)=(l/4~)

COS’+

de d$/J(r,

sin 8 de d4/J(r,

z),

(56)

z),

(57)

z),

(58)

de d@/.Z(r, z),

(59)

n

u(r,

4=0 0=0

~(t, I, 8,&) cos 8 sin%

cos ip

$=O e=o n/z

fie(K z) = (1/4x)

* u(R, s #=-n/Z s e=o

Z, l-l,+)

sin28 cos

4

dtI d$/J(R,

z),

(60)

Disk geometry solution of the radiation transport equation

&(r, Z) = (1/4x)

2n n/Zu(r, Z, 8,4) sin 8 cos 8 dB d4/J(r, s o=o s 8=0

ss

513

Z),

(61)

42

.&(ro z) = (1/4x)

II u(ro

+=-n/2

j&r, z,) = (1/4x)

Z,

8,q5)

sin28 cos 4 de d4/J(r,,

z),

(62)

e=o

2n n’2 u (r, s 4=0 s e=o

4. SOLUTION

z,,

8, 4) sin 0 cos 8 de d4 /J(r,

OF THE

MOMENT

z,).

(63)

EQUATIONS

As stated earlier, for a given grain temperature distribution T(r, z) and under the diffusion approximation (or with some assumed values for the anisotropy factors), the CME [equation (2711 becomes linear in J and hence can be solved directly using the method of finite differences. Because of the assumed axial symmetry we need to solve the CME only for a single (first) quadrant of the (r-z)-plane. We discretize the first quadrant of the (r-z)-plane by a 2-D grid characterized by a set of r-points [{ri}, ranging from r, = 0 to rNR= R, where R is the disk radius and NR is the number of r-points] and a set of z-points [{zj}, ranging from z, = 0 to z,,,z = Z, where Z is the disk half-thickness and NZ is the number of z-points]. The CME is then discretized using a finite difference scheme based on the trapezoidal rule and is replaced by a set of algebraic equations. In these finite difference equations, the geometric factors, which depend on the grid spacings (involving only ri and Zj) and independent of frequency, are calculated only once and stored as a set of spatial quadrature weights. The boundary equations are obtained by incorporating the corresponding first-order equations into a second-order Taylor series at each boundary in a manner described by Fox.‘~ Compared to using first-order boundary equations, this method provides improved accuracy and numerical stability. Thus, at the disk surface (z = Z), the product fzJ is expanded in a second-order Taylor series in z as follows: firJIz_ti

= frrJI, -

WWJ>/~zllz

+ (Az2/2)[a2(f,,J)/az211~

(64)

The first-order term comes from the first-order moment equation for H, [equation (26)] and the second order term from the CME [equation (27)]. Similarly, at the disk edge (r = R), the product (f,,(J)is expanded in a Taylor series in r, with first-order term coming from the first-order moment equation for H, [equation (25)] and the second-order term from the CME. Thus at the disk edge we have

_L~JIR-~r=f,t&

- ~~WXO/~~llR + (Ar2/2)[a2(f,,5J)lar211R.

(65)

(z = 0), we havethe additional symmetry condition that all first derivatives z must vanish. Therefore the Taylor series in z involves only the second-order term taken from the CME, with the first derivatives in z, and the cross-derivative terms in r and z set equal to zero. Similarly, at the figure axis (r = O), we require the first derivative in r, all cross derivatives in r and z, and the functionS, to vanish. This latter condition implies that the first z-derivative of JJ will vanish as well. For models with a central core at the inner boundary, a set of boundary equations similar to equations (64) and (65) would apply to the inner core with coordinates (rc, zJ. The complete set of finite difference equations for both boundary and interior points is then arranged in order of increasing zj, and within each j ordered by increasing ri. The resulting system of equations is block-tridiagonal in form and may be represented in general as

At the midplane

$=AjJi_,+BjJj+CjJj+,-Dj=O,

W)

forj=2,... , NZ - 1. Here the Aj, Bj, and Cj are (NR x NR) tridiagonal matrices, and the J’s and Dis are vectors with NR elements. This block-tridiagonal system of equations can easily be solved by the standard Gaussian elimination method.i9

574

GEORGE

F. SPAGNAJR and CHUN MINGLEUNG

5. SOLUTION

OF THE

RAY

EQUATIONS

In principle, the angular distribution of the radiation field can be determined by solving the CRE in either 1-D form [equation (50)] or 2-D form [equation (51)] on a series of impact planes. In practice, some caution must be exercised in choosing which form to use in the numerical solution by the method of finite differences, Specifically, a ray passing through a given grid point at an arbitrary angle is unlikely to intercept the grid at any other grid point, but will instead cross the grid at some intermediate points. For problems in which the characteristics of the medium and the source function do not vary too rapidly along the grid lines, the ray equations may be solved in its 2-D form [equation (51)] on a 2-D grid directly using the method of finite differences. In such cases, even though all derivatives are taken along grid lines (instead of along the rays), the formulation still provides sufficient accuracy and numerical stability. Furthermore, the overall storage needed for the program is somewhat smaller. Boundary equations for the 2-D representation of the CRE are developed in complete analogy with those for the CME. For example, at the disk face in the first quadrant, we expand u as a Taylor series in z. We substitute the finite difference form of equation (48) for the first derivative terms and for the second derivative terms we use the finite difference form of equation (51). For problems in which the opacity and source function vary rapidly with position, the 2-D formulation may present a computational difficulty in that interpolation of the solution along the grid lines may not give a good representation of the behavior of the functions u and v along a given ray. In this case, the 2-D formulation has been found to generate non-physical results for some combinations of model parameters. Acccordingly, we have also solved the CRE in 1-D form [equation (SO)] following the approach discussed by Mihalas et al.’ In this method we interpolate the solution of the 1-D CRE at each point on an impact plane onto the grid of surrounding points. This allows us to treat those variables whose behavior is well-described by interpolation along the grid lines (e.g. the opacity) and still treat the ray variables accurately along their trajectories. Specifically, at each point in the (y, z)-grid we trace a ray at an angle 0 which will not in general intercept the surrounding grid at any of the discrete grid points. We determine the forward and backward intervals for each point, dr + and d.s-, respectively, and discretize the ray equation along these ray segments using Lagrange (parabolic) interpolation. This finite difference equation is then interpolated onto the stencil of nine points (see Fig. 3) centered on (vi, zj), using parabolic interpolation. The procedure is carried out for each point of an impact plane. As shown in Mihalas

I

I

I

p

I

---

I -y

I

I

l

I

I

I

I

YI-1

yi

Y1+1

Fig. 3. In the solution of the ray equations on an impact plane, a ray through the grid-point 0 does not, in general, intersect any of the other grid-points defined by the intersections of the mesh lines. Furthermore, the spacing of points of intersections of the ray with the mesh lines is extremely irregular. To overcome this numerical problem, values of all variables at points hi and P are obtained by parabolic interpolation on the periphery of a nine-point stencil surrounding point 0 in a (y-z) impact plane.

515

Disk geometry solution of the radiation transport equation

et al.‘, this scheme retains the overall block-tridiagonal forms of the equations, preserving the accuracy of the interpolations of those variables best described along the grid while providing accurate interpolations for those variables best described along the ray. The weights for these interpolations are stored sequentially in a set of 1-D arrays, in order of their use by the program. Second-order boundary equations for the 1-D CRE are analogous to those in the 2-D case. One sided finite differences are written at the boundary, with derivatives expressed along the ray and with parabolic interpolation onto the interior grid points. The parabolic interpolation weights for the boundary are defined in exactly the same manner as for the interior, with the simplification that one need only consider forward or backward rays, but not both. As in the case of the CME, the finite difference weights for solving the CRE are computed only once and stored for subsequent use by the program. To determine the anisotropy and bounday factors, angular integrations are carried out by discretizing the angular variables and assigning quadrature weights for the various integrals. We choose a Gaussian quadrature set for the 8 integrations. A three-point quadrature on 0 from 0 to n is adequate, but we find that a five-point quadrature provides a higher accuracy and more detailed information on the angular variation of the emergent intensity. Clearly the 4 integration is more complicated. Not only does the discretized variable change in a non-uniform manner along a given ray for a given impact parameter, but the discrete values will be different for each impact plane. This means that a different quadrature set in 4 must be used at each r value in the grid for each impact plane. A simple scheme involving the trapezoidal rule is accordingly used, which provides sufficient accuracy and stability even near the axis where there are fewer points in the quadrature set. Solution to the ray equation is also used to determine the emergent intensity as a function of both position on the surface and viewing angle. 6. SOLUTION OF TEMPERATURE

THE LINEARIZED AND INTENSITY

EQUATIONS CORRECTIONS

FOR

A self-consistent determination of the grain temperature problem amounts to simultaneously solving the CME and the equation of radiative equilibrium. If the diffusion approximation is not used, the CRE must also be solved for each frequency. These coupled, nonlinear, partial differential equations are solved numerically by the method of finite differences using an iteration procedure similar to the Newton-Raphson method. Below we describe the details of the linearization scheme and the computational procedure. In applying the Newton-Raphson method, we may either completely linearize the system around the current values of the mean intensity J and the temperature T, or we may partially linearize the system in one variable (such as the temperature T) if J and T are not strongly coupled. For the grain temperature problem involving frequency-dependent opacities, Tin general depends only weakly on J. In such a case we can use a partial linearization method in which we only linearize the energy balance equation in T (and not in both T and J), with J obtained by solving the CME using the updated temperature. Thus for NF discrete frequencies, represented by the index n, the discretized form of the energy balance equation [equation (42)] is Ej = F WF, (Qabsxa*)n[Ji,n - Sl;, n=l

.I= 0,

(67)

where WF,, and (Qab~nu2),are respectively the frequency quadrature weight and the absorption cross section for the n th frequency, and BTj_ is the Planck function corresponding to temperature q. If we linearize this equation in the temperature T only (keeping J fixed), we obtain AT = -Ej

F WF, (Qabsm2),, (aBT&YTi> n=l

1 .

(68)

The mean intensities are then obtained by solving the CME using the updated temperature (q + A?). Thus we solve the problem by iterating among the CME, the linearized energy balance equation, and the CRE (to update the anisotropy factors), with most of the computing time spent in solving the ray equations. This scheme, though less elegant and general than complete linearization, reduces the overall computing effort considerably since it avoids the necessity of

GEORGEF. SPAGNAJR and Cm

516

MING LEUNG

inverting large full matrices which occur in the complete linearization approach as discussed below. This computational procedure is justified physically by the relatively weak dependence of the temperature on the mean intensity and details of the radiation field anisotropy. Experience shows this partial linearization method to be stable, with convergence usually within 10 iterations for reasonable initial guesses of the temperature distribution. We have successfully used this approach to construct all the models reported here. For problems in which the temperature and the radiation field are strongly coupled, the partial linearization method may converge very slowly, if at all. To overcome this difficulty, one must then use the complete linearization method in which the CME and the energy balance equation are both linearized and solved simultaneously. The linearized equations of our system of (NR x NZ) x (NE + 1) nonlinear difference equations, I;;;,(Ji,,,, rj) = 0 and ,!$(Ji,+, q) = 0 for j=l,..., NZandn=l,... , NF, are obtained by evaluating the partial derivatives of I;;;” and Ej with respect to each of the variables & and Tj. Thus from the discretized form of the CME [equation (66)], we can write its linearized form at each frequency n as Aj,nAJI-,~,+Bj,RJj,“+Cj,“AJi+l,,-D~,ATi=

-q,;;;,,

(69)

forn=l,... , NF. Here D;+ = aD,./aT,. This system of equations must be solved simultaneously with the equation of energy balance, which is now linearized in both T and J as f n=l

WFn(Q o~s~a2),,[AJj,n - (aBT,./aTj)ATj]

= - Ej.

(70)

The system of linearized CME and energy balance equation in equations (69) and (70) may be organized in frequency blocks and rewritten in a compact form as W,AJ,-&AT=

-F,,

(71)

forn=l,...,NF,and nE1(V, AJ, - G, AT) = - E,

(72)

where W, is a N x N block-tridiagonal matrix with elements A,,, B,,, and Cj,.; U,, V,, and G,, are N x N block diagonal matrices, while the others are either column or row vectors of dimension N. Here N is equal to NR x NZ. Eliminating AJ, between equations (71) and (72) to solve for AT, we have “E, (V, W,’ U, - G, )AT= -E+

E V,W,‘F,. n=l

(73)

Essentially we have arranged the system of linearized equations into frequency blocks instead of the usual spatial blocks as proposed by Rybicki*’ and discussed in more details in Leung.‘O In principle, this equation is readily solved for AT (from which AJ,, for all frequencies can be obtained), although the solution involves the inversion of a large array, of dimension NF x (NR x NZ)‘. Iteration then proceeds as outlined in the partial linearization method. The complete linearization method becomes particularly simple in the case of gray (frequency independent) opacities. In this case the discretized equation of radiative equilibrium [equation (4211 reduces to a simple algebraic equation, after integrating over all frequencies E,=Ji--BTj=O,

where BTj is the vector for the frequency-integrated energy balance equation is then AJj - (aBq/aT,)ATj

(74)

Planck function. The linearized form of this = - Ej.

(75)

Similarly the linearized form of the CME [equation (72)] can be written as AjAJi_,+BjAJi+Cj~~,,-DIATj=-Fi,

(76)

where D( = aDj/aTj. Using equation (75) to eliminate the AYs from equation (76), we have Aj(BB~_,/la~_,)A~_,+(-Dol+Bj(aBl;./aTi)ATi +Cj(aBl;+,/a?+,)Aq+,

= -Fi+AjEj_I+BjEj+C,Ej+I.

(77)

Disk geometry solution of the radiation transport equation

We may combine terms and rewrite this equation in block-tridiagonal QjAT,_,+RjAqi-SjAq+,=

Wj,

577

form as (78)

which may be solved for the temperature corrections at each grid point using the standard Gaussian elimination scheme. The AT’s are then used in equation (75) to find the corrections to the mean intensity. Again iteration takes place between solving these linearized equations (to update the temperature and mean intensity) and the ray equations (to update the anisotropy factors). 7.

Selection

PRACTICAL

of spatial, angular, and frequency

CONSIDERATIONS

grids

As with most finite difference methods, the success of this iteration scheme depends on the proper choice of the spatial, angular, and frequency grids. The same considerations, as discussed in Leung”, apply in selecting the spatial and frequency grids. The reader is referred to that paper for details. It is found that the solution to the CRE is highly sensitive to the details of the spatial grid, while the CME solution is more sensitive to the variations in the source function. In selecting the radial grid, care must be taken to assure that the azimuthal angle 4 does not vary too drastically, especially near the core. One must also exercise care to keep the successive variations in grid spacing from being too drastic, which can lead to numerical instability or even non-physical results. Computer

storage and time requirements

Consider a typical model using NR r-points, NZ z-points, NTH e-angles, and NF frequency points. The storage requirement for the computer code comes from three components. First, for the main part of the code, including general arrays, parameters, constants, etc. the storage required is roughly (NR* x NZ + 40 x NR* + 16 x NF + 100) words. Second, temporary storage (scratch arrays) for coefficient matrices to solve the CME takes up about (3 x NR* x NZ) words. Scratch arrays for solution of the CRE require approx. (40 x NR3 x NZ x NTH) words. To run a model using a 22 x 22 grid, with NTH = 5, the total storage allocation required is about 22 Mbytes when using double precision. The actual output file takes up about 10% of this. In terms of computing time, the effort required to solve the block-tridiagonal system of the CME scales as (NR x NZ)* for each frequency. The solution of the CRE involves systems of the same structure. However, the solution for each impact plane involves NTH angles, and each impact plane has (2 x NZ - 1) z-grid points and up to (2 x NR - 1) y-grid points (the CRE must be solved for all four quadrants of an impact plane). Since the solution for each impact plane requires computing time which scales as [(NR x NZ)* x NTH], to solve the CRE for all NR impact planes requires a computing effort which scales as (NR3 x NZ* x NTH). Extension to multiple frequencies increases this requirement by a factor of NF. Hence most of the computer time is taken up by the solution of the CRE. For example, using FORTRAN and double precision, a single (one frequency) solution of the CME on a 30 x 30 grid takes about 1.5 min on our Ridge 32C computer (which has a speed comparable to a VAX-l l/780 and operates at about 3 million instructions per second), while a single solution of the CRE requires about 5 h.

8.

TESTS

AND

APPLICATIONS

There are two internal checks applied to a converged solution of the grain temperature problem. First, we require that the heating and cooling rates for the dust grains be equal locally at each point in the grid, i.e. check the constraint of radiative equilibrium as given by equation (42). Second, we require that the net flux, when integrated over all frequencies, satisfies both locally and globally the flux conservation condition as given by equation (43). For example, for an externally heated cloud without internal energy source, the integrated net flux should be zero. In general, the fist requirement for energy conservation is expected to be satisfied more accurately and in all cases this condition is satisfied to machine round-off accuracy. The degree to which flux conservation is satisfied depends critically on the number of grid points and the details of the grid spacing. For most problems using reasonable grid size, a few percent error is tolerable.

GEORGEF.

578

SPAGNA JR and CHUNMINGLEUNG

Table 2. Test results for the analytic gray solution Model

A

i!

No. of iterations

la lb 2a 2b 3a 3b 4a 4b

2.4048

1.0 0.5 1.0 0.5 1.0 0.5 1.0 0.5

9 9 5 5 4 4 5 4

2.01 1.5 1.0

Maximum error 6.5(-3) 5.8(-3) 5.6(-3) 4.6(-3) 3.7(-3) 2.9(-3) 1.8(-3) 1.4(-3)

A vital test of any numerical code is its ability to reproduce analytic solutions when they exist. Accordingly, we have examined some limiting cases where analytic solutions are available. In particular, we consider the conservative gray case: a dusty medium with uniform density and frequency-independent opacity (gray opacity) under the condition of radiative equilibrium. With these properties (i.e. constant x and xJ = q) and the additional assumption of diffusion approximation (i.e.fil = f, = f++ = 1/3,f;, = 0, and 5 = l), the CME [equation (27)] and the constraint of radiative equilibrium [equation (42)] can be combined to give (d’J/iTz’) + (a2J/Jr2) + (l/r)(aJ/c%)

= 0,

(79)

which we recognize as Laplace’s equation in cylindrical geometry with azimuthal symmetry. This equation has the following known solution J(r, z) = A, cash (&z) Jo@, r) + A, cash (n,z) Z&,r),

(80)

where J,, is the zero-order Bessel function and I0 is the zero-order modified Bessel function. Once the value of J(r, z) is specified at the boundary, the overall solution is uniquely determined. For simplicity, we selected A, = 10 and A2 = 0 so that J(r, z) = 10 cash (AZ)J&r). For our test we fixed the value of J at the boundary, forcing it to match the analytic solution after one iteration. The program was then allowed to iterate from a uniform J until a solution was obtained. The results of several test runs, which are summarized in Table 2, indicate that the computer code reproduces with good accuracy the analytic solution. All models were computed with a radial grid of 20 points and setting R = 1. Axial grids were either 20 points (for 2 = l), or 12 points (for 2 = l/2). While the numerical solution to the conservative gray case serves as a useful check on the accuracy of the computer code and the validity of the numerical method, extension to models with nongray opacity (hence to multi-frequency case) is necessary to study the effects of disk geometry on the dust temperature distribution and the emergent energy spectrum. Following the procedure outlined in Section 6, we have constructed realistic models for a variety of astronomical sources with disk geometry, e.g. protoplanetary nebulae, circumstellar disks and interstellar molecular clouds. As an example to illustrate the effects of disk geometry, we have computed a set of nongray models (using realistic opacities for graphite grains) for interstellar dust clouds which are heated externally by the ambient interstellar radiation field. For these models we have assumed the diffusion approximation in calculating the dust temperature distribution. This is a reasonable assumption for externally heated models since the temperature does not depend strongly on the angular distribution of the radiation field. From the dust temperature distribution we then solve the ray equation at selected wavelengths to determine the emergent intensity and i.r. flux spectrum. The disk models are characterized by an aspect ratio (ratio of radius to half-thickness, R :Z) and an optical depth at 0.55 pm (measured from the cloud surface to the disk center in the midplane). We have asswned a Gaussian density distribution such that the ratio of central to surface density is 100. To illustrate the effects of source geometry, we compare results from models with spherical and disk geometry. Results for the spherical models were obtained using the computer code as documented by Spagna and Leung. ‘I Both models have the same radius (taken to be 1 parsec) and central optical depth (equal to loo), the disk models having 1: 1 aspect ratio. These objects are similar in size and optical depth to quiescent dark clouds known as Bok globules. While the dust temperature distributions in the two cases are very similar as shown in Fig. 4, the excess flux for the disk model depends sensitively on the viewing angle. Here the excess flux is defined as the emergent flux minus the flux from the ambient radiation as measured by an observer when the cloud is not resolved. It represents the energy emitted by the cloud above that

Disk

579

geometrysolution of the radiation transport equation

---

Midplane Axis

Fig. 4. The dust temperature distributions for graphite grains in spherical and disk (with a 1: 1 aspect ratio) models of the same central optical depth (r. = 100) are compared. The dotted line denotes the radial temperature distribution for the spherical model while the dashed and solid lines denote respectively the temperature profiles along the symmetry axis and along the midplane for the disk model.

of the ambient radiation, as received by a distant observer. For clouds which are unresolved, one would naively expect, since the thermal emission is isotropic in the neighborhood of an emitting grain and the emission (in the far i.r.) is optically thin, that the observed flux spectrum should be characteristic only of the dust temperature and independent of viewing angle. However, because of the lack of spherical symmetry in disk geometry and the resulting radiation anisotropy, this is not the case for disk geometry. Figure 5 shows the excess flux spectrum for the disk model viewed

,/ _-..-... I’

_-s

,;-

,*,‘A

:’ :‘I

1

/I

\

.‘--..._ -y %> ‘l\

P; ,,*‘, 3

:f $ i;

7

%.

I’ ,’

‘,\ \\ \ 1’

\

\

‘\ \ \ \

\\

I

- -e

=

8.40 ---8 = 41.50 ______ e = 90.00

Fig. 5. The excess flux spectra for the spherical and disk models of Fig. 4 are compared. The solid line is for the spherically symmetric model. For the disk model, because of the lack of complete symmetry and the resulting radiation anisotropy, the excess flux profile depends on the viewing angle 0 with respect to the symmetry axis. For edge-on viewing (~9= 90”) the excess flux is strongest. Q.S.R.T.

37,6E

GEORGE F. SPAGNA JR and CHUNMINGLEUNG

580

at three different viewing angles (0 is measured with respect to the disk symmetry axis). The brightness of the disk clearly varies with 0, being somewhat brighter seen edge-on than face-on, and with both extremes brighter than when seen at an intermediate angle. The edge-on (0 = 907 spectrum is similar in shape to that of the spherical model, but the flux is greater by a factor of 5. The flux at other viewing angles are weaker and the long wavelength side of the emission peak falls off faster, indicating that we are seeing less of the cooler interior dust. There are essentially two contributions to the variation in observed flux with viewing angle. The first is a purely geometric projection effect. For a disk of radius R, half-thickness Z, and having a uniform surface brightness Z,, the observed flux as seen by an observer at distance D (D % R or Z) depends on 0 as I;disk= (nR2 cos 8 + 4RZ sin 0) IO/D’.

(81)

This implies that for a disk with a 1: 1 aspect ratio, the observed flux is larger seen edge-on than face-on and it reaches a maximum when 8 N 52”. Compared to a sphere with the same uniform surface brightness and radius, the ratio of observed fluxes for the two geometries is just given by the ratio of their projected areas as seen by the distant observer (FdiJFq,,)

= (aR2 cos 8 + 4RZ sin tl)/nR2,

(82)

which indicates that for R = Z, Fdisk2 Fsphere for all viewing angles. The second contribution to the variation comes from radiative transfer effects of disk geometry, which affect the radiation field anisotropy, the dust temperature distributions, and hence the surface brightness. The relative importance of these two contributions depend on such factors as the degree of disk flattening, the cloud optical depth, and the dust density distribution, An important consequence of this variation in excess flux with viewing angle is that it implies a large uncertainty in estimating the radiating dust mass for disk-shaped clouds which are unresolved. In particular, if spherical geometry is incorrectly assumed, the estimated dust mass may easily be off by over an order of magnitude. Other physical effects of radiation transport in disk geometry and further observational implications related to the interpretation of infrared observations of unresolved disk-shaped objects are discussed in Spagna and Leung.” Acknowledgement-This

research was partially supported by NASA under Grant NAGW-577. REFERENCES

1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20.

D. Mihalas, L. H. Auer and B. W. Mihalas, Astrophys. J. 220, 1001 (1975). H. P. Jones and A. Skumanich, Astrophys. J. Suppl. 2, 221 (i980). T. M. Shih, Namer. Heut Transfer 8. 1 (1985). J. Lefevre, J. Y. Daniel and J. Bergeat. >str&. Astrophys. 121, 51 (1983). L. H. Auer and D. Mihalas, Mon. Not. R. a&r. Sot. 149, 65 (1970). D. G. Hummer and G. B. Rybicki, Mon. Not. R. astr. Sot. 152, 1 (1971). D. G. Hummer, C. V. Kunasz and P. B. Kunasz, Comp. Phys. Commun. 6, 38, (1973). D. Mihalas and D. G. Hummer, Astrophys. J. Suppl. 28, 343 (1974). C. M. Leung, Astrophyk. J. 199, 340 (1975). C. M. Leung, JQSRT 16, 559 (1976). G. F. Spagna Jr and C. M. Leung, Camp. Phys. Commun. 28, 337 (1983). M. P. Egan, C. M. Leung and G. F. Spagna Jr, submitted to Camp. Phys. Commun. (1987). D. Mihalas, Stellar Atmospheres, 2nd edition. Freeman, New York (1978). L. H. Auer, JQSRT 11, 573 (1971). G. F. Spagna Jr and C. M. Leung, Astrophys. J. In press (1987). S. Chandrasekhar, Rudiutiue Transfer. Dover, New York (1960). C. D. Levermore,. JQSRT 31, 149-(1984). L. Fox, Numerical Solution of Ordinary and Partial Differential Equations. Addison-Wesley, Reading, Mass. (1962). E. Isaacson and H. B. Keller, Analysis of Numerical Methods. Wiley, New York (1966). G. B. Rybicki, JQSRT 11, 589 (1971).