- Email: [email protected]
On laboratory-frame radiation hydrodynamics
Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97 www.elsevier.com/locate/jqsrt
On laboratory-frame radiation hydrodynamics Dimitri Mihalas ∗ , L.H. Auer Group X-6, MS-D409, Applied Physics Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA Received 28 September 2000
Abstract We derive and discuss the equations of radiation hydrodynamics with all radiation quantities expressed in the laboratory frame. Relativistically exact transformations from the comoving to laboratory frame are used. We obtain simple, fully covariant, expressions for the radiation energy and momentum source=sink terms using lab-frame radiation quantities. The mathematical simplicity of having lab-frame radiation quantities in conservative operators on the left-hand side is maintained. The resulting system can be solved e5ciently using an Eddington tensor formulation for angular dependences, and an iterative frequency-splitting, in the spirit of the multifrequency-gray method. Use of lab-frame radiation quantities can present numerical di5culties in making an accurate connection to the di8usion limit for extremely optically thick media; but it may have no di5culty for neutron-transport problems where very large scattering-thicknesses are generally not encountered (for non-critical systems). Published by Elsevier Science Ltd. Keywords: Lab-frame radiation hydrodynamics; Opacity distribution functions
1. Introduction In this paper, we discuss radiation hydrodynamics when all radiation quantities are evaluated in the laboratory frame, i.e. the
0022-4073/01/$ - see front matter. Published by Elsevier Science Ltd. PII: S 0 0 2 2 - 4 0 7 3 ( 0 1 ) 0 0 0 1 3 - 9
62
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
energy and momentum equations. In Section 4 we combine the radiation energy and momentum equations with nonrelativistic hydrodynamical equations to get a system for nonrelativistic radiating =ow. We o8er alternative expressions for the momentum and energy equations (see also Refs. [2,3]). In Section 5, we quote Thomas’s transformation laws [4] for the radiation
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
63
departures on the micro-scale from the strict isotropy of a Maxwellian distribution usually can be accounted for with transport coe5cients (conductivity, viscosity) calculated from, say, Chapman–Enskog theory. In this limit one can analytically integrate over phase space to derive explicit expressions, depending only on local values of variables and their gradients, for all macroscopic =ow quantities appearing in the hydro equations. Thus for pure hydro (no radiation) the system is at most four-dimensional. (For free molecule :ow, where particle mean-free-paths become comparable or large compared to gradient lengths, it is necessary to include nonlocal e8ects in the distribution function by solving the full Boltzmann equation. We do not discuss such problems this paper.) Therefore for material, the comoving frame is natural because all physical properties can be de
64
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
density because of advection. The mathematical consequence of these physical e8ects is that photons, moving with >xed laboratory-frame angle and frequency, are observed to move on a curved trajectory with varying frequency in the comoving frame; hence partial derivatives with respect to these variables appear in the di8erential operator describing the transport process. (In geometric language, the di8erential operator must account for changes in the basis vectors of the tangent spaces attached to the space-time manifold along the null geodesic of a photon.) The resulting transport equation is seven-dimensional (4 space-time variables, 2 angle cosines of a photon’s direction, and frequency). Explicit expressions for the comoving-frame transport operator may not be available in the literature for a particular geometry or symmetry. The comoving-frame transport equation, correct to O(v=c) for nonrelativistic 1D =ows in spherical symmetry is derived in Ref. [6]. An exact comoving equation for relativistic =ows in 1D spherical symmetry is given in Ref. [7]. The comoving equation, correct to O(v=c) is written for general curvilinear, orthogonal, coordinates in Ref. [8], and generalized to relativistic =ows in Ref. [9]. A general, coordinate-free, relativistic, expression for the comoving equation is given in Ref. [10] for arbitrary geometries. But to obtain explicit expressions for a particular geometry from the results in Refs. [8–10] may require considerable algebra. Explicit expressions are given in Ref. [10] for spherically symmetric, 1D, relativistic =ows, and to O(v=c) in: (1) coordinate-free form, (2) Cartesian coordinates and 1D slab symmetry, (3) 1D spherical coordinates and symmetry, (4) 2D cylindrical coordinates and symmetry, and (5) in conservative coordinate-free form. Monochromatic and total radiation energy and momentum equations, obtained by integrating over angles or both angles and frequency, are given in all of the papers just cited; see particularly [11]. In contrast, in lab-frame transport equation, photons are seen to move on straight-line paths with
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
65
di5culty for a laboratory-frame calculation is that in the di8usion regime, the di8usive =ux, which depends on the local temperature gradient, may be dwarfed (by many orders of magnitude) in the total radiant =ux seen transported by a
66
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
spatial derivatives, and the Lagrangian time-derivative, which tracks the motion of the =uid, namely: D(•) @(•) ≡ + v · ∇(•): Dt @t
(1)
Many writers would call this description “Lagrangian”. Pomraning calls it a “modixed in the :uid, cf. Ref. [6, Sections 9:7 and 9:8]. As shown rigorously in Ref. [14], on an adaptive mesh, Eq. (1) is replaced by d(•) @(•) ≡ + vgrid · ∇(•): dt @t
(2)
Clearly when vgrid ≡ 0 one recovers the Eulerian time derivative, and when vgrid ≡ v=uid one recovers the Lagrangian time derivative. We stress again that these transformations are kinematical, and Galilean. Alternatively, one can write the transport equation and its moments (dynamical equations) for comoving-frame radiation quantities, on aeld, especially in transparent parts of the =ow. In practice one splits inner iterations on (1) the angular distribution of the radiation, from (2) those for the spectra of the radiation energy density and =ux, and both from an outer iteration of the radiation-hydro equations, which yield total energy and momentum balance. This process should be iterated to convergence at a given time-level. The determination of the radiation spectra and angular distribution is hardly “subsidiary” in terms of the computation required. In the total energy and momentum equations one has only 9 variables (; T; v; E; and F) per grid-point; but for the angle-frequency dependence of the radiation
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
67
3. The radiation equations 3.1. The transport equation One can describe the radiationc intensity I (x; t; n; ), a function of position x, time t, the direction n (a unit vector) of a ray, and the frequency of the radiation. The units of I are (erg=cm2 =s1 =Hz1 =sr 1 ), where “sr” denotes solid angle measured in steradians. Let = (x; t; n; ) denote the total opacity per unit pathlength of the material, and = (x; t; n; ), denote the total emissivity per unit volume of the material. In the comoving frame both and are functions (in LTE and for isotropic scattering) of frequency, density, and temperature only. For brevity we will suppress explicit reference to and T because both are functions of x and t). As described in Section 2, for moving material seen in the lab frame and also depend on the direction of propagation n of the radiation, and the =uid velocity v. Note that our de
(3)
Here (I=s) is an absolute (or intrinsic) derivative, which is de
(4)
But Rx = n ds, and the time to cross the element is Rt = ds=c, hence @I 1 @I + ni i = − I: @x c @t
(5)
We use the summation convention on repeated indices throughout this paper. Eq. (5) holds only in Cartesian coordinates. In general, the transport equation is a six- or seven-dimensional system: a four-dimensional spacetime, two angles, and perhaps frequency. A direct solution of such dimensionality is daunting, so one seeks simpli
68
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
3.2. The radiation energy equation De
E = monochromatic radiation energy density in lab frame ≡ c I d!; F i = monochromatic radiation energy =ux in lab frame ≡ I ni d!; ij −1 P = monochromatic radiation pressure in lab frame ≡ c Ini nj d!: Integration of the transport equation over solid angle ! gives the lab-frame monochromatic radiation energy equation @E(x; t; ) @F i (x; t; ) = [(x; t; n; ) − (x; t; n; )I (x; t; n; )] d! (6) + @t @xi and integration over frequency gives the total radiation energy equation @E(x; t) @F i (x; t) + @t @xi ∞ = d d![(x; t; n; ) − (x; t; n; )I (x; t; n; )] ≡ −cG 0 : 0
(7)
Here cG 0 is the (rate of radiation energy absorbed by the material) − (rate of radiation energy emission by the material) per unit volume. All quantities in Eqs. (6) and (7) are measured in the lab frame. 3.3. The radiation momentum equation Similarly, integration of the transport equation against ni d! over solid angle gives the lab-frame monochromatic radiation momentum equation 1 @F i (x; t; ) @P ij (x; t; ) + c2 @t @xj 1 = [(x; t; n; ) − (x; t; n; )I (x; t; n; )]ni d! c
(8)
and integrating over frequency, gives the total radiation momentum equation 1 @F i (x; t) @P ij (x; t) + c2 @t @xj ∞ 1 = d d![(x; t; n; ) − (x; t; n; )I (x; t; n; )]ni ≡ −G i : c 0
(9)
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
69
Here G i is the rate of (deposition of the ith component of the radiation momentum into the material) − (emission of the ith component of the radiation momentum by the material) per unit volume. 3.4. Conservation laws The di8erential operators in the lab-frame radiation and momentum equations are exact, and are in conservation form. Their apparent simplicity is, however, compromised by the di5culty of evaluating accurately the integrals on the right-hand sides in the lab frame. As described in Section 2 the principal di5culty is that for moving material, radiation moving with a
(10)
(11)
Because of symmetry f has 1, 3, and 6 nonredundant components in 1D, 2D, and 3D, respectively. The elements of f depend upon (x; t; ). We suppress explicit mention of these variables to emphasize that, at a given stage of the calculation, f is a collection of angular factors computed from an angle-frequency-dependent transport equation with given source=sink terms, as discussed in Section 9. In practice this splitting works because the elements of f are ratios, which re=ect primarily the presence of open boundaries. Beyond one optical depth from a boundary, values of fij can be accurate even when current estimates of the source=sink terms used in the transport equation are not yet well determined by the outer iteration procedure which gives material properties, and E and F, that satisfy the radiation hydro equations at the advanced time. But we emphasize that the solution of the radiation-dynamical equations can be no better than the accuracy of the tensor elements used to close the system. Thus if one merely assumes
70
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
the lab-frame Eddington tensor is isotropic in a moving =uid, then the radiation stress is forced to act like a hydrostatic pressure (like the gas pressure). This assumption would be valid in the comoving frame if "p =‘1, (where "p is a photon mean free path and ‘ is a characteristic gradient length in the =uid). But one needs to Lorentz transform to the lab frame. And if there is an open boundary, or the opacity of the material is small at some frequencies, this approximation may not be valid even in the comoving frame. When f in the lab frame has nonzero o8-diagonal elements, radiation stresses can induce shear in the =ow. A physically correct solution requires Eddington-tensor elements that account fully for the problem’s geometry, the frequency dependence of photon mean-free-paths, and the e8ects of the =uid’s velocity; we address some aspects of this problem in Section 9. 4. Lab-frame radiation hydrodynamics equations In nonrelativistic laboratory =ows, where $ ≡ v=c1, the
(12)
(13)
(14)
(c) the material energy equation @(e0 ) @vi @ + i (e0 vi + q0i ) + (p0 ij − '0ij ) j = *0 + cG 0 − vi G i @t @x @x
(15)
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
or the radiating =uid energy equation @ @ @vi (e0 + E) + i (e0 vi + q0i + F i ) + (p0 ij − '0ij ) j = *0 − vi G i @t @x @x or the material total energy equation @ @ [( 12 v2 + e0 )] + i {[(( 12 v2 + e0 ) + p0 ]vi − '0ij vj + q0i } = vi fi + *0 + cG 0 ; @t @x or the radiating =uid total energy equation @ @ [( 12 v2 + e0 ) + E] + i {[(( 12 v2 + e0 ) + p0 ]vi − '0ij vj + q0i + F i } = vi fi + *0 ; @t @x (d) the radiation energy equation @E @F i + i = −cG 0 @t @x and (e) the radiation momentum equation 1 @F i @P ij + j = −G i : c2 @t @x
71
(16)
(17)
(18)
(19)
(20)
In Eqs. (12) – (18), p0 is the gas pressure, e0 the =uid’s internal energy, '0ij its viscous stress tensor, q0i the thermal heat conduction vector, and *0 the rate of nuclear energy release per unit mass. By solving these equations one gets ; v; T0 ; E, and F. The di8erential operators are all in conservation form, which makes them attractive computationally. The two basic problems still to be faced are computation of elements of the Eddington tensor f, and lab-frame values of the four-force G + . We address these issues in Sections 9 and 6, respectively. 5. Transformation laws 5.1. Transformation of the photon four-momentum Denote comoving-frame quantities with a subscript “0” and leave lab-frame quantities unadorned. The connection between and 0 , and n and n0 , is given by Lorentz transformation of the photon four-momentum (h=c)(1; n). Transforming from the =uid frame to the lab frame one
n =
0
ni0
%$j nj0 + 1+ %+1
(21a)
i
%$ :
(21b)
The inverse transformation is 0 = %(1 − $i ni )
(22a)
72
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
and
%$j nj = n − 1− %$i : %+1 Here $i ≡ vi =c; $2 ≡ $i $i , and % ≡ 1= 1 − $2 .
ni0
0
i
(22b)
5.2. L.H. Thomas’s transformation laws In a seminal paper, Thomas [4] showed by simple, but genuinely covariant, Gedanken experiments, I (n; ) = (3 =30 )I0 (n0 ; 0 ) ≡%
−3
i −3
(1 − $i n )
(n; ) = (2 =20 )0 (0 ) ≡
I0
n − %(1 − %$i ni =(% + 1))R i ; %(1 − $i n ) ; %(1 − $i ni )
0 [%(1 − $i ni )] %2 (1 − $i ni )2
(23) (24)
and (n; ) = (0 =)0 (0 ) ≡ %(1 − $i ni )0 [%(1 − $i ni )]:
(25)
Notice that isotropy of the material properties 0 and 0 in the comoving frame is assumed. The second form of Eqs. (23) – (25) displays explicitly the angular dependences that result from transforming from the comoving frame to the lab frame. Eqs. (23) – (25) say that I=3 ; =2 , and , are Lorentz invariants. In addition, for any particle, d 3 p= e; ˜ where d 3 p = dpx dpy dpz = p2 dp d! is the three-space di8erential of the particle’s momentum, and e˜ is the particle’s energy, is also a Lorentz invariant (see [3, Eq. (43.17)]; [5, Eq. (A.65)]). Hence for photons d d! = 0 d0 d!0
(26)
d! = (20 =2 ) d!0 :
(27)
and We will use these results repeatedly. 5.3. Covariance of radiation four-force and stress-energy tensor Although I; E; F; P; , and , have all been taken to be in the lab frame, it should be noted that G + ≡ (G 0 ; G) ≡ (G 0 ; G 1 ; G 2 ; G 3 )
(28)
is a four-vector, which can be transformed from frame to frame by Lorentz transformation. One can verify this statement by using the relations given in Section 5.2. One sees that each
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
73
component of G + as given by Eqs. (7) and (9) is (c=h) times a component of the photon four-momentum M + ≡ (M 0 ; M) ≡ (h=c)(1; n)
(29)
times either the Lorentz invariant emissivity =2 , or the Lorentz invariant opacity times the invariant speci
%$j G j i =G −% G − $: %+1 The inverse transformation is
G0i
i
0
G 0 = %(G00 + $i G0i ) and
j G %$ j 0 G i = G0i + % G00 + $i : %+1
(30a) (30b) (31a)
(31b)
Likewise, the left-hand sides of these equations are c times the four-divergence of the 4 × 4 radiation stress-energy tensor R, where ∞ R+$ ≡ c−1 d d! I (n; )n+ n$ : (32) 0
Here one uses the convention that n0 = 1. In tensor form E c−1 F† R≡ : c−1 F P
(33)
R is a symmetric four-tensor, which can be transformed into the comoving frame by Lorentz transformation. One can verify this statement by noting each element of R is a component of the
outer product of c=h times the photon four momentum with itself, times the invariant speci
(34a)
F0i =c = %{(F i =c) − %$i E − $j P ij + $i [%=(% + 1)][(2% + 1)($j F j =c) − %($j $k P jk )]}
(34b)
and P0ij = P ij − %[($i F j + $j F i )=c] + %2 $i $j E + [%2 =(% + 1)][$j $k P ik + $i $k P kj − 2%$i $j ($k F k =c)] + [%2 =(% + 1)]2 $i $j ($k $l P kl ):
(34c)
74
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
The same results are found by using Eqs. (21a), (21b), and (26) in the
(35a)
F i =c = %{(F0i =c) + %$i E0 + $j P0ij + $i [%=(% + 1)] j
jk
×[(2% + 1)($j F0 =c) + %($j $k P0 )]}
(35b)
and P ij = P0ij + %[($i F0j + $j F0i )=c] + %2 $i $j E0 + [%2 =(% + 1)][$j $k P0ik + $i $k P0kj + 2%$i $j ($k F0k =c)] + [%2 =(% + 1)]2 $i $j ($k $l P0kl ):
(35c)
Covariant (both in coordinates and relativistically) radiation dynamical equations are @R+$ = −G + : (36) @x$ The total radiation energy and momentum equations are direct analogues, for radiation, of the energy and momentum equations for the material part of the =uid. As mentioned before, they are exact. But to apply them one still needs Eddington tensor elements, and usable expressions for the four-force. 6. Lab-frame radiation source= sink terms 6.1. Comoving-frame source=sink terms Consider
(37)
Likewise for the emissivity, one discriminates between (the isotropic) radiation emitted thermally, and radiation scattered from the incident radiation
(38)
For brevity, explicit mention of the dependence of 0 ; 10 ; '0 , and 0 on x and t has been suppressed.
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
75
Assuming LTE, one invokes Kirchho8’s law to write the thermal emissivity as t0 (0 ) = 10 (0 )B(0 ; T0 );
(39)
where B(0 ; T0 ) is the Planck function, which is isotropic. The scattering source is more complicated. Let the redistribution function R(n0 ; 0 ; n0 ; 0 ) be the probability that an incoming photon at (n0 ; 0 ) in the comoving frame emerges in the comoving frame at (n0 ; 0 ). R is normalized such that ∞ d!0 d!0 ∞ d0 d0 (40) R(n0 ; 0 ; n0 ; 0 ) = 1: 43 0 43 0 Then the scattering emissivity is ∞ d!0 s0 (n0 ; 0 ) = '0 d0 R(n0 ; 0 ; n0 ; 0 )I0 (n0 ; 0 ): 43 0
(41)
Here '0 is the total scattering cross-section. Redistribution functions can be calculated analytically for various models of redistribution in angle and frequency by the scattering process in an atom’s frame. Scattering shows another aspect of the di5culty in solving the transfer equation. Not only is it a 7D partial di8erential equation in spacetime, angles, and frequencies, but it is an integro-partial di8erential equation in both the input and output angle-frequency spaces. One must therefore seek simpli
(42)
(43)
Thomson scattering has a mild dipole phase factor, which, for simplicity, we take as isotropic, g ≡ 1, and coherent, R(0 ; 0 ) ≡ (0 − 0 ). Hence d!0 c'0 (0 ) c'0 (0 ) s 0 (0 ) = d0 (0 − 0 ) (44) I0 (n0 ; 0 ) = E0 (0 ): 43 43 43 In this simple case s0 (0 ) is also isotropic. The total comoving-frame emissivity is 0 (0 ) = t0 (0 ) + s0 (0 ) = 10 (0 )B(0 ; T0 ) + '0 (0 )[cE0 (0 )=43]:
(45)
We adopt this simpli
76
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
6.2. Comoving-frame energy and momentum input terms Calculating the net radiant energy input and momentum input into the material, measured in the comoving frame, one
and G0i
≡c
−1
∞
0
d0
d!0 [0 (x; t; 0 )I0 (x; t; n0 ; 0 ) − 0 (x; t; 0 )]ni0
and taking advantage of the isotropy of 0 and 0 in this frame, one has ∞ 0 cG0 = 43 10 (x; t; 0 )[(c=43)E0 (x; t; 0 ) − B(x; t; 0 ; T0 )] d0 0
and G0i = c−1
0
∞
(x; t; 0 )F0i (x; t; 0 ) d0 :
(46b)
(47a)
(47b)
Both G00 and G0i are simpler than G 0 and G i because (a) they contain only moments of the radiation
(49)
Here E0 (0 ) ≡ E0 (0 )=E0 is the comoving-frame radiation energy spectrum, which is normalized: 0∞ E0 (0 ) d0 = 1. Then cG00 = c10E E0 − 4310P B(T0 ) = c(10E E0 − 10P aR T04 ):
(50)
Here aR is Stefan’s radiation constant. Both 10E and 10P are scalars. Whereas 10P is a function of density and temperature only, evaluation of 10E requires knowledge of the comoving-frame
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
77
radiation energy spectrum E0 (0 ). If one can actually evaluate this spectrum, then no approximation is made by using 10E and 10P in G00 . (i) Similarly, for G i0 , de
where F(i) (0 ) ≡ F0i (0 )=F0i is the comoving-frame radiation :ux spectrum in the ith direction, which also is normalized: 0∞ F(i) 0 (0 ) d0 = 1. Then (i) i G0i = 0F F0 =c:
(52)
(i) is not a function of (0 ; T0 ) only; nor is it, in general, a scalar, because the =ux Like 10E ; 0F (i) spectrum can be di8erent along each axis. One should regard 0F as a 3 × 3 diagonal matrix. (i) i Note that F0 along a given axis may be zero, so that 0F in that direction is unde
6.4. Lab-frame energy and momentum input terms In the laboratory frame, the net radiant energy and momentum input to the material is given by Eqs. (7) and (9). It is not obvious that one can express G 0 and G i solely in terms of radiation moments in that frame. And, as mentioned above, accurate evaluation of the double integral in Eq. (9) can be di5cult in opaque material because of the near cancellation of the intensity along rays in opposite directions. One can deal with these problems as follows:
and i
G =
(i) i (0F F0 =c)
+ % 10E E0 −
10P aR T04
% (j) j + $j (0F F0 =c) $i : %+1
(53a)
(53b)
The terms containing the temperature result from thermal emission. We thus have lab-frame energy and momentum input terms expressed in terms of comovingframe moments, and frequency integrals evaluated in the comoving frame. These expressions are exact to all orders of (v=c); the only simpli
78
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
Eqs. (53a) and (53b) are still not suitable for use in Eqs. (7) and (9) because lab-frame radiation variables would appear on the left-hand sides, and comoving-frame quantities on the right. However, because only frequency-integrated moments appear in Eqs. (53a) and (53b) one can remove the discordance by expressing E0 and F0 in the lab frame by using Eqs. (34a) – (34c). (Note that T0 and the opacities do not require re-transformation because they are world scalars (see Refs. [17, p. 29]; [18]). Carrying out these substitutions, and grouping terms in ascending orders of $, one
− [%4 =(% + 1)]$i 0F $i ($j $k fjk )E:
(54a)
(i) can be replaced Or, assuming that the =ux spectrum is the same in all directions, so that 0F by a scalar 0F ,
G 0 = %[%2 10E + (1 − %2 )0F ]E − %10P aR T04 − %($i F i =c)[0F − 2%2 (0F − 10E )] − %3 (0F − 10E )($i $j fij )E:
(54b)
Similarly, (i) i (i) G i = %(0F F =c) − {[%0F ($j fij + %$i ) − %2 10E $i ]E − %10P aR T04 $i }
+
%2 (i) (j) {[(2% + 1)0F − %10E ]($j F j =c) + %($j 0F F j =c)}$i %+1
+
%3 (i) (j) {[(% + 1)10E − 0F ]($j $k fjk ) − $j 0F (%$j + $k fjk )}$i E %+1
+
%4 ( j) j ($j 0F $ )[(2% + 1)($k F k =c) − %$k $l fkl E]$i : (% + 1)2
(54c)
Or, again assuming a scalar 0F can be used, G i = %0F (F i =c) − %10P aR T04 $i − [%3 (0F − 10E )$i + %0F $j fij ]E + %3 (0F − 10E )[2($j F j =c) − ($j $k fjk )E]$i :
(54d)
Eqs. (54a) – (54d) are covariant, both geometrically and relativistically, and provide workable expressions for (G 0 ; G) correct to all orders of v=c. All the radiation variables, including Eddington tensor elements, are measured in the lab frame. In contrast, the material variables (i) 0 ; T0 ; 10P ; 10E ; and 0F are all measured in the comoving frame.
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
79
6.5. Check of formulae by using thermal radiation One can check the algebra required to reduce Eqs. (53a) and (53b) to Eqs. (54a) – (54d), by considering an in
(55a)
F i =c = − 43 %2 $i aR cT04
(55b)
P ij = 13 (ij + 4%2 $i $j )aR T04 :
(55c)
and
Eq. (55b) shows that the radiation =ux has a dipole component in the lab frame, proportional to the negative of the lab’s velocity relative to the =uid frame. The o8-diagonal components of the lab-frame radiation-pressure tensor are O(v2 =c2 ). Even if radiation is isotropic in the =uid frame, it is not in the lab frame. In fact, the observed dipole component of the cosmic microwave background has allowed determination of the motion of our galaxy relative to the “rest frame of the Universe”. From Eqs. (55a) and (55c) the lab-frame Eddington tensor elements for thermal radiation are fij ≡ P ij =E = (ij + 4%2 $i $j )=(4%2 − 1);
(56a)
whence one
(56b)
$i $j fij = (4%2 − 3)$2 =(4%2 − 1):
(56c)
and
As noted above, the o8-diagonal terms of the Eddington tensor are only O(v2 =c2 ). Nevertheless, at high temperatures, radiation pressure dominates gas pressure by orders of magnitude, hence even a small anisotropy in the radiation pressure may have important consequences for the =ow dynamics. If one assumes strict thermal equilibrium in the rest frame of the =uid, one knows that G00 ≡ 0, and G0 ≡ 0. Then the principle of equivalence between inertial frames states that an observer in the lab frame must obtain the same result: i.e. the =uid is observed to be in thermal equilibrium and is not accelerated by radiation forces. Hence one must
80
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
6.6. Lab-frame radiation dynamical equations Use of Eqs. (54a) – (54d) in Eq. (36) gives a
0 @R+$ (1=c)(@E=@t + @F j =@xj ) G E = =− =A · + B: 1=c2 @F i =@t + @(fij E)=@xj Gi F j =c @x$
(57)
Here + runs from 0 to 3, i and j run from 1 to 3, A is a 4 × 4 matrix, and B is a vector of length 4. This system has an exact lab-frame (four-) divergence on the left-hand side, and a linear combination of the lab-frame energy density, =ux components, and Eddington tensor, correct, in principle, to all orders of (v=c) on the right. The actual accuracy of the system will depend, of course, on the accuracy of the spectral pro
7. Evaluation of mean opacities and radiation spectra Consider now the evaluation of mean opacities. The opacity spectrum of elements other than hydrogenic or helium-like ions is complicated, and a million frequency sample-points may be needed to characterize it adequately for mixtures of elements. Direct calculation of the integrals in Eqs. (49) and (51) is then prohibitive. A number of approximations, ranging from over-simpli
(58)
G i = 1[(F i =c) − (aR T04 + 13 E)$i ];
(59)
and
which agree with the third and fourth un-numbered equations following Eq. (36) in Ref. [19]. Other approximate source terms based on less restrictive assumptions can also be deduced from Eqs. (54a) – (54d). But no real material is gray, and even if this assumption is useful in studying the mathematical behavior of the equations it is an inadequate representation of reality.
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
81
We can still gain some important insight from Eq. (59); using it in the radiation momentum equation we have, to O(v=c), "p @F i @P ij (60) + c"p j = −F i + 43 Evi : c @t @x There are three characteristic timescales to be considered: the radiation-=ow time tr ∼ "p =c, the =uid-=ow time tf ∼ ‘=v, and the di8usion time td ∼ (‘="p )2 ("p =c) = ‘2 ="p c. The ratio of the time derivative of the =ux to the =ux itself is O("p =‘); O("p v=‘c), and O("p2 =‘2 ) on the radiation-=ow, =uid-=ow, and di8usion timescales, respectively. For di8usion theory to apply, "p =‘1, hence the time derivative can always be neglected in a very opaque material. Therefore, in the di8usion regime F i = −c"p
@P ij 4 i + 3 Ev = F0i + 43 Evi ; @xj
(61)
which is Eq. (35b) to
82
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
Fig. 1. Upper panel: Planck function, normalized to its maximum value, at T0 = 0:02 KeV, versus frequency in KeV. Lower panel: Spectrum of iron at T0 = 0:02 KeV and = 10−4 gm=cm3 , versus frequency in KeV. Vertical dotted lines illustrate division of frequency into
of the Planck function B(; T0 ). As shown in Fig. 1, this variation is easily resolved. From tables of the Planck function (e.g. Ref. [20, pp. 106 –107]) one sees that a range of fewer than three decades in u ≡ h=kT0 , namely 0:2 6 u 6 25, is su5cient to cover the full variation of B(; T0 ). One then divides the spectrum into bands, each bounded by frequencies (b ; b+1 ); b = 1; : : : ; B. The band sizes :b ≡ b+1 − b need not be equal, but, once chosen, must be the same for all values (0 ; T0 ) in a table. Each band should be narrow enough that the Planck function varies only mildly within it, so that b+1 Bb (T0 ) ≡ B(0 ; T0 ) d=:b (62) b
is well de
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
83
Fig. 2. Planck function normalized to its maximum value (upper panel), and spectrum of iron (lower panel) at T0 = 0:02 KeV and = 10−4 gm=cm3 , versus frequency on the range 0:04 . . 0:075 KeV.
The second source of variation is from the variation of the material’s opacity with energy. As seen in Fig. 1, which is for Iron at T0 = 0:02 KeV and = 10−4 g=cm, this variation is extremely complicated and much more di5cult to treat than the smooth variation of the Planck function. A subset of these data on the range 0:04 . . 0:75 is shown in Fig. 2. Although the Planck function varies by only 10% over this range, the opacity varies by orders of magnitude over the width of each spectral line. To resolve this variation fully we would need O(:b =:D ) frequency points in a band, where :D is the Doppler width of a typical line. At this resolution we would require hundreds of thousands of frequency points per band. The treatment of a frequency spectrum in this much detail in time-dependent problems is beyond the capability of any existing computer. Yet, the frequency variation of the opacity is a de
84
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
Currently, the technique proven to be the best method for treating complex spectra in astrophysical work [21,22] is the opacity distribution function. Once the bands have been de
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
85
Fig. 3. Opacity distribution function with nine ODF pickets for the Iron spectrum shown in Fig. 2.
7.3. Full transport Consider
(64a)
or, if the scattering terms are canceled out analytically, i @Epb @Fpb = 1pb [43Bb (T0 ) − cEpb ] + @t @xi
(64b)
and from Eq. (8) the radiation momentum equation is ij i pb i 1 @Fpb @(fpb Epb ) =− F ; + j 2 @x c pb c @t
(65)
86
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
(p = 1; : : : ; Pb ); (b = 1; : : : ; B). If we know the Eddington tensor elements, from a separate i . Suppose now the medium ray-by-ray calculation, these equations can be solved for Epb and Fpb is time-independent. At a point far from a boundary, and where ("p =‘1) in all pickets, Epb → Bb (T0 ), and strict energy balance is guaranteed in Eq. (64b). Fortunately, the ODF procedure is strongly convergent with increasing number of ODF pickets. i for all p and b, one can compute Solving Eqs. (64b) and (65) to determine Epb and Fpb approximate energy and =ux means needed in the integrated moment equations: B Pb b=1 Rb p=1 wpb 1pb Epb (66) 10E = B Pb b=1 Rb p=1 wpb Epb and (i) 0F
B
=
Pb i b=1 Rb p=1 wpb pb Fpb : B Pb i b=1 Rb p=1 wpb Fpb
(67)
Having 10E and 0F one can evaluate the source terms in Eqs. (54a) and (54c), and solve system (12) – (20). (i) These means are vastly superior to use of 10P for 10E and 0R for 0F , or use of two-temperature means, because they account for: (1) the huge di8erences in opacity between lines and continua, (2) decoupling of the local monochromatic energy density from the Planck function, (3) anisotropy of the radiation
8. Velocity gradient e/ects Using an ODF representation of opacity one can obtain accurate comoving-frame energy and =ux spectra, and compute the mean opacities in equations (54). But the e8ects of Doppler shifts, resulting from velocity gradients, have thus far been ignored. To get a rough estimate, one might compare the Doppler shift over a gradient length ‘, or a photon mean-free-path "p , namely ∼ (=c)(dv=ds) min(‘; "p ), with :sb , the interval between sample points in the opacity spectrum. If :sb , true in very opaque material, one might argue that Doppler shifts are negligible. But caution is needed even then, because they produce red- or blue-shifting of the radiation’s spectrum under adiabatic expansion or compression [23,3, p. 474]. Near boundaries (where "p → ∞), or where the velocity gradient is large, Doppler shifts can be quite large. The monochromatic comoving-frame radiation energy equation is
@F i (0 ) D E0 (0 ) vi @F i (0 ) + 0 i = 10 (0 )[43B(0 ; T0 ) − cE0 (0 )] − 2 0 Dt @x c @t ij @[f0 E0 (0 )] @vi : (68) + 0 @0 @xj
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
87
The dominant terms of Eq. (68) are written on the
@[f0ij E0 (0 )] vj @[f0ij E0 (0 )] D F0i (0 ) 0 (0 ) i 1 @vi j + = − ( ) − − F (0 ) F 0 0 c2 Dt @xj c c2 @t c2 @xj 0 ijk 1 @[0 Q0 (0 )] @vj : (69) + 2 c @0 @xk Eqs. (68) and (69) are accurate only to O(v=c), which is su5cient to determine spectra. In Eq. (69), Q0ijk (0 ) ≡ I0 (n0 ; 0 ) ni0 nj0 nk0 d!0 : (70) One can write an approximate closure relation for this three-index tensor by assuming a P1 approximation for the radiation
(71)
Q0ijk (0 ) ≈ ( 15 + 25 ij )jk F0i (0 ):
(72)
Using Eq. (72) in Eq. (69), one
@[f0ij E0 (0 )] vj @[f0ij E0 (0 )] D F0i (0 ) 0 (0 ) i 1 @vi j + = − ( ) − − F (0 ) F 0 0 c2 Dt c c2 c2 @xj 0 @xj @t 1 @[0 F0i (0 )] @v j @v(i) + 2 + 2 : (73) 5c @0 @xj @x(i) In Eq. (73), the parentheses on the index i in the last term mean no summation on i. The dominant terms are the divergence of P0ij on the left-hand side and the =ux-attenuation term on the right-hand side. Compared to the latter, the time derivative of F0i on the left-hand side is O("p v=‘c) on a =uid-=ow timescale, and O("p =‘) on a radiation-=ow timescale, hence it is important only in regions where the medium is transparent. On a radiation-=ow timescale, the term containing the time derivative of P0ij on the right-hand side is O(v=c) compared to the divergence of P0ij on the left-hand side, and is O(v2 =c2 ) on a =uid-=ow timescale. The term containing the =ux times the velocity-gradient =ux on the right-hand side is always O("p v=‘c) compared to the =ux attenuation term. The frequency derivative of F0i is O("p v=‘c) compared to the =ux-attenuation term, thus the frequency redistribution of the =ux can also be treated iteratively.
88
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
In writing
ij ij @(fpb Epb ) pb i vj @(fpb Epb ) D Fpb 1 @vi j + = − − − F F c2 Dt @xj c pb c2 @t c2 @xj pb b i i 1 ( p vpp Xp b Fp b − Xpb Fpb ) @v j @v(i) + 2 + 2 : (75) 5c Rs @xj @x(i) Here it is assumed that time and frequency di8erencing is done in the usual way. As pointed out above, the frequency derivatives on the right-hand side of Eqs. (74) and (75) are small compared to terms on the left-hand side. Thus these equations need not be treated as coupled in frequency, but can be solved for each frequency, in parallel, using previous estimates of Epb i i , and Fpb in the frequency di8erence terms. This process yields new values for Epb and Fpb which can then be used in the frequency di8erences, and the procedure iterated to convergence. i , the Our expectation is that the process will converge rapidly. Given values for Epb and Fpb mean opacities can be computed from Eqs. (66) and (67). We mention in passing that the structure of the r and v matrices can provide guidance about the quality of the ODF discretization. For example, if the original spectrum is over-resolved, then most opacity sample-points binned into a given ODF picket will likely couple to other points binned into the same picket, and both r and v will be nearly (i.e. too strongly) diagonal. Ideally these matrices should be nearly bidiagonal. On the other hand, if the sampling width :sb in the opacity spectrum is too large, then sampling of line pro
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
89
to construct correlation matrices connecting the ODF picket representations on one side of the interface with those on the other. 9. The lab-frame Eddington tensor for a moving 1uid 9.1. Analytical lab-frame Eddington tensors for simple cases It is possible to calculate the lab-frame Eddington tensor analytically for simple problems. The results may give insight into the size of velocity e8ects, and provide benchmarks for testing numerical codes. 9.1.1. The di
(76a)
F i =c = %{− 43 (@ ln T0 =@xi )= + 43 %$i − [( 43 %$i =)(2% + 1)=(% + 1)] ×(@ ln T0 =@xj )$j }aR T04
(76b)
and P ij = ( 13 ij + 43 %2 $i $j − ( 43 %=){$i (@ ln T0 =@xj ) + $j (@ ln T0 =@xi ) + [2%2 $i $j =(% + 1)][$k (@ ln T0 =@xk )]})aR T04
(76c)
hence ij
f =
(ij + 4%2 $i $j ) −
4% i 2%2 $i $j $k j j i k [$ (@ ln T0 =@x ) + $ (@ ln T0 =@x ) + (%+1) (@ ln T0 =@x )] : (4%2 − 1) − (8%2 =)(@ ln T0 =@xi )$i
(77)
The dominant o8-diagonal terms of the Eddington tensor are only O(v2 =c2 ) or O("p v=‘c) (from the di8usive =ux) relative to the diagonal, hence are generally negligible to O(v=c). The dominant di8erence between the lab-frame =ux F and the comoving frame =ux F0 is 43 v times the radiation energy density; the latter (advective) term can dominate F0 by orders of magnitude in optically thick material, hence must always be taken into account for a moving =uid. 9.1.2. 1D semi-in>nite grey medium Another problem for which an exact solution is known is a 1D, semi-in
(78a)
90
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
and F0 ≡ F0z = constant;
F0x = F0y ≡ 0:
(78b)
The radiation pressure tensor is diagonal, but anisotropic, with cPzz 0 = F0 [7 + q(∞)];
(78c)
yy 3 1 cPxx 0 = cP0 = F0 [7 + 2 q(7) − 2 q(∞)];
(78d)
xy
yz
P0 = Pxz 0 = P0 = 0:
(78e)
Here 7 is the optical depth from the boundary, q(7) is the Hopf function, and F0 is the constant =ux in the outgoing direction. An exact expression for the Hopf function was found√by Mark [24] by Laplace transform techniques. It is a monotonic function ranging from q(0)=1= 3=0:577350, to q(∞) = 0:710446, which can be represented accurately by the variational
(79)
where En (x) denotes the exponential integral of order n. Thus the exact value of the Eddington tensor in the comoving frame of a semi-in
(80a)
f0xx (7) = f0yy (7) = 13 [7 + 32 q(7) − 12 q(∞)]=[7 + q(7)]
(80b)
with all o8-diagonal entries zero. Clearly lim fii 7→∞ 0
= 13 ;
(80c)
showing that the radiation
(81a)
F z = %F0 {1 + 3%$z [7 + q(7)] + $z [7 + q(∞)] + [%$z2 (2% + 1)=(% + 1)] + [%2 $z =(% + 1)]{$2 [7 + q(∞)] + 32 ($x2 + $y2 )[q(7) − q(∞)]}}
(81b)
F x = %F0 {3%$x [7 + q(7)] + $x [7 + 32 q(7) − q(∞)] + %$x $z (2% + 1)=(% + 1) + [%2 $x =(% + 1)]{$2 [7 + q(∞)] + 32 ($x2 + $y2 )[q(7) − q(∞)]}}
(81c)
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
with an identical expression for F y with x and y interchanged, and P zz = F0 7 + q(∞) + 2%$z {1 + [%2 $z2 =(% + 1)]} + 4%2 $z2 7 %2 ($x2 + $y2 ) 3%2 ($x2 + $y2 ) 2 2 2 2 + 3% $z 1 + q(7) + % $z 1 − q(∞) : 2(% + 1)2 2(% + 1)2
91
(81d)
It is not di5cult to calculate other components of P, but these are not useful for the problem being considered, a 1D medium with constant velocity, because in this case (@=@x) and (@=@y) of all quantities are identically zero. Thus a velocity
I () d, c−1 F = c−1 I () d n (n − n) d! = (nx i + ny j + nk k)c−1 I
and −1
(n ⊗ n )(n − n) d! n2x ii nx ny ij nx nz ik = nx ny ji n2y jj ny nz jk c−1 I: nx nz ki nz ny kj n2z kk
P=c
(82b)
I () d
(82c)
Then in the =uid, moving with velocity v in the lab frame, one
(83a)
F0i =c = %(1 − $i ni )(ni − %$i {1 − [%$i ni =(% + 1)]})c−1 I
(83b)
P0ij = (ni − %$i {1 − [%$i ni =(% + 1)]})(nj − %$j {1 − [%$j nj =(% + 1)]})c−1 I:
(83c)
and
92
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
Or, recalling Eqs. (22a) and (22b), one can write, equivalently, E0 = (0 =)2 (I=c); F0i =c = (0 =)2 (I=c)ni0 , and P0ij = (0 =)2 ni0 nj0 (I=c). These expressions are valid only for a collimated beam (i.e. a delta-function in angle in the lab frame). The Eddington tensor in the moving frame has O(v=c) terms which account for the aberration of the beam’s direction as seen by the moving observer. In summary, for the problems considered above, velocity-induced e8ects in the radiation pressure tensor are only O(v2 =c2 ) or O(v"p =c‘) at depth. They are at most O(v=c) near open boundaries, but even then below the angular resolution of typical methods of solving the transport equation. For such problems velocity e8ects are thus negligible in the computation of the Eddington tensor. O(1) departures of the radiation pressure tensor from isotropy re=ect transport e8ects caused by optical geometry, such as boundaries (internal or external), voids in the material, and=or steep gradients in temperature or the radiation
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
93
For thermal emission, t0 (0 ) = 10 (2h30 =c2 )[exp(h0 =kT0 ) − 1]−1 , and from Eqs. (24) and (22a) one has 3 4 4 ∞ 3 3 ∞ 2k T0 10 x dx aR c t (n; ) d = 10 T 4 = 0 c2 h3 ex − 1 0 43 0 0 a c 1 R = 3 10 T04 ≡ t (n): (86) % (1 − $i ni )3 43 For coherent isotropic scattering, Eqs. (44) and (22a) imply 2 c' s (0 ) E0 (0 ) 0 s (n; ) = s0 (0 ) = 2 0 = : i 2 2 0 % (1 − $i n ) 43 % (1 − $i ni )2 Integrating over frequency, and using Eq. (34a), ∞ c' E0 0 s (n; ) d = 3 (1 − $ ni )3 43 % i 0 c' [(1 + $ $ fij )E − (2$ F i =c)] i j i 0 = ≡ s (n): 43 %(1 − $i ni )3
(87)
(88)
Collecting these results in Eq. (84) one
(89)
To check the correctness of the analysis just presented, integrate the absorbtion=emission terms in Eqs. (85), (86), and (88) by taking their zeroth and
d! ni (1 − $i nj ) I (n) = 0 %[(F i =c) − $j P ij ]:
(90b)
For thermal emission, the inner product $i ni appears inside the integral. But because t (n) is a scalar, and $i ni is invariant under a rotation of the coordinate axes, one can choose the polar axis to lie along the velocity vector R. Taking ? to be the polar angle between n and R, de
94
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
To take the
d! ni s (n) = %3 '0 [(1 + $i $j fij )E − (2$i F i =c)]$i :
(92b)
Assembling Eqs. (90a), (91a), and (92a), one gets G 0 = %(0 E − 10 aR T04 ) − %0 ($i F i =c)] − %3 '0 [(1 + $i $j fij )E − (2$i F i =c)]:
(93a)
Likewise, assembling Eqs. (90b), (91b) and (92a), one
(93b)
Eqs. (93a) and (93b) are identical to Eqs. (54a) – (54d) if 0 is interpreted as 0F , 10 as 10P , and '0 as (0F − 10E ), which is correct because we have assumed the opacity is gray. Thus: (a) one can con
(94)
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
95
Doing so would account for the major e8ects of boundaries, and low opacity in some bands versus high opacity in others. 10. Overall iteration procedure Suppose the time-evolution of a radiating =ow has been followed to time-level t n , and all the hydrodynamic variables, radiation energy density and =ux, Eddington tensor elements, and radiation spectra, at all grid points, are known. Using radiation spectra and Eddington tensor elements at t n as
(i); n+1 n+1 and 0F new estimates of the radiation energy and =ux spectra, hence new values for 10E via Eqs. (66) and (67). These are used in Eqs. (54a) – (54d) to get new source=sink terms G 0; n+1 and G i; n+1 for the radiation-hydro equations, which can then be updated for a better estimate of the solution at t n+1 . In principle this process should be sub-cycled to convergence.
11. Conclusions We have outlined the major steps needed to perform a physically realistic calculation of radiation hydrodynamics using lab-frame radiation quantities. The greatest computational e8ort goes into computing the angle-frequency information needed to: (1) evaluate elements of the Eddington tensor, and (2) to evaluate energy and =ux spectra needed to calculate mean opacities source=sink terms. We argue that angular information can be generated for all angles and frequencies in parallel, neglecting velocity e8ects. In contrast, realistic radiation energy and =ux spectra require a solution of the frequency-dependent comoving-frame radiation energy and momentum equations. For this work we propose an algorithm based on opacity distribution functions, and an approximate closure of the comoving-frame radiation momentum equation. Use of lab-frame radiation variables in the radiation hydro equations makes sense only if one uses a
96
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
this approach. Use of lab-frame radiation quantities can also encounter numerical di5culties in the di8usion limit for extremely optically-thick media because the di8usive =ux gets swamped by the advective =ux, even at small values of v=c. However, we note that this di5culty may not occur in neutron-transport problems where very large scattering-thicknesses are not encountered (for non-critical systems). Comoving-frame variables are most often used in Lagrangian-mesh calculations. One then avoids interface tracking, though zones with subscale mixing still present problems. In the past, Lagrangian-mesh computations had trouble with energy and momentum conservation. But recently Caramana and his associates have developed methods [28–31] that give precise total energy conservation (computer-word accuracy), and resist deformation and=or tangling (producing “boomerang” and “bowtie”, or “hourglass”, cells), on Lagrangian meshes. The computational e8ort for a radiation hydro step is about the same whichever radiation frame is used, being dominated by construction of the radiation
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
97
[10] Munier A, Weaver R. Radiation transfer in the =uid frame: a covariant formulation. Part II: the radiation transfer equation. Comput Phys Rep 1986;3:165–208. [11] Munier A, Weaver R. Radiation transfer in the =uid frame: a covariant formulation. Part I: radiation hydrodynamics. Comput Phys Rep 1986;3:125–64. [12] Dor< E, Drury L. Simple adaptive grids for 1-D initial value problems. J Comput Phys 1987;69:175–95. [13] Winkler K-H, Norman ML, Mihalas D. Implicit adaptive-grid radiation hydrodynamics. In Multiple timescales: computational techniques. vol. 2. New York: Academic Press, 1985 [Chapter 6] p. 145 –84. [14] Winkler K-H, Norman ML, Mihalas D. Adaptive-mesh radiation by hydrodynamics. I. The radiation transport equation in a completely adaptive coordinate system. JQSRT 1984;31:473–8. [15] Mihalas D, Winkler K-H, Norman ML. Adaptive-mesh radiation hydrodynamics. II. The radiation and =uid equations in relativistic =ows. JQSRT 1984;31:479–89. [16] Lowrie R, Mihalas D, Morel JE. Comoving-frame radiation transport for nonrelativistic =uid velocities. JQSRT 2001;69:291–304. [17] Synge JL. The relativistic gas. Amsterdam: North-Holland, 1957. [18] Hazlehurst J, Sargent WLW. Hydrodynamics in a radiation
On laboratory-frame radiation hydrodynamics Dimitri Mihalas ∗ , L.H. Auer Group X-6, MS-D409, Applied Physics Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA Received 28 September 2000
Abstract We derive and discuss the equations of radiation hydrodynamics with all radiation quantities expressed in the laboratory frame. Relativistically exact transformations from the comoving to laboratory frame are used. We obtain simple, fully covariant, expressions for the radiation energy and momentum source=sink terms using lab-frame radiation quantities. The mathematical simplicity of having lab-frame radiation quantities in conservative operators on the left-hand side is maintained. The resulting system can be solved e5ciently using an Eddington tensor formulation for angular dependences, and an iterative frequency-splitting, in the spirit of the multifrequency-gray method. Use of lab-frame radiation quantities can present numerical di5culties in making an accurate connection to the di8usion limit for extremely optically thick media; but it may have no di5culty for neutron-transport problems where very large scattering-thicknesses are generally not encountered (for non-critical systems). Published by Elsevier Science Ltd. Keywords: Lab-frame radiation hydrodynamics; Opacity distribution functions
1. Introduction In this paper, we discuss radiation hydrodynamics when all radiation quantities are evaluated in the laboratory frame, i.e. the
0022-4073/01/$ - see front matter. Published by Elsevier Science Ltd. PII: S 0 0 2 2 - 4 0 7 3 ( 0 1 ) 0 0 0 1 3 - 9
62
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
energy and momentum equations. In Section 4 we combine the radiation energy and momentum equations with nonrelativistic hydrodynamical equations to get a system for nonrelativistic radiating =ow. We o8er alternative expressions for the momentum and energy equations (see also Refs. [2,3]). In Section 5, we quote Thomas’s transformation laws [4] for the radiation
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
63
departures on the micro-scale from the strict isotropy of a Maxwellian distribution usually can be accounted for with transport coe5cients (conductivity, viscosity) calculated from, say, Chapman–Enskog theory. In this limit one can analytically integrate over phase space to derive explicit expressions, depending only on local values of variables and their gradients, for all macroscopic =ow quantities appearing in the hydro equations. Thus for pure hydro (no radiation) the system is at most four-dimensional. (For free molecule :ow, where particle mean-free-paths become comparable or large compared to gradient lengths, it is necessary to include nonlocal e8ects in the distribution function by solving the full Boltzmann equation. We do not discuss such problems this paper.) Therefore for material, the comoving frame is natural because all physical properties can be de
64
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
density because of advection. The mathematical consequence of these physical e8ects is that photons, moving with >xed laboratory-frame angle and frequency, are observed to move on a curved trajectory with varying frequency in the comoving frame; hence partial derivatives with respect to these variables appear in the di8erential operator describing the transport process. (In geometric language, the di8erential operator must account for changes in the basis vectors of the tangent spaces attached to the space-time manifold along the null geodesic of a photon.) The resulting transport equation is seven-dimensional (4 space-time variables, 2 angle cosines of a photon’s direction, and frequency). Explicit expressions for the comoving-frame transport operator may not be available in the literature for a particular geometry or symmetry. The comoving-frame transport equation, correct to O(v=c) for nonrelativistic 1D =ows in spherical symmetry is derived in Ref. [6]. An exact comoving equation for relativistic =ows in 1D spherical symmetry is given in Ref. [7]. The comoving equation, correct to O(v=c) is written for general curvilinear, orthogonal, coordinates in Ref. [8], and generalized to relativistic =ows in Ref. [9]. A general, coordinate-free, relativistic, expression for the comoving equation is given in Ref. [10] for arbitrary geometries. But to obtain explicit expressions for a particular geometry from the results in Refs. [8–10] may require considerable algebra. Explicit expressions are given in Ref. [10] for spherically symmetric, 1D, relativistic =ows, and to O(v=c) in: (1) coordinate-free form, (2) Cartesian coordinates and 1D slab symmetry, (3) 1D spherical coordinates and symmetry, (4) 2D cylindrical coordinates and symmetry, and (5) in conservative coordinate-free form. Monochromatic and total radiation energy and momentum equations, obtained by integrating over angles or both angles and frequency, are given in all of the papers just cited; see particularly [11]. In contrast, in lab-frame transport equation, photons are seen to move on straight-line paths with
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
65
di5culty for a laboratory-frame calculation is that in the di8usion regime, the di8usive =ux, which depends on the local temperature gradient, may be dwarfed (by many orders of magnitude) in the total radiant =ux seen transported by a
66
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
spatial derivatives, and the Lagrangian time-derivative, which tracks the motion of the =uid, namely: D(•) @(•) ≡ + v · ∇(•): Dt @t
(1)
Many writers would call this description “Lagrangian”. Pomraning calls it a “modi
(2)
Clearly when vgrid ≡ 0 one recovers the Eulerian time derivative, and when vgrid ≡ v=uid one recovers the Lagrangian time derivative. We stress again that these transformations are kinematical, and Galilean. Alternatively, one can write the transport equation and its moments (dynamical equations) for comoving-frame radiation quantities, on a
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
67
3. The radiation equations 3.1. The transport equation One can describe the radiation
(3)
Here (I=s) is an absolute (or intrinsic) derivative, which is de
(4)
But Rx = n ds, and the time to cross the element is Rt = ds=c, hence @I 1 @I + ni i = − I: @x c @t
(5)
We use the summation convention on repeated indices throughout this paper. Eq. (5) holds only in Cartesian coordinates. In general, the transport equation is a six- or seven-dimensional system: a four-dimensional spacetime, two angles, and perhaps frequency. A direct solution of such dimensionality is daunting, so one seeks simpli
68
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
3.2. The radiation energy equation De
E = monochromatic radiation energy density in lab frame ≡ c I d!; F i = monochromatic radiation energy =ux in lab frame ≡ I ni d!; ij −1 P = monochromatic radiation pressure in lab frame ≡ c Ini nj d!: Integration of the transport equation over solid angle ! gives the lab-frame monochromatic radiation energy equation @E(x; t; ) @F i (x; t; ) = [(x; t; n; ) − (x; t; n; )I (x; t; n; )] d! (6) + @t @xi and integration over frequency gives the total radiation energy equation @E(x; t) @F i (x; t) + @t @xi ∞ = d d![(x; t; n; ) − (x; t; n; )I (x; t; n; )] ≡ −cG 0 : 0
(7)
Here cG 0 is the (rate of radiation energy absorbed by the material) − (rate of radiation energy emission by the material) per unit volume. All quantities in Eqs. (6) and (7) are measured in the lab frame. 3.3. The radiation momentum equation Similarly, integration of the transport equation against ni d! over solid angle gives the lab-frame monochromatic radiation momentum equation 1 @F i (x; t; ) @P ij (x; t; ) + c2 @t @xj 1 = [(x; t; n; ) − (x; t; n; )I (x; t; n; )]ni d! c
(8)
and integrating over frequency, gives the total radiation momentum equation 1 @F i (x; t) @P ij (x; t) + c2 @t @xj ∞ 1 = d d![(x; t; n; ) − (x; t; n; )I (x; t; n; )]ni ≡ −G i : c 0
(9)
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
69
Here G i is the rate of (deposition of the ith component of the radiation momentum into the material) − (emission of the ith component of the radiation momentum by the material) per unit volume. 3.4. Conservation laws The di8erential operators in the lab-frame radiation and momentum equations are exact, and are in conservation form. Their apparent simplicity is, however, compromised by the di5culty of evaluating accurately the integrals on the right-hand sides in the lab frame. As described in Section 2 the principal di5culty is that for moving material, radiation moving with a
(10)
(11)
Because of symmetry f has 1, 3, and 6 nonredundant components in 1D, 2D, and 3D, respectively. The elements of f depend upon (x; t; ). We suppress explicit mention of these variables to emphasize that, at a given stage of the calculation, f is a collection of angular factors computed from an angle-frequency-dependent transport equation with given source=sink terms, as discussed in Section 9. In practice this splitting works because the elements of f are ratios, which re=ect primarily the presence of open boundaries. Beyond one optical depth from a boundary, values of fij can be accurate even when current estimates of the source=sink terms used in the transport equation are not yet well determined by the outer iteration procedure which gives material properties, and E and F, that satisfy the radiation hydro equations at the advanced time. But we emphasize that the solution of the radiation-dynamical equations can be no better than the accuracy of the tensor elements used to close the system. Thus if one merely assumes
70
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
the lab-frame Eddington tensor is isotropic in a moving =uid, then the radiation stress is forced to act like a hydrostatic pressure (like the gas pressure). This assumption would be valid in the comoving frame if "p =‘1, (where "p is a photon mean free path and ‘ is a characteristic gradient length in the =uid). But one needs to Lorentz transform to the lab frame. And if there is an open boundary, or the opacity of the material is small at some frequencies, this approximation may not be valid even in the comoving frame. When f in the lab frame has nonzero o8-diagonal elements, radiation stresses can induce shear in the =ow. A physically correct solution requires Eddington-tensor elements that account fully for the problem’s geometry, the frequency dependence of photon mean-free-paths, and the e8ects of the =uid’s velocity; we address some aspects of this problem in Section 9. 4. Lab-frame radiation hydrodynamics equations In nonrelativistic laboratory =ows, where $ ≡ v=c1, the
(12)
(13)
(14)
(c) the material energy equation @(e0 ) @vi @ + i (e0 vi + q0i ) + (p0 ij − '0ij ) j = *0 + cG 0 − vi G i @t @x @x
(15)
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
or the radiating =uid energy equation @ @ @vi (e0 + E) + i (e0 vi + q0i + F i ) + (p0 ij − '0ij ) j = *0 − vi G i @t @x @x or the material total energy equation @ @ [( 12 v2 + e0 )] + i {[(( 12 v2 + e0 ) + p0 ]vi − '0ij vj + q0i } = vi fi + *0 + cG 0 ; @t @x or the radiating =uid total energy equation @ @ [( 12 v2 + e0 ) + E] + i {[(( 12 v2 + e0 ) + p0 ]vi − '0ij vj + q0i + F i } = vi fi + *0 ; @t @x (d) the radiation energy equation @E @F i + i = −cG 0 @t @x and (e) the radiation momentum equation 1 @F i @P ij + j = −G i : c2 @t @x
71
(16)
(17)
(18)
(19)
(20)
In Eqs. (12) – (18), p0 is the gas pressure, e0 the =uid’s internal energy, '0ij its viscous stress tensor, q0i the thermal heat conduction vector, and *0 the rate of nuclear energy release per unit mass. By solving these equations one gets ; v; T0 ; E, and F. The di8erential operators are all in conservation form, which makes them attractive computationally. The two basic problems still to be faced are computation of elements of the Eddington tensor f, and lab-frame values of the four-force G + . We address these issues in Sections 9 and 6, respectively. 5. Transformation laws 5.1. Transformation of the photon four-momentum Denote comoving-frame quantities with a subscript “0” and leave lab-frame quantities unadorned. The connection between and 0 , and n and n0 , is given by Lorentz transformation of the photon four-momentum (h=c)(1; n). Transforming from the =uid frame to the lab frame one
n =
0
ni0
%$j nj0 + 1+ %+1
(21a)
i
%$ :
(21b)
The inverse transformation is 0 = %(1 − $i ni )
(22a)
72
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
and
%$j nj = n − 1− %$i : %+1 Here $i ≡ vi =c; $2 ≡ $i $i , and % ≡ 1= 1 − $2 .
ni0
0
i
(22b)
5.2. L.H. Thomas’s transformation laws In a seminal paper, Thomas [4] showed by simple, but genuinely covariant, Gedanken experiments, I (n; ) = (3 =30 )I0 (n0 ; 0 ) ≡%
−3
i −3
(1 − $i n )
(n; ) = (2 =20 )0 (0 ) ≡
I0
n − %(1 − %$i ni =(% + 1))R i ; %(1 − $i n ) ; %(1 − $i ni )
0 [%(1 − $i ni )] %2 (1 − $i ni )2
(23) (24)
and (n; ) = (0 =)0 (0 ) ≡ %(1 − $i ni )0 [%(1 − $i ni )]:
(25)
Notice that isotropy of the material properties 0 and 0 in the comoving frame is assumed. The second form of Eqs. (23) – (25) displays explicitly the angular dependences that result from transforming from the comoving frame to the lab frame. Eqs. (23) – (25) say that I=3 ; =2 , and , are Lorentz invariants. In addition, for any particle, d 3 p= e; ˜ where d 3 p = dpx dpy dpz = p2 dp d! is the three-space di8erential of the particle’s momentum, and e˜ is the particle’s energy, is also a Lorentz invariant (see [3, Eq. (43.17)]; [5, Eq. (A.65)]). Hence for photons d d! = 0 d0 d!0
(26)
d! = (20 =2 ) d!0 :
(27)
and We will use these results repeatedly. 5.3. Covariance of radiation four-force and stress-energy tensor Although I; E; F; P; , and , have all been taken to be in the lab frame, it should be noted that G + ≡ (G 0 ; G) ≡ (G 0 ; G 1 ; G 2 ; G 3 )
(28)
is a four-vector, which can be transformed from frame to frame by Lorentz transformation. One can verify this statement by using the relations given in Section 5.2. One sees that each
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
73
component of G + as given by Eqs. (7) and (9) is (c=h) times a component of the photon four-momentum M + ≡ (M 0 ; M) ≡ (h=c)(1; n)
(29)
times either the Lorentz invariant emissivity =2 , or the Lorentz invariant opacity times the invariant speci
%$j G j i =G −% G − $: %+1 The inverse transformation is
G0i
i
0
G 0 = %(G00 + $i G0i ) and
j G %$ j 0 G i = G0i + % G00 + $i : %+1
(30a) (30b) (31a)
(31b)
Likewise, the left-hand sides of these equations are c times the four-divergence of the 4 × 4 radiation stress-energy tensor R, where ∞ R+$ ≡ c−1 d d! I (n; )n+ n$ : (32) 0
Here one uses the convention that n0 = 1. In tensor form E c−1 F† R≡ : c−1 F P
(33)
R is a symmetric four-tensor, which can be transformed into the comoving frame by Lorentz transformation. One can verify this statement by noting each element of R is a component of the
outer product of c=h times the photon four momentum with itself, times the invariant speci
(34a)
F0i =c = %{(F i =c) − %$i E − $j P ij + $i [%=(% + 1)][(2% + 1)($j F j =c) − %($j $k P jk )]}
(34b)
and P0ij = P ij − %[($i F j + $j F i )=c] + %2 $i $j E + [%2 =(% + 1)][$j $k P ik + $i $k P kj − 2%$i $j ($k F k =c)] + [%2 =(% + 1)]2 $i $j ($k $l P kl ):
(34c)
74
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
The same results are found by using Eqs. (21a), (21b), and (26) in the
(35a)
F i =c = %{(F0i =c) + %$i E0 + $j P0ij + $i [%=(% + 1)] j
jk
×[(2% + 1)($j F0 =c) + %($j $k P0 )]}
(35b)
and P ij = P0ij + %[($i F0j + $j F0i )=c] + %2 $i $j E0 + [%2 =(% + 1)][$j $k P0ik + $i $k P0kj + 2%$i $j ($k F0k =c)] + [%2 =(% + 1)]2 $i $j ($k $l P0kl ):
(35c)
Covariant (both in coordinates and relativistically) radiation dynamical equations are @R+$ = −G + : (36) @x$ The total radiation energy and momentum equations are direct analogues, for radiation, of the energy and momentum equations for the material part of the =uid. As mentioned before, they are exact. But to apply them one still needs Eddington tensor elements, and usable expressions for the four-force. 6. Lab-frame radiation source= sink terms 6.1. Comoving-frame source=sink terms Consider
(37)
Likewise for the emissivity, one discriminates between (the isotropic) radiation emitted thermally, and radiation scattered from the incident radiation
(38)
For brevity, explicit mention of the dependence of 0 ; 10 ; '0 , and 0 on x and t has been suppressed.
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
75
Assuming LTE, one invokes Kirchho8’s law to write the thermal emissivity as t0 (0 ) = 10 (0 )B(0 ; T0 );
(39)
where B(0 ; T0 ) is the Planck function, which is isotropic. The scattering source is more complicated. Let the redistribution function R(n0 ; 0 ; n0 ; 0 ) be the probability that an incoming photon at (n0 ; 0 ) in the comoving frame emerges in the comoving frame at (n0 ; 0 ). R is normalized such that ∞ d!0 d!0 ∞ d0 d0 (40) R(n0 ; 0 ; n0 ; 0 ) = 1: 43 0 43 0 Then the scattering emissivity is ∞ d!0 s0 (n0 ; 0 ) = '0 d0 R(n0 ; 0 ; n0 ; 0 )I0 (n0 ; 0 ): 43 0
(41)
Here '0 is the total scattering cross-section. Redistribution functions can be calculated analytically for various models of redistribution in angle and frequency by the scattering process in an atom’s frame. Scattering shows another aspect of the di5culty in solving the transfer equation. Not only is it a 7D partial di8erential equation in spacetime, angles, and frequencies, but it is an integro-partial di8erential equation in both the input and output angle-frequency spaces. One must therefore seek simpli
(42)
(43)
Thomson scattering has a mild dipole phase factor, which, for simplicity, we take as isotropic, g ≡ 1, and coherent, R(0 ; 0 ) ≡ (0 − 0 ). Hence d!0 c'0 (0 ) c'0 (0 ) s 0 (0 ) = d0 (0 − 0 ) (44) I0 (n0 ; 0 ) = E0 (0 ): 43 43 43 In this simple case s0 (0 ) is also isotropic. The total comoving-frame emissivity is 0 (0 ) = t0 (0 ) + s0 (0 ) = 10 (0 )B(0 ; T0 ) + '0 (0 )[cE0 (0 )=43]:
(45)
We adopt this simpli
76
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
6.2. Comoving-frame energy and momentum input terms Calculating the net radiant energy input and momentum input into the material, measured in the comoving frame, one
and G0i
≡c
−1
∞
0
d0
d!0 [0 (x; t; 0 )I0 (x; t; n0 ; 0 ) − 0 (x; t; 0 )]ni0
and taking advantage of the isotropy of 0 and 0 in this frame, one has ∞ 0 cG0 = 43 10 (x; t; 0 )[(c=43)E0 (x; t; 0 ) − B(x; t; 0 ; T0 )] d0 0
and G0i = c−1
0
∞
(x; t; 0 )F0i (x; t; 0 ) d0 :
(46b)
(47a)
(47b)
Both G00 and G0i are simpler than G 0 and G i because (a) they contain only moments of the radiation
(49)
Here E0 (0 ) ≡ E0 (0 )=E0 is the comoving-frame radiation energy spectrum, which is normalized: 0∞ E0 (0 ) d0 = 1. Then cG00 = c10E E0 − 4310P B(T0 ) = c(10E E0 − 10P aR T04 ):
(50)
Here aR is Stefan’s radiation constant. Both 10E and 10P are scalars. Whereas 10P is a function of density and temperature only, evaluation of 10E requires knowledge of the comoving-frame
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
77
radiation energy spectrum E0 (0 ). If one can actually evaluate this spectrum, then no approximation is made by using 10E and 10P in G00 . (i) Similarly, for G i0 , de
where F(i) (0 ) ≡ F0i (0 )=F0i is the comoving-frame radiation :ux spectrum in the ith direction, which also is normalized: 0∞ F(i) 0 (0 ) d0 = 1. Then (i) i G0i = 0F F0 =c:
(52)
(i) is not a function of (0 ; T0 ) only; nor is it, in general, a scalar, because the =ux Like 10E ; 0F (i) spectrum can be di8erent along each axis. One should regard 0F as a 3 × 3 diagonal matrix. (i) i Note that F0 along a given axis may be zero, so that 0F in that direction is unde
6.4. Lab-frame energy and momentum input terms In the laboratory frame, the net radiant energy and momentum input to the material is given by Eqs. (7) and (9). It is not obvious that one can express G 0 and G i solely in terms of radiation moments in that frame. And, as mentioned above, accurate evaluation of the double integral in Eq. (9) can be di5cult in opaque material because of the near cancellation of the intensity along rays in opposite directions. One can deal with these problems as follows:
and i
G =
(i) i (0F F0 =c)
+ % 10E E0 −
10P aR T04
% (j) j + $j (0F F0 =c) $i : %+1
(53a)
(53b)
The terms containing the temperature result from thermal emission. We thus have lab-frame energy and momentum input terms expressed in terms of comovingframe moments, and frequency integrals evaluated in the comoving frame. These expressions are exact to all orders of (v=c); the only simpli
78
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
Eqs. (53a) and (53b) are still not suitable for use in Eqs. (7) and (9) because lab-frame radiation variables would appear on the left-hand sides, and comoving-frame quantities on the right. However, because only frequency-integrated moments appear in Eqs. (53a) and (53b) one can remove the discordance by expressing E0 and F0 in the lab frame by using Eqs. (34a) – (34c). (Note that T0 and the opacities do not require re-transformation because they are world scalars (see Refs. [17, p. 29]; [18]). Carrying out these substitutions, and grouping terms in ascending orders of $, one
− [%4 =(% + 1)]$i 0F $i ($j $k fjk )E:
(54a)
(i) can be replaced Or, assuming that the =ux spectrum is the same in all directions, so that 0F by a scalar 0F ,
G 0 = %[%2 10E + (1 − %2 )0F ]E − %10P aR T04 − %($i F i =c)[0F − 2%2 (0F − 10E )] − %3 (0F − 10E )($i $j fij )E:
(54b)
Similarly, (i) i (i) G i = %(0F F =c) − {[%0F ($j fij + %$i ) − %2 10E $i ]E − %10P aR T04 $i }
+
%2 (i) (j) {[(2% + 1)0F − %10E ]($j F j =c) + %($j 0F F j =c)}$i %+1
+
%3 (i) (j) {[(% + 1)10E − 0F ]($j $k fjk ) − $j 0F (%$j + $k fjk )}$i E %+1
+
%4 ( j) j ($j 0F $ )[(2% + 1)($k F k =c) − %$k $l fkl E]$i : (% + 1)2
(54c)
Or, again assuming a scalar 0F can be used, G i = %0F (F i =c) − %10P aR T04 $i − [%3 (0F − 10E )$i + %0F $j fij ]E + %3 (0F − 10E )[2($j F j =c) − ($j $k fjk )E]$i :
(54d)
Eqs. (54a) – (54d) are covariant, both geometrically and relativistically, and provide workable expressions for (G 0 ; G) correct to all orders of v=c. All the radiation variables, including Eddington tensor elements, are measured in the lab frame. In contrast, the material variables (i) 0 ; T0 ; 10P ; 10E ; and 0F are all measured in the comoving frame.
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
79
6.5. Check of formulae by using thermal radiation One can check the algebra required to reduce Eqs. (53a) and (53b) to Eqs. (54a) – (54d), by considering an in
(55a)
F i =c = − 43 %2 $i aR cT04
(55b)
P ij = 13 (ij + 4%2 $i $j )aR T04 :
(55c)
and
Eq. (55b) shows that the radiation =ux has a dipole component in the lab frame, proportional to the negative of the lab’s velocity relative to the =uid frame. The o8-diagonal components of the lab-frame radiation-pressure tensor are O(v2 =c2 ). Even if radiation is isotropic in the =uid frame, it is not in the lab frame. In fact, the observed dipole component of the cosmic microwave background has allowed determination of the motion of our galaxy relative to the “rest frame of the Universe”. From Eqs. (55a) and (55c) the lab-frame Eddington tensor elements for thermal radiation are fij ≡ P ij =E = (ij + 4%2 $i $j )=(4%2 − 1);
(56a)
whence one
(56b)
$i $j fij = (4%2 − 3)$2 =(4%2 − 1):
(56c)
and
As noted above, the o8-diagonal terms of the Eddington tensor are only O(v2 =c2 ). Nevertheless, at high temperatures, radiation pressure dominates gas pressure by orders of magnitude, hence even a small anisotropy in the radiation pressure may have important consequences for the =ow dynamics. If one assumes strict thermal equilibrium in the rest frame of the =uid, one knows that G00 ≡ 0, and G0 ≡ 0. Then the principle of equivalence between inertial frames states that an observer in the lab frame must obtain the same result: i.e. the =uid is observed to be in thermal equilibrium and is not accelerated by radiation forces. Hence one must
80
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
6.6. Lab-frame radiation dynamical equations Use of Eqs. (54a) – (54d) in Eq. (36) gives a
0 @R+$ (1=c)(@E=@t + @F j =@xj ) G E = =− =A · + B: 1=c2 @F i =@t + @(fij E)=@xj Gi F j =c @x$
(57)
Here + runs from 0 to 3, i and j run from 1 to 3, A is a 4 × 4 matrix, and B is a vector of length 4. This system has an exact lab-frame (four-) divergence on the left-hand side, and a linear combination of the lab-frame energy density, =ux components, and Eddington tensor, correct, in principle, to all orders of (v=c) on the right. The actual accuracy of the system will depend, of course, on the accuracy of the spectral pro
7. Evaluation of mean opacities and radiation spectra Consider now the evaluation of mean opacities. The opacity spectrum of elements other than hydrogenic or helium-like ions is complicated, and a million frequency sample-points may be needed to characterize it adequately for mixtures of elements. Direct calculation of the integrals in Eqs. (49) and (51) is then prohibitive. A number of approximations, ranging from over-simpli
(58)
G i = 1[(F i =c) − (aR T04 + 13 E)$i ];
(59)
and
which agree with the third and fourth un-numbered equations following Eq. (36) in Ref. [19]. Other approximate source terms based on less restrictive assumptions can also be deduced from Eqs. (54a) – (54d). But no real material is gray, and even if this assumption is useful in studying the mathematical behavior of the equations it is an inadequate representation of reality.
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
81
We can still gain some important insight from Eq. (59); using it in the radiation momentum equation we have, to O(v=c), "p @F i @P ij (60) + c"p j = −F i + 43 Evi : c @t @x There are three characteristic timescales to be considered: the radiation-=ow time tr ∼ "p =c, the =uid-=ow time tf ∼ ‘=v, and the di8usion time td ∼ (‘="p )2 ("p =c) = ‘2 ="p c. The ratio of the time derivative of the =ux to the =ux itself is O("p =‘); O("p v=‘c), and O("p2 =‘2 ) on the radiation-=ow, =uid-=ow, and di8usion timescales, respectively. For di8usion theory to apply, "p =‘1, hence the time derivative can always be neglected in a very opaque material. Therefore, in the di8usion regime F i = −c"p
@P ij 4 i + 3 Ev = F0i + 43 Evi ; @xj
(61)
which is Eq. (35b) to
82
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
Fig. 1. Upper panel: Planck function, normalized to its maximum value, at T0 = 0:02 KeV, versus frequency in KeV. Lower panel: Spectrum of iron at T0 = 0:02 KeV and = 10−4 gm=cm3 , versus frequency in KeV. Vertical dotted lines illustrate division of frequency into
of the Planck function B(; T0 ). As shown in Fig. 1, this variation is easily resolved. From tables of the Planck function (e.g. Ref. [20, pp. 106 –107]) one sees that a range of fewer than three decades in u ≡ h=kT0 , namely 0:2 6 u 6 25, is su5cient to cover the full variation of B(; T0 ). One then divides the spectrum into bands, each bounded by frequencies (b ; b+1 ); b = 1; : : : ; B. The band sizes :b ≡ b+1 − b need not be equal, but, once chosen, must be the same for all values (0 ; T0 ) in a table. Each band should be narrow enough that the Planck function varies only mildly within it, so that b+1 Bb (T0 ) ≡ B(0 ; T0 ) d=:b (62) b
is well de
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
83
Fig. 2. Planck function normalized to its maximum value (upper panel), and spectrum of iron (lower panel) at T0 = 0:02 KeV and = 10−4 gm=cm3 , versus frequency on the range 0:04 . . 0:075 KeV.
The second source of variation is from the variation of the material’s opacity with energy. As seen in Fig. 1, which is for Iron at T0 = 0:02 KeV and = 10−4 g=cm, this variation is extremely complicated and much more di5cult to treat than the smooth variation of the Planck function. A subset of these data on the range 0:04 . . 0:75 is shown in Fig. 2. Although the Planck function varies by only 10% over this range, the opacity varies by orders of magnitude over the width of each spectral line. To resolve this variation fully we would need O(:b =:D ) frequency points in a band, where :D is the Doppler width of a typical line. At this resolution we would require hundreds of thousands of frequency points per band. The treatment of a frequency spectrum in this much detail in time-dependent problems is beyond the capability of any existing computer. Yet, the frequency variation of the opacity is a de
84
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
Currently, the technique proven to be the best method for treating complex spectra in astrophysical work [21,22] is the opacity distribution function. Once the bands have been de
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
85
Fig. 3. Opacity distribution function with nine ODF pickets for the Iron spectrum shown in Fig. 2.
7.3. Full transport Consider
(64a)
or, if the scattering terms are canceled out analytically, i @Epb @Fpb = 1pb [43Bb (T0 ) − cEpb ] + @t @xi
(64b)
and from Eq. (8) the radiation momentum equation is ij i pb i 1 @Fpb @(fpb Epb ) =− F ; + j 2 @x c pb c @t
(65)
86
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
(p = 1; : : : ; Pb ); (b = 1; : : : ; B). If we know the Eddington tensor elements, from a separate i . Suppose now the medium ray-by-ray calculation, these equations can be solved for Epb and Fpb is time-independent. At a point far from a boundary, and where ("p =‘1) in all pickets, Epb → Bb (T0 ), and strict energy balance is guaranteed in Eq. (64b). Fortunately, the ODF procedure is strongly convergent with increasing number of ODF pickets. i for all p and b, one can compute Solving Eqs. (64b) and (65) to determine Epb and Fpb approximate energy and =ux means needed in the integrated moment equations: B Pb b=1 Rb p=1 wpb 1pb Epb (66) 10E = B Pb b=1 Rb p=1 wpb Epb and (i) 0F
B
=
Pb i b=1 Rb p=1 wpb pb Fpb : B Pb i b=1 Rb p=1 wpb Fpb
(67)
Having 10E and 0F one can evaluate the source terms in Eqs. (54a) and (54c), and solve system (12) – (20). (i) These means are vastly superior to use of 10P for 10E and 0R for 0F , or use of two-temperature means, because they account for: (1) the huge di8erences in opacity between lines and continua, (2) decoupling of the local monochromatic energy density from the Planck function, (3) anisotropy of the radiation
8. Velocity gradient e/ects Using an ODF representation of opacity one can obtain accurate comoving-frame energy and =ux spectra, and compute the mean opacities in equations (54). But the e8ects of Doppler shifts, resulting from velocity gradients, have thus far been ignored. To get a rough estimate, one might compare the Doppler shift over a gradient length ‘, or a photon mean-free-path "p , namely ∼ (=c)(dv=ds) min(‘; "p ), with :sb , the interval between sample points in the opacity spectrum. If :sb , true in very opaque material, one might argue that Doppler shifts are negligible. But caution is needed even then, because they produce red- or blue-shifting of the radiation’s spectrum under adiabatic expansion or compression [23,3, p. 474]. Near boundaries (where "p → ∞), or where the velocity gradient is large, Doppler shifts can be quite large. The monochromatic comoving-frame radiation energy equation is
@F i (0 ) D E0 (0 ) vi @F i (0 ) + 0 i = 10 (0 )[43B(0 ; T0 ) − cE0 (0 )] − 2 0 Dt @x c @t ij @[f0 E0 (0 )] @vi : (68) + 0 @0 @xj
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
87
The dominant terms of Eq. (68) are written on the
@[f0ij E0 (0 )] vj @[f0ij E0 (0 )] D F0i (0 ) 0 (0 ) i 1 @vi j + = − ( ) − − F (0 ) F 0 0 c2 Dt @xj c c2 @t c2 @xj 0 ijk 1 @[0 Q0 (0 )] @vj : (69) + 2 c @0 @xk Eqs. (68) and (69) are accurate only to O(v=c), which is su5cient to determine spectra. In Eq. (69), Q0ijk (0 ) ≡ I0 (n0 ; 0 ) ni0 nj0 nk0 d!0 : (70) One can write an approximate closure relation for this three-index tensor by assuming a P1 approximation for the radiation
(71)
Q0ijk (0 ) ≈ ( 15 + 25 ij )jk F0i (0 ):
(72)
Using Eq. (72) in Eq. (69), one
@[f0ij E0 (0 )] vj @[f0ij E0 (0 )] D F0i (0 ) 0 (0 ) i 1 @vi j + = − ( ) − − F (0 ) F 0 0 c2 Dt c c2 c2 @xj 0 @xj @t 1 @[0 F0i (0 )] @v j @v(i) + 2 + 2 : (73) 5c @0 @xj @x(i) In Eq. (73), the parentheses on the index i in the last term mean no summation on i. The dominant terms are the divergence of P0ij on the left-hand side and the =ux-attenuation term on the right-hand side. Compared to the latter, the time derivative of F0i on the left-hand side is O("p v=‘c) on a =uid-=ow timescale, and O("p =‘) on a radiation-=ow timescale, hence it is important only in regions where the medium is transparent. On a radiation-=ow timescale, the term containing the time derivative of P0ij on the right-hand side is O(v=c) compared to the divergence of P0ij on the left-hand side, and is O(v2 =c2 ) on a =uid-=ow timescale. The term containing the =ux times the velocity-gradient =ux on the right-hand side is always O("p v=‘c) compared to the =ux attenuation term. The frequency derivative of F0i is O("p v=‘c) compared to the =ux-attenuation term, thus the frequency redistribution of the =ux can also be treated iteratively.
88
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
In writing
ij ij @(fpb Epb ) pb i vj @(fpb Epb ) D Fpb 1 @vi j + = − − − F F c2 Dt @xj c pb c2 @t c2 @xj pb b i i 1 ( p vpp Xp b Fp b − Xpb Fpb ) @v j @v(i) + 2 + 2 : (75) 5c Rs @xj @x(i) Here it is assumed that time and frequency di8erencing is done in the usual way. As pointed out above, the frequency derivatives on the right-hand side of Eqs. (74) and (75) are small compared to terms on the left-hand side. Thus these equations need not be treated as coupled in frequency, but can be solved for each frequency, in parallel, using previous estimates of Epb i i , and Fpb in the frequency di8erence terms. This process yields new values for Epb and Fpb which can then be used in the frequency di8erences, and the procedure iterated to convergence. i , the Our expectation is that the process will converge rapidly. Given values for Epb and Fpb mean opacities can be computed from Eqs. (66) and (67). We mention in passing that the structure of the r and v matrices can provide guidance about the quality of the ODF discretization. For example, if the original spectrum is over-resolved, then most opacity sample-points binned into a given ODF picket will likely couple to other points binned into the same picket, and both r and v will be nearly (i.e. too strongly) diagonal. Ideally these matrices should be nearly bidiagonal. On the other hand, if the sampling width :sb in the opacity spectrum is too large, then sampling of line pro
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
89
to construct correlation matrices connecting the ODF picket representations on one side of the interface with those on the other. 9. The lab-frame Eddington tensor for a moving 1uid 9.1. Analytical lab-frame Eddington tensors for simple cases It is possible to calculate the lab-frame Eddington tensor analytically for simple problems. The results may give insight into the size of velocity e8ects, and provide benchmarks for testing numerical codes. 9.1.1. The di
(76a)
F i =c = %{− 43 (@ ln T0 =@xi )= + 43 %$i − [( 43 %$i =)(2% + 1)=(% + 1)] ×(@ ln T0 =@xj )$j }aR T04
(76b)
and P ij = ( 13 ij + 43 %2 $i $j − ( 43 %=){$i (@ ln T0 =@xj ) + $j (@ ln T0 =@xi ) + [2%2 $i $j =(% + 1)][$k (@ ln T0 =@xk )]})aR T04
(76c)
hence ij
f =
(ij + 4%2 $i $j ) −
4% i 2%2 $i $j $k j j i k [$ (@ ln T0 =@x ) + $ (@ ln T0 =@x ) + (%+1) (@ ln T0 =@x )] : (4%2 − 1) − (8%2 =)(@ ln T0 =@xi )$i
(77)
The dominant o8-diagonal terms of the Eddington tensor are only O(v2 =c2 ) or O("p v=‘c) (from the di8usive =ux) relative to the diagonal, hence are generally negligible to O(v=c). The dominant di8erence between the lab-frame =ux F and the comoving frame =ux F0 is 43 v times the radiation energy density; the latter (advective) term can dominate F0 by orders of magnitude in optically thick material, hence must always be taken into account for a moving =uid. 9.1.2. 1D semi-in>nite grey medium Another problem for which an exact solution is known is a 1D, semi-in
(78a)
90
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
and F0 ≡ F0z = constant;
F0x = F0y ≡ 0:
(78b)
The radiation pressure tensor is diagonal, but anisotropic, with cPzz 0 = F0 [7 + q(∞)];
(78c)
yy 3 1 cPxx 0 = cP0 = F0 [7 + 2 q(7) − 2 q(∞)];
(78d)
xy
yz
P0 = Pxz 0 = P0 = 0:
(78e)
Here 7 is the optical depth from the boundary, q(7) is the Hopf function, and F0 is the constant =ux in the outgoing direction. An exact expression for the Hopf function was found√by Mark [24] by Laplace transform techniques. It is a monotonic function ranging from q(0)=1= 3=0:577350, to q(∞) = 0:710446, which can be represented accurately by the variational
(79)
where En (x) denotes the exponential integral of order n. Thus the exact value of the Eddington tensor in the comoving frame of a semi-in
(80a)
f0xx (7) = f0yy (7) = 13 [7 + 32 q(7) − 12 q(∞)]=[7 + q(7)]
(80b)
with all o8-diagonal entries zero. Clearly lim fii 7→∞ 0
= 13 ;
(80c)
showing that the radiation
(81a)
F z = %F0 {1 + 3%$z [7 + q(7)] + $z [7 + q(∞)] + [%$z2 (2% + 1)=(% + 1)] + [%2 $z =(% + 1)]{$2 [7 + q(∞)] + 32 ($x2 + $y2 )[q(7) − q(∞)]}}
(81b)
F x = %F0 {3%$x [7 + q(7)] + $x [7 + 32 q(7) − q(∞)] + %$x $z (2% + 1)=(% + 1) + [%2 $x =(% + 1)]{$2 [7 + q(∞)] + 32 ($x2 + $y2 )[q(7) − q(∞)]}}
(81c)
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
with an identical expression for F y with x and y interchanged, and P zz = F0 7 + q(∞) + 2%$z {1 + [%2 $z2 =(% + 1)]} + 4%2 $z2 7 %2 ($x2 + $y2 ) 3%2 ($x2 + $y2 ) 2 2 2 2 + 3% $z 1 + q(7) + % $z 1 − q(∞) : 2(% + 1)2 2(% + 1)2
91
(81d)
It is not di5cult to calculate other components of P, but these are not useful for the problem being considered, a 1D medium with constant velocity, because in this case (@=@x) and (@=@y) of all quantities are identically zero. Thus a velocity
I () d, c−1 F = c−1 I () d n (n − n) d! = (nx i + ny j + nk k)c−1 I
and −1
(n ⊗ n )(n − n) d! n2x ii nx ny ij nx nz ik = nx ny ji n2y jj ny nz jk c−1 I: nx nz ki nz ny kj n2z kk
P=c
(82b)
I () d
(82c)
Then in the =uid, moving with velocity v in the lab frame, one
(83a)
F0i =c = %(1 − $i ni )(ni − %$i {1 − [%$i ni =(% + 1)]})c−1 I
(83b)
P0ij = (ni − %$i {1 − [%$i ni =(% + 1)]})(nj − %$j {1 − [%$j nj =(% + 1)]})c−1 I:
(83c)
and
92
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
Or, recalling Eqs. (22a) and (22b), one can write, equivalently, E0 = (0 =)2 (I=c); F0i =c = (0 =)2 (I=c)ni0 , and P0ij = (0 =)2 ni0 nj0 (I=c). These expressions are valid only for a collimated beam (i.e. a delta-function in angle in the lab frame). The Eddington tensor in the moving frame has O(v=c) terms which account for the aberration of the beam’s direction as seen by the moving observer. In summary, for the problems considered above, velocity-induced e8ects in the radiation pressure tensor are only O(v2 =c2 ) or O(v"p =c‘) at depth. They are at most O(v=c) near open boundaries, but even then below the angular resolution of typical methods of solving the transport equation. For such problems velocity e8ects are thus negligible in the computation of the Eddington tensor. O(1) departures of the radiation pressure tensor from isotropy re=ect transport e8ects caused by optical geometry, such as boundaries (internal or external), voids in the material, and=or steep gradients in temperature or the radiation
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
93
For thermal emission, t0 (0 ) = 10 (2h30 =c2 )[exp(h0 =kT0 ) − 1]−1 , and from Eqs. (24) and (22a) one has 3 4 4 ∞ 3 3 ∞ 2k T0 10 x dx aR c t (n; ) d = 10 T 4 = 0 c2 h3 ex − 1 0 43 0 0 a c 1 R = 3 10 T04 ≡ t (n): (86) % (1 − $i ni )3 43 For coherent isotropic scattering, Eqs. (44) and (22a) imply 2 c' s (0 ) E0 (0 ) 0 s (n; ) = s0 (0 ) = 2 0 = : i 2 2 0 % (1 − $i n ) 43 % (1 − $i ni )2 Integrating over frequency, and using Eq. (34a), ∞ c' E0 0 s (n; ) d = 3 (1 − $ ni )3 43 % i 0 c' [(1 + $ $ fij )E − (2$ F i =c)] i j i 0 = ≡ s (n): 43 %(1 − $i ni )3
(87)
(88)
Collecting these results in Eq. (84) one
(89)
To check the correctness of the analysis just presented, integrate the absorbtion=emission terms in Eqs. (85), (86), and (88) by taking their zeroth and
d! ni (1 − $i nj ) I (n) = 0 %[(F i =c) − $j P ij ]:
(90b)
For thermal emission, the inner product $i ni appears inside the integral. But because t (n) is a scalar, and $i ni is invariant under a rotation of the coordinate axes, one can choose the polar axis to lie along the velocity vector R. Taking ? to be the polar angle between n and R, de
94
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
To take the
d! ni s (n) = %3 '0 [(1 + $i $j fij )E − (2$i F i =c)]$i :
(92b)
Assembling Eqs. (90a), (91a), and (92a), one gets G 0 = %(0 E − 10 aR T04 ) − %0 ($i F i =c)] − %3 '0 [(1 + $i $j fij )E − (2$i F i =c)]:
(93a)
Likewise, assembling Eqs. (90b), (91b) and (92a), one
(93b)
Eqs. (93a) and (93b) are identical to Eqs. (54a) – (54d) if 0 is interpreted as 0F , 10 as 10P , and '0 as (0F − 10E ), which is correct because we have assumed the opacity is gray. Thus: (a) one can con
(94)
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
95
Doing so would account for the major e8ects of boundaries, and low opacity in some bands versus high opacity in others. 10. Overall iteration procedure Suppose the time-evolution of a radiating =ow has been followed to time-level t n , and all the hydrodynamic variables, radiation energy density and =ux, Eddington tensor elements, and radiation spectra, at all grid points, are known. Using radiation spectra and Eddington tensor elements at t n as
(i); n+1 n+1 and 0F new estimates of the radiation energy and =ux spectra, hence new values for 10E via Eqs. (66) and (67). These are used in Eqs. (54a) – (54d) to get new source=sink terms G 0; n+1 and G i; n+1 for the radiation-hydro equations, which can then be updated for a better estimate of the solution at t n+1 . In principle this process should be sub-cycled to convergence.
11. Conclusions We have outlined the major steps needed to perform a physically realistic calculation of radiation hydrodynamics using lab-frame radiation quantities. The greatest computational e8ort goes into computing the angle-frequency information needed to: (1) evaluate elements of the Eddington tensor, and (2) to evaluate energy and =ux spectra needed to calculate mean opacities source=sink terms. We argue that angular information can be generated for all angles and frequencies in parallel, neglecting velocity e8ects. In contrast, realistic radiation energy and =ux spectra require a solution of the frequency-dependent comoving-frame radiation energy and momentum equations. For this work we propose an algorithm based on opacity distribution functions, and an approximate closure of the comoving-frame radiation momentum equation. Use of lab-frame radiation variables in the radiation hydro equations makes sense only if one uses a
96
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
this approach. Use of lab-frame radiation quantities can also encounter numerical di5culties in the di8usion limit for extremely optically-thick media because the di8usive =ux gets swamped by the advective =ux, even at small values of v=c. However, we note that this di5culty may not occur in neutron-transport problems where very large scattering-thicknesses are not encountered (for non-critical systems). Comoving-frame variables are most often used in Lagrangian-mesh calculations. One then avoids interface tracking, though zones with subscale mixing still present problems. In the past, Lagrangian-mesh computations had trouble with energy and momentum conservation. But recently Caramana and his associates have developed methods [28–31] that give precise total energy conservation (computer-word accuracy), and resist deformation and=or tangling (producing “boomerang” and “bowtie”, or “hourglass”, cells), on Lagrangian meshes. The computational e8ort for a radiation hydro step is about the same whichever radiation frame is used, being dominated by construction of the radiation
D. Mihalas, L.H. Auer / Journal of Quantitative Spectroscopy & Radiative Transfer 71 (2001) 61–97
97
[10] Munier A, Weaver R. Radiation transfer in the =uid frame: a covariant formulation. Part II: the radiation transfer equation. Comput Phys Rep 1986;3:165–208. [11] Munier A, Weaver R. Radiation transfer in the =uid frame: a covariant formulation. Part I: radiation hydrodynamics. Comput Phys Rep 1986;3:125–64. [12] Dor< E, Drury L. Simple adaptive grids for 1-D initial value problems. J Comput Phys 1987;69:175–95. [13] Winkler K-H, Norman ML, Mihalas D. Implicit adaptive-grid radiation hydrodynamics. In Multiple timescales: computational techniques. vol. 2. New York: Academic Press, 1985 [Chapter 6] p. 145 –84. [14] Winkler K-H, Norman ML, Mihalas D. Adaptive-mesh radiation by hydrodynamics. I. The radiation transport equation in a completely adaptive coordinate system. JQSRT 1984;31:473–8. [15] Mihalas D, Winkler K-H, Norman ML. Adaptive-mesh radiation hydrodynamics. II. The radiation and =uid equations in relativistic =ows. JQSRT 1984;31:479–89. [16] Lowrie R, Mihalas D, Morel JE. Comoving-frame radiation transport for nonrelativistic =uid velocities. JQSRT 2001;69:291–304. [17] Synge JL. The relativistic gas. Amsterdam: North-Holland, 1957. [18] Hazlehurst J, Sargent WLW. Hydrodynamics in a radiation