Radiative transfer in protoplanetary disks under irradiation by the protostar

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Radiative transfer in a protoplanetary accretion disk, which is irradiated by the protostar, is examined under the non-gray frequency-dependent treatment, including the anisotropic scattering effect. Due to irradiation heating, the disk is divided into an inner optically thick flat disk, a middle optically thick flared disk, and an outer optically thin flared disk. In each regime, emergent intensities as well as other radiative quantities are analytically obtained under the Eddington approximation. In the inner flat disk, where the disk is optically thick and viscous heating is dominant, we obtain solutions for the case of vertically uniform heating. In the middle flared disk, where the disk is optically thick and geometrically flared due to irradiation heating, we obtain solutions for the local thermodynamic equilibrium (LTE) case under the assumption of a linear Planck function. In the outer flared disk, where the disk is optically thin and flared due to irradiation heating, we obtain solutions for the LTE case of a vertically isothermal disk. The temperature structure is consistent with that seen in previous studies, while disk intensities are different from those in previous cases without anisotropic scattering. We found that anisotropic scattering enhances the emergent intensity in the poleward direction, when forward scattering dominates. As a result, in addition to the usual limb-darkening effect, the protoplanetary disk also becomes bright in the poleward direction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .


Introduction
Accretion disks are believed to be luminous energy sources in the universe, i.e., protostars and protoplanetary disks, cataclysmic variables and supersoft X-ray sources including white dwarfs, microquasars, and galactic X-ray binaries including neutron stars and black holes, ultraluminous X-ray sources, quasars, and active galaxies, and so on.Accretion-disk models have been extensively studied during the past four decades; many researchers have investigated their internal structure, observational spectra, stability, formation and evolution, and oscillations (see Refs. [22,23] for reviews).
Accretion disks around a protostar-protoplanetary disks-are different from those around compact stars in two points; irradiation of the central star is usually important, and scattering of dust grains is also important, since dust is the main absorber of stellar radiation in protoplanetary disks.
Hence, the study of protoplanetary accretion disks started, with particular attention to the irradiation effect of the central protostar [8,9,14,16,24,25].In some cases irradiation heating dominates over viscous heating; such an irradiated disk is often called a passive disk, in contrast to an active disk, where viscous heating is dominant [22,23,32].In these classical studies, in an irradiated passive PTEP 2013, 053E02 J. Fukue disk with an infinitesimal thickness, the surface temperature T varies as T ∝ r −3/4 for large radius r , whereas, in an irradiated disk with a flaring thickness, T ∝ r −3/7 for large r .
The effect of irradiation on the disk structure has been extensively studied in recent years [4,5,29].In these current studies, however, the attention of researchers has focused mainly on the irradiation effect on the vertical and radial structures of the irradiated disk and the resultant spectra.
For example, D'Alessio et al. [6,7] studied the vertical structure of irradiated protoplanetary disks, and obtained detailed vertical and radial structures.However, they treated the radiative transfer only approximately.That is, they separately calculated the thermal, scattered, and incident components, and adopted semi-infinite solutions, although the disk's total optical depth is finite.
Jang-Condell and Sasselov [20] and Jang-Condell [18,19] considered radiative transfer in protoplanetary disks, related to the planet shadows, under a gray approximation.However, their boundary conditions for the transfer problem were not appropriate for those under irradiation, and they did not give analytical solutions for, e.g., emergent intensities.That is, Jang-Condell and Sasselov [20] adopted the boundary condition for a Lambertian surface without irradiation, which is inconsistent with the irradiated disk (see Eq. (44) below).Furthermore, they separately calculated the thermal and scattered components.
While Jang-Condell and Turner [21] gave several analytical expressions, they adopted the solutions of Jang-Condell and Sasselov [20], and, therefore, their boundary conditions would also not be appropriate.Furthermore, they also treated the scattered and thermal lights separately.
Thus, analytical studies of radiative transfer in protoplanetary disks including scattering are few, and often the boundary conditions are not appropriate.Furthermore, the anisotropic properties of scattering are not examined in the current analytical studies.
On the other hand, the radiative transfer problem in accretion disks around compact objects under irradiation has recently been examined [10,11,13].For example, Fukue [11] examined radiative transfer in an isothermal disk with finite optical depth under irradiation by the central compact star.However, the situations are rather different in protoplanetary disks; e.g., the scattering process is electron scattering in accretion disks around compact objects, while in protoplanetary disks it is Mie scattering [6].
Thus, in this paper we analytically re-examine radiative transfer in protoplanetary disks irradiated by the central protostar using a frequency-dependent treatment with adequate boundary conditions, focusing our attention on the (anisotropic) scattering process [12].Although the temperature structure obtained in the present study is similar to that of the previous studies, the disk intensity could be different from that in the previous cases due to the anisotropic scattering effect.
In the next section we describe the model and basic equations.We obtain the vertical and radial structures of the present analytical model in Sects.3 and 4, respectively.We then derive analytical solutions of the radiative transfer problem for the frequency-dependent case in Sect. 5.The final section is devoted to concluding remarks.

Basic equations
In this paper we assume the following situations: (i) The disk is steady and axisymmetric.(ii) It is also geometrically thin and locally plane-parallel, although there exists a flaring effect.(iii) The viscous PTEP 2013, 053E02 J. Fukue heating per mass is uniform in the vertical direction, if it is dominant over the irradiation heating.(iv) The radiation field in the disk is sufficiently isotropic, and we adopt the Eddington approximation.(v) The disk atmosphere is under local thermodynamic equilibrium (LTE) or vertically uniform heating.Under these situations, we analytically consider external irradiation and radiative transfer in the irradiated disk with the anisotropic scattering effect (see Refs. [13] for gray treatment; [10] for scattering; and [11] for irradiation).

For radiation
The radiative transfer equations are given in several pieces of work [3,[10][11][12]23,27,28,31].For planeparallel geometry in the vertical direction (z), the frequency-dependent transfer equation, the zeroth moment equation, the first moment equation, and the Eddington approximation become, respectively, where μ is the direction cosine (= cos θ ), I ν the specific intensity, J ν the mean intensity (= cE ν /4π , E ν being the radiation energy density), H ν the Eddington flux (H ν = F ν /4π , F ν the vertical component of the radiative flux), K ν the mean radiation stress (K ν = cP ν /4π , P ν the zz-component of the radiation stress tensor), ρ the gas density, and c the speed of light.The mass emissivity j ν , the absorption opacity κ ν , and the scattering opacity σ ν generally depend on the frequency.
Here, we approximately take into account anisotropy in scattering [12], bearing in mind the dust opacity; i.e., the parameter a in these equations means the fraction of anisotropy from the isotropic case (a = 0) to the extremely anisotropic one (a = ∞).These equations including the anisotropic scattering terms were derived from a model phase function [12], which approximately expresses the Henyey-Greenstein phase function.
Introducing the frequency-dependent optical depth, defined by dτ ν ≡ −ρ (κ ν + σ ν ) dz, we rewrite the radiative transfer equations in the familiar frequency-dependent form: where ε ν is the photon destruction probability: 3/17 and δ ν the anisotropic parameter: Finally, the disk's total optical depth under the frequency-dependent form becomes where H is the disk's half-thickness.Integrating Eqs. ( 1)-( 3) over the frequency, we obtain frequency-integrated equations for radiation: where j ≡ j ν dν, [23].Introducing the frequency-integrated optical depth, defined by dτ ≡ −ρ (κ I + σ I ) dz, we rewrite the radiative transfer equations as where then Eqs. ( 15)-( 17) become, respectively, where 4/17 It should be noted that this anisotropic parameter is generally a function of z or τ .The physical properties of dust particles and the size distribution of particles change in the vertical direction, and hence, the photon destruction probability also varies in the vertical direction.In the present model, however, we assume that this anisotropic parameter is constant in the vertical direction.
Using the energy conservation between radiation and matter, where q + vis is the internal viscous heating rate, these frequency-integrated equations are finally expressed as where is the viscous heating rate in unit mass, and is assumed to be constant in this paper (uniform heating).
The following should be noted on this assumption of uniform heating.In a realistic accretion disk, the viscous heating rate would be a function of the optical depth.For example, in the recent 3D radiation magnetohydrodynamics (MHD) simulations using the shearing-box approximation [15], the disk interior is colder than the classical viscous model, while the disk atmosphere is ionized by stellar X-rays and is sufficiently turbulent for the dissipation per unit mass to be greatest.However, in the present paper we assume uniform heating per unit mass for simplicity, and concentrate on the effect of anisotropic scattering.

For matter
The basic equations of standard accretion disks around a central object of mass M are given in, e.g., Ref. [23].In protoplanetary accretion disks, the gas pressure dominates over the radiation pressure, and the opacity is mainly that of dust.The disk is assumed to be steady and axisymmetric.The physical quantities are vertically averaged.For simplicity, we drop the boundary factor.
The basic equations for matter are then: 5/17 where (= 2ρ H ) is the surface density, ρ the density, H the disk thickness, Ṁ the mass-accretion rate, v the radial velocity, v ϕ the azimuthal velocity, T r ϕ the r ϕ-component of the stress tensor, p the gas pressure, T c the central temperature at the disk equator, T eff the effective temperature at the disk surface, Q vis the viscous heating rate, Q irr the irradiation heating rate, and α the alpha parameter.
It should be noted that the effective temperature is meaningful when the disk is optically thick in the vertical direction.When the disk is optically thin, on the other hand, we assume that the disk is isothermal in the vertical direction, as below.
For the present purpose, after some manipulations, we have equations for , H , T c , T eff , and τ b : where κ and σ are the mean absorption and scattering opacities, respectively; both are a few cm 2 g −1 for dust in protoplanetary disks at temperatures of 200-1200 K [2,7,19].In the outer disk, where the temperature is much lower, thermal emission occurs at long wavelengths and the dust is more transparent, while, in the inner disk, the temperature is so high that silicate grains evaporate and the mean opacity should be reduced.In the following, therefore, the mean opacities are left as a parameterized form.
The viscous heating rate and that of irradiation by the central star (with radius R * , temperature T * , luminosity L * ) are given, respectively, by where we also drop the boundary factor in Eq. ( 40).In the irradiation heating rate, the first term comes from the size effect of the central star, which is important in the inner region, while the second term is the flaring effect of the disk thickness, which is important in the outer region [9].When H ∝ r h+1 , (d H/dr − H/r ) = h H/r .

Vertical disk structure
In general, there is internal heating in the viscous accretion disk.If the viscous heating is concentrated at the disk equator, then the vertical flux in the disk atmosphere may be constant [13].If, on the other hand, the viscous heating is distributed in the disk, then the vertical flux depends on the optical depth [1,10,13,16,17,26].In this paper we consider an irradiated accretion disk with internal viscous heating, which is assumed to be uniform as stated in Eq. (29).

Radiative moments
When the viscous heating rate q is constant, Eqs. ( 27) and ( 28) are easily integrated to give where H (0) and J (0) are the surface values of H and J , respectively.Now, we impose the boundary conditions.At the disk equator (τ = τ b ), the flux is zero: therefore, q = H (0)/τ b .At the disk surface (τ = 0), on the other hand, we should impose the boundary condition under irradiation: where A is the gray albedo, W the total dilution factor defined as W J and W H being the dilution factors for J and H , respectively, and I * (= σ T 4 * /π ) the irradiation intensity [11].Here, the dilution factors W J and W H are, respectively, calculated as where It should be stressed that the boundary conditions used in the previous studies were not adequate [18][19][20][21].For example, the boundary condition imposed in Ref. [20] was J (0) − 2H (0) = 0 (Eq.( 13) in Ref. [20]), which is the simple Lambertian condition and does not include the irradiation effect as in Eq. (45).Jang-Condell and Turner [21] used the solutions of Jang-Condell and Sasselov [20] for scattered light, and, therefore, their boundary conditions do not include the irradiation effect.Although we assume a Lambertian surface for simplicity, the boundary condition (45) takes the irradiation effect appropriately into consideration.
Under these boundary conditions, the solutions of the radiative moments become We then express W and H (0) by system parameters such as R * and T * .As derived in Ref. [11], for r R * , the total dilution factor is approximately expressed as Hence, PTEP 2013, 053E02 J. Fukue On the other hand, from the heating rates (40) and ( 41), H (0) is expressed as

Temperatures
Here we summarize several expressions for temperatures in and on the irradiated accretion disk.At first, from Eq. ( 54), the effective temperature T eff (0) at the surface is This effective temperature is often used in the traditional irradiated-disk model.If we assume that the radiation energy density in the disk is expressed by a blackbody one, from Eq. ( 51), the temperature distribution T (τ ) in the disk becomes Hence, the surface temperature T (0) is expressed as In Eqs. ( 55) and (57), the first, second, and third terms on the right-hand side of Eq. ( 55) correspond to the second, third, and fourth terms on the right-hand side of Eq. (57), respectively, beside the factor 2. These terms just satisfy the usual Milne-Eddington atmosphere, where the relation between the surface temperature T (0) and the effective temperature T eff is T 4 (0) = 1 2 T 4 eff .On the other hand, the first term on the right-hand side of Eq. ( 57) comes from the irradiated boundary condition, which does not appear in the usual Milne-Eddington atmosphere.
Finally, from Eq. ( 56), the central temperature T c at the equator is expressed as 8/17 and the stellar luminosity L * (= 4π R 2 * σ T 4 * ) these temperatures are finally expressed as where we introduce the normalized radius r = r/R * .

Radial disk structure
Using the temperature distribution and equations of matter, we can obtain solutions of radial disk structure analytically using several combinations of optical depth and disk flaring.In what follows, we adopt the standard parameters for protoplanetary disks [18,19]: We first assume that the disk is optically thin in the vertical direction; τ b < 1.In this case, we assume that the disk is isothermal in the vertical direction, and do not distinguish the effective temperature and the internal temperature; T c = T (0) ∼ T eff (0) = T , and from Eq. (37) we can obtain analytical solutions for T and H as usual.We further divide the disk into the inner flat disk, where the first term of the right-hand side of Eq. ( 55) is dominant, and the outer flared disk, where the second term is dominant.We then obtain analytical solutions for the inner flat disk and outer flared one, but the inner flat disk and some part of the outer flared one turn out to be optically thick.Only outside of r ∼ 10 4 (∼ 100 AU) are the obtained solutions optically thin.Hence, for an outer flared optically thin disk of r > 10 4 , we have = 125 400 g cm We next assume that the disk is optically thick in the vertical direction; τ b > 1 (rigorously speaking, ( Ĥ /r )δτ b > 1).In this case, Eq. (62) becomes and with the help of Eqs.(36), (37), and (39) we can obtain analytical solutions for T c , H , and .We further divide the disk into the inner flat disk and the outer flared disk.We then obtain analytical solutions for the inner flat disk and outer flared one, but some part of the outer flared disk turns out to be optically thin.Inside of r ∼ 10 4 (∼ 100 AU) of the outer flared disk, the obtained solutions are optically thick.Hence, for the middle flared optically thick disk of 70 (∼ 1 AU) < r < 10 4 , we have where h = 1/6.In addition, the inner flat disk of r < 70 (∼ 1 AU) is also optically thick.For this inner flat optically thick disk of r < 70, we have The boundary between the inner flat and middle flared disks (r ∼ 70 R * ∼ 1 AU) is consistent with the result in Ref. [7], where the boundary is about 2 AU.In addition, they also obtained similar PTEP 2013, 053E02 J. Fukue results for the temperature distribution, although they did not give the analytical expressions.That is, the temperature distribution and the radial disk structure are essentially the same as those in the previous studies.In Fig. 1 the temperature distributions in the present study are shown [6,7].
It should be noted that in the outer disk the temperature is so low that the thermal emission occurs at long wavelengths and the dust is more transparent, as already stated.Even in such a case, the optically thin state is still satisfied.On the other hand, in the innermost region the temperature becomes so high that the silicate grains evaporate and the mean opacity is reduced.In such a case the optically thick condition in the innermost region would be violated.
In addition, in this paper we assume that dust grains are uniformly distributed in the disk gas for simplicity.As the protoplanetary disks evolve, however, the dust grains settle into the interior of the disk to form planetesimals.As a result, the dust abundance decreases, and the disk opacities would also be reduced.Such a case is beyond the scope of the present paper.

Radiative transfer
Now, we examine the radiative transfer problem of protoplanetary accretion disks in each region; the inner flat optically thick disk, the middle flared optically thick disk, and the outer flared optically thin disk.We assume that the disk is geometrically thin even in the outer flared part, and use the plane-parallel approximation.
It should be noted that the present formalism including anisotropic scattering assumes an axisymmetry on the surface, while the incident stellar radiation field is not axisymmetric, and depends on the azimuthal angle.Thermal and scattered light is not influenced by this non-axisymmetric effect, since it is determined by the diluted radiation field, while direct light should retain its original direction.

Inner flat optically thick disk
In the inner flat optically thick disk, viscous heating dominates over irradiated heating.Hence, we consider internal heating, similar to the frequency-integrated case.That is, in any frequency we assume vertically uniform heating: PTEP 2013, 053E02 J. Fukue As already stated, the incident stellar radiation depends on the azimuthal angle as well as the polar angle.Hence, in the boundary condition (116) the angle should be replaced by more general angles.This boundary condition determines the direct light, the first terms on the right-hand sides of Eqs. ( 117) and (118).In other words, these first terms only exist in the stellar directions.
It should be mentioned that the consistency of the Lambertian surface for the disk under irradiation and anisotropic scattering.As is well known, the Lambertian surface assumes that the radiation emits isotropically, but does not make any assumption on the scattering process.In classical cases without irradiation and anisotropic scattering, after solving the transfer equation, the surface radiation is no longer isotropic, but often becomes slightly anisotropic, like so-called limb-darkening effects.Hence, the obtained solutions are slightly inconsistent with the original condition.However, as the lowest order approximation, the Lambertian surface is often assumed.In the present case with irradiation and anisotropic scattering, the situation is similar; the obtained solutions of the transfer equation slightly violate the boundary condition imposed.However, in the present study, we adopted the Lambertian surface at the lowest order approximation, as in the traditional studies.

Concluding remarks
In this paper we have analytically solved the radiative transfer problem of a geometrically thin protoplanetary disk, which is irradiated by the central protostar, and obtained analytical expressions for the emergent intensity as well as other radiative quantities.In particular, we solved the frequency-integrated transfer equations consistently with the hydrodynamical ones, and the frequency-dependent transfer equations under the appropriate boundary conditions for irradiation.We also focused our attention on anisotropic scattering.
The main results are as follows.First, we solved the frequency-integrated transfer equations, and then we derived analytical expressions for several temperature distributions (55)-(58) under the irradiated boundary conditions.Second, we also showed that the radial disk structure is divided into three parts analytically.Third, we next solved the frequency-dependent transfer equations, and obtained analytical expressions for the emergent intensity and moment quantities, including both irradiated and reprocessed light, for the inner flat optically thick disk, the middle flared optically thick one, and the outer flared optically thin one.Fourth, for the anisotropic scattering effect, we confirmed the results of the previous study [12], that it enhances the intensity and other quantities in the poleward (upward) direction, when forward scattering dominates.These properties are common to all three regimes.Of these, the first and second conclusions are not new but essentially the same as those obtained in the previous studies.On the other hand, the third and fourth results are new, as we are the first to solve the transfer equation including anisotropic scattering.
In the present study, we have not applied the present analytical solutions to observational results such as spectral energy distributions (SEDs), since there are many other complicated factors in protoplanetary disks.Instead, the present analytical expressions for various quantities will help us to interpret the observational data as well as numerical simulations including many complicated ingredients.

Fig. 1 .
Fig. 1.Radial distribution of the midplane T c and surface T (0) temperatures in the typical cases.Thick curves represent the midplane temperatures T c of the inner flat disk (solid curve), the middle flared disk (dashed), and the outer flared disk (dotted).Thin curves denote the surface temperature T (0) of the inner disk (solid curve), and the middle disk (dashed).
Downloaded from https://academic.oup.com/ptep/article-abstract/2013/5/053E02/1511133 by guest on 10 December 2018 Downloaded from https://academic.oup.com/ptep/article-abstract/2013/5/053E02/1511133 by guest on 10 December 2018 As already stated, when H ∝ r h+1 , d H/dr − H/r = h H/r .In this case, the first term on the righthand side of Eq. (58) comes from the irradiated boundary condition, but other terms are the usual viscous and irradiated heating cases.Using the accretion luminosity