Analytical Models of Exoplanetary Atmospheres. VI. Full Solutions for Improved Two-stream Radiative Transfer Including Direct Stellar Beam

Two-stream radiative transfer is a workhorse in the Earth, planetary and exoplanetary sciences due to its simplicity and ease of implementation. However, a longstanding limitation of the two-stream approximation is its inaccuracy in the presence of medium-sized or large aerosols. This limitation was lifted by the discovery of the improved two-stream technique, where the accuracy of the scattering greenhouse effect is matched to that of multi-stream calculations by construction. In the present study, we derive the full solutions for improved two-stream radiative transfer, following its introduction by Heng&Kitzmann (2017), and include contributions from the direct stellar beam. The generalization of the original two-stream flux solutions comes in the form of a correction factor, traditionally set to unity, which is the ratio of a pair of first Eddington coefficients. We derive an analytical expression for this correction factor and also provide a simple fitting function for its ease of use by other workers. We prove that the direct stellar beam is associated with a second Eddington coefficient that is on the order of unity. Setting this second Eddington coefficient to 2/3 and $1/\sqrt{3}$ reproduces the Eddington and quadrature closures, respectively, associated with the direct beam. We use our improved two-stream solutions for the fluxes to derive two-stream source function solutions for both the intensity and fluxes.


Introduction
A workhorse of radiative transfer is the two-stream method, which approximates the passage of radiation through an atmosphere as a pair of incoming and outgoing fluxes (Schuster 1905). Instead of solving a partial differential equation for the angle-dependent intensity, one solves an ordinary differential equation for the angle-integrated flux. The price to pay is mathematical underdetermination. To mathematically close the system, one has to assume that ratios of moments of the intensity are constants known as Eddington coefficients.
There is a rich literature on two-stream radiative transfer. Joseph et al. (1976) and Wiscombe (1977) introduced ad hoc modifications to the scattering phase function to treat forward-peaked scattering of radiation due to large aerosols, known as the "delta-Eddington approximation." Meador & Weaver (1980) presented a unified theoretical framework that described different ways of specifying two-stream closures. Toon et al. (1989) presented a specific formulation of the two-stream method that has been influential in the planetary and exoplanetary sciences (e.g., Marley & McKay 1999), partly because it allows for fast computation via the inversion of a tridiagonal matrix. Toon et al. (1977Toon et al. ( , 1989 introduced the "two-stream source function" method, which reduces the radiative transfer equation from an integro-differential equation to a differential equation for the intensity by utilizing a sleight: the two-stream solution is inserted into the integral involving the scattering phase function and intensity, by assuming that the intensity is related to the two-stream flux by a constant factor. Heng et al. (2014) introduced a formalism that unified the two-stream solutions with analytical solutions for temperature-pressure profiles, based partly on generalizing the relevant parts of the monograph of Pierrehumbert (2010). Kitzmann et al. (2013) realized that the original two-stream method performs poorly in the presence of medium-sized or large aerosols, as it underestimates the backscattered radiation by ∼10%. Kitzmann (2016) demonstrated that this shortcoming of the two-stream method translates into an overestimation of the scattering greenhouse effect for early Mars by about 50 K. Motivated by the work of Kitzmann et al. (2013) and Kitzmann (2016), Heng & Kitzmann (2017) discovered a simple improvement to the two-stream method that allows it to match the accuracy of 32-stream calculations (at the ∼0.01%-1% level or better, depending on the optical depth) in the presence of medium-sized or large aerosols.
The overarching intention of this study is to fully flesh out the improved two-stream method introduced by Heng & Kitzmann (2017) and also use it as input for the two-stream source function method of Toon et al. (1989). The entirety of the paper is devoted to these novel derivations. At the heart of the improved two-stream flux solutions is a correction factor E, shown in Figure 1, which is usually set to unity in the original two-stream method. We provide a convenient fitting function in Equation (31) for the reader to compute E.

Basic Theory
The radiative transfer equation may generally be written as (Chandrasekhar 1960;Mihalas 1970Mihalas , 1978Meador & Weaver 1980;Toon et al. 1989) Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
where I is the wavelength-, frequency-, or wavenumberdependent intensity 1 (which is why we have omitted subscripts on all quantities that commit us to any of these choices, including for the Planck function), τ is the optical depth, ω 0 is the single-scattering albedo, and B is the Planck function. The scattering phase function,  m f m f ¢ ¢ ( ) , , , , relates the incident polar (θ′) and azimuthal (f′) angles to the emergent polar (θ) and azimuthal (f) angles. We further define m q º cos , Physically, one has to account for the scattered intensity coming from all incident directions.
In addition to the diffuse radiation field, we assume the presence of a direct beam of radiation due to the star, which impinges on the atmosphere with μ′=−μ å and f′=f å . Some portion of this incident stellar beam is absorbed, while the rest is scattered into the diffuse radiation field, which is represented by where F å is the incident stellar flux at the top of the atmosphere, and the pair of delta functions enforces the directionality of the beam.
The radiative transfer equation is multiplied by some arbitrary function  m ( ) and integrated over all emergent angles (Pierrehumbert 2010; Heng et al. 2014;Heng 2017), Equation (3) is the starting point for the suite of derivations we will perform in the current study. The minus sign associated with the gradient terms comes about because we have rescaled the optical depth to be always positive, as has been previously elucidated in Heng et al. (2014) and Chapter 3.5.1 of Heng (2017).
Mathematically, the contribution of the direct beam to  may be evaluated in two ways: either by evaluating the delta functions first (over dΩ′) or by evaluating  first. Since these mathematical operations commute, the resulting expressions must agree. The former approach corresponds to the derivation for isotropic scattering, while the latter is for non-isotropic (or anisotropic) scattering.

Improved Two-stream with Direct Beam and Isotropic Scattering
In the isotropic limit, Equation (1) becomes The total intensity is and it is related to the total flux by ò=F + /J. The ratio of the flux to the total intensity in one hemisphere only is denoted by ò′, where we have assumed that this ratio is the same in both hemispheres.
By integrating over the incoming () and outgoing () hemispheres, we obtain the incoming and outgoing fluxes, 2 The coefficients of the various terms are given by The derivation technique was established in Heng et al. (2014) and also Chapter 3 of Heng (2017), and we will not repeat it here. Instead, we will point out key, novel aspects of the current derivation that involve both the improvement to the two-stream method and the inclusion of the direct beam. Previously, Heng (2017) designated both ò and ò′ as the "first Eddington coefficients." By demanding that the limiting case of an opaque, purely absorbing (ω 0 =0) atmosphere emits π B of flux in each hemisphere, one obtains ò=1/2 (Toon et al. 1989;Pierrehumbert 2010). By demanding that a purely scattering (ω 0 =1) atmosphere remains in radiative equilibrium, one sets ò′=ò (Toon et al. 1989). However, since the two-stream solutions are formally invalid in the limit of ω 0 =1, the second demand is academic, as has already been realized by Heng & Kitzmann (2017). Following Heng & Kitzmann (2017) and proceed to obtain an analytical expression for E in terms of ω 0 and the scattering asymmetry factor (g 0 ).
We have written Since the integration involves only μ and f, we expect that  has the same mathematical behavior as   with respect to this procedure. Physically, for isotropic scattering we expect c c = =   1 2. We keep these quantities distinct, because we will explicitly prove in Section 4 that c c =   is a natural outcome of setting g 0 =0 for the more general expressions involving the gradients of the fluxes. When scattering is non-isotropic, we generally expect that c c ¹   , but because the scattering phase function needs to be properly normalized, we have In both Meador & Weaver (1980) and Toon et al. (1989), these authors have written g c We will defer seeking the solutions to F + and F − , and hence  F and  F , until we deal with the non-isotropic treatment, because the mathematical forms of these governing equations are almost identical.

Governing Equation for Net Flux
For non-isotropic scattering, we return to using Equation (3) as our starting point. If we set  = 1, then we obtain  = 1 (Pierrehumbert 2010; Heng et al. 2014;Heng 2017). The equation governing the net flux becomes There is an unsatisfactory ambiguity if   ¹ ¢, because one could convert J to the total flux using . We resolve this ambiguity by demanding that the equation reduces to its counterpart in the isotropic limit, which yields 14 0 0 0

Governing Equation for Total Flux
To obtain the equation governing the total flux, we set  m = , which yields  m = ¢ g 0 (Pierrehumbert 2010; Heng et al. 2014;Heng 2017), from which it follows that K K K and  K and  K are the second moments of the intensity. Following Heng et al. (2014) and Heng (2017), we define a second Eddington coefficient, For the diffuse emission (which is usually dominant in the infrared range of wavelengths), we demand that  = + -F EF 2 2 , in order for the governing equation for the total flux to reduce to its counterpart in the limit of isotropic scattering. The governing equation becomes

Coefficients in Governing Equations
For the direct beam (which is usually dominant in the optical/visible range of wavelengths), we intuitively expect that ~1 2 . We will now demonstrate this claim. By comparing Equations (14) and (17) with the pair of equations in (12), we may infer that g g w g g w In the limit of g 0 =0 (isotropic scattering), the beam term in Equation (17) vanishes. By comparing to the first equation in (12), we obtain c c =   . Another useful comparison is between Equations (14) and (17) and the pair of equations in (7), which yields When we compare c  with the g 3 values listed in Table 1 of both Meador & Weaver (1980) and Toon et al. (1989), we find that  = 2 3 2 reproduces the Eddington closure, while  = » 1 3 0.58 2 reproduces the quadrature closure. Regardless of the choice of value for the second Eddington coefficient, we always have c c + =   1, which ensures that the scattering phase function is normalized properly.
The physical interpretation of the second Eddington coefficient is that, in the limit of   m = 1 2 , an extra fraction g 2 0 of the stellar direct beam gets scattered into the forward/ incoming direction, such that when

Solutions for the Two-stream Fluxes
The recipe to solve for the fluxes has previously been established by both Heng et al. (2014) and Chapter 3 of Heng (2017). We will only highlight the key points and novel generalizations. One first merges Equations (14) and (17) into a single second-order ordinary differential equation for + F . The solution is s a s and the coefficients A 1 and A 2 are determined by applying boundary conditions to the solution. A linear expansion in the Planck function is performed to allow for the two-stream solutions to treat nonisothermal atmospheric layers described by a gradient t ¢ º ¶ ¶ B B , which has been shown to significantly increase computational efficiency (Malik et al. 2017). The index i=1, 2 will be explained shortly. Note the correction factor in front of the Planck function associated with ¹ E 1. The coefficient  C is derived using the method of undetermined coefficients, From the expressions for  F and  F , one realizes that a more compact way of expressing the algebra is to define the coupling coefficients, These are termed coupling coefficients, because they couple the outgoing and incoming fluxes in the presence of scattering (see Section 3.2 of Heng 2017). As already described in Toon et al. (1989),  C becomes indeterminate when the denominator in Equation (23) becomes zero, which occurs when Toon et al. (1989) suggest that, "In practice, if the equality happens to occur, this problem can be eliminated by simply choosing a slightly different value of  m ." The final two-stream solutions for the incoming and outgoing fluxes are expressed in the convention of a pair of atmospheric layers, subscripted by "1" and "2" with the former residing above the latter. In this convention,  F 2 and  F 1 are the boundary conditions, which are used to eliminate A 1 and A 2 . This procedure is documented in Section 3.8.2 of Heng (2017).
When the dust has settled on the algebra, the two-stream solutions are  where, for compactness of notation, we have defined  The optical depth of the atmospheric layer is given by Δ. The transmissivity of the layer is represented by  . Extinction of the direct stellar beam above and through the layer are given by  above and   , respectively. It is instructive to examine the limit of = g 1 0 and E=1, because in the two-stream approximation it behaves as if one is in the limit of pure absorption. In this limit, the coupling coefficients become z = -0 and z = + 1, implying that several of the direct beam terms in Equations (25) and (26) . One may verify that regardless of whether is always positive.

The E-factor
The ratio of the first Eddington coefficients, E, lies at the heart of our improved two-stream method. In the optically thick limit, the fraction of flux reflected by an atmospheric layer is z z -+ , a fact that may be verified by setting  = 0 in Equations (25) and (26) and ignoring the blackbody and direct beam terms. The expression for z z = ¥ -+ R may be manipulated to obtain an analytical expression for E. In other words, E is constructed to reproduce the reflectivity in the optically thick limit ( ¥ R ), as long as one has a way to numerically evaluate ¥ R . Previously, Heng & Kitzmann (2017)  and reasoned that it generalized to g w = - With the formal, complete derivation performed in the present study, we obtain a different expression for E, 1 . The two-stream fluxes depend on w 0 and g 0 only via the coupling coefficients zand z + . When one inserts either expression for E into its corresponding expression for z  , one obtains an identical result, In other words, E is constructed such that the coupling coefficients only depend on ¥ R . Similarly, because of the ambiguity with relating J,  J , and  J to the fluxes via ò and  ¢, there are multiple ways of expressing the coupling coefficients in terms of E. Nevertheless, they all reduce to the same expression for z  as stated above. This property implies that the results shown in Figures2, 3 and 4 of Heng & Kitzmann (2017) are identical 3 even when the expression for E in Equation (29) is used, which is why we have not reproduced them in the current study. The conclusions of Heng & Kitzmann (2017) remain qualitatively and quantitatively unchanged. Figure 1 shows the calculations of E using Equation (29).
¥ R is numerically evaluated using a brute-force, 32-stream discrete ordinates method (Stamnes et al. 1988) as implemented by the DISORT code (Hamre et al. 2013). In these numerical calculations, one has to specify an explicit form of the scattering phase function. A Henyey-Greenstein scattering phase function (Henyey & Greenstein 1941) is used when , the error may be as high as 2%.

Improved Two-stream Source Function Method
Starting from Equation (1), it is always possible to rewrite the radiative transfer equation as The trick employed by the two-stream source function method of Toon et al. (1989) is to argue that µ   I F F , in the integral. The scattering phase function is assumed to be p in the forward and reverse directions, respectively. Since the two-stream fluxes do not depend on W¢, the integral becomes trivial to evaluate. The direct beam term is ignored, because the two-stream source function method was originally designed to treat diffuse thermal emission.
If we focus only on the integral term involving the diffuse emission, where we need p =  I F in order for the preceding expression to reduce correctly to its isotropic limit of w p J 4 0 when = g 0 0 . This is a departure from Toon et al. (1989), who assumed  =  I F . The minus and plus signs depend on whether one is writing the integral term for the outgoing () or incoming () flux.
We proceed somewhat differently from previous studies. In Toon et al. (1977), who first introduced the two-stream source function approximation, the form of their governing equations and the absence of g 0 in them suggest that their formalism does not generally treat ¹ g 0 0 scenarios. In Toon et al. (1989), who generalized their formalism to also consider ¹ g 0 0 values, their Equations (53) and (54) imply that their two-stream fluxes may be written as exponentials and linear terms involving the optical depth. The coefficients of these terms do not appear to depend on the optical depth. This is contrary to what we find for our two-stream fluxes in Equations (25) and (26), where the denominator contains the transmission function.
In our approach, we consider the two-stream fluxes to already have undergone the integration from t 1 to t 2 , such that Note that  F 1 and  F 2 are our two-stream fluxes given by Equations (25) and (26), respectively.
If we again define t t D º -  where we have used the symbols  F 1 and  F 2 to distinguish the two-stream-source function fluxes from the two-stream fluxes. We have defined where E 3 is the exponential integral of the third order (Abramowitz & Stegun 1970;Arfken & Weber 1995). The reversal in the minus and plus signs for the two-streamsource function fluxes accounts for the change in the forward versus reverse directions for the incoming versus outgoing fluxes. For example, when = g 1 0 , there should be no correction due to backscattering and all of the scattered flux should be in the forward direction.

Summary and User's Manual
The computational implementation of the two-stream solutions of radiative transfer has previously been described elsewhere (e.g., Malik et al. 2017). Instead, we focus on summarizing the different variants of the two-stream method from a practical standpoint.
1. If the reader is interested in modeling only diffuse infrared emission in the absence of medium-sized or large aerosols, the original two-stream formulation may suffice. For a concise statement of the two-stream fluxes with the hemispheric closure, see Equation (3.58) of Heng (2017). Alternatively, use Equations (25) and (26) of the present study with E=1 and discard/ignore the direct beam terms. 2. If the reader is interested in including the direct stellar beam in the original two-stream solutions, see Meador & Weaver (1980) and Toon et al. (1989). 3. If the reader is interested in accurately modeling diffuse infrared emission in the presence of medium-sized or large aerosols, then our improved two-stream formulation in Equations (25) and (26) should be used. Consult Table1 of Meador & Weaver (1980) or Toon et al. (1989) for the various choices of closure associated with the direct stellar beam (i.e., value of  2 ). The E-factor in Equation (29) should be used in tandem with the coupling coefficients in Equation (24). For convenience, we provide a simple fitting function for E in Equation (31). 4. Our improved two-stream solutions may be used to derive the two-stream source function solutions, which we present in Equations (36) and (37) for the intensity and fluxes, respectively.
For completeness, the summary of key symbols used and their correspondence to the symbols used in Toon et al. (1989) are listed in Table 1. Appendix lists some typographical errors we found in Meador & Weaver (1980) and Toon et al. (1989). Typographical Errors in Meador & Weaver (1980) and Toon et al. (1989) A.1. Toon et al. (1989) We use the same notation as Toon et al. (1989) and m f ( ) , , for the incident and emergent directions, respectively. Equation (1) of Toon et al. (1989) concisely states the functional dependence of their scattering phase function: . The scattering phase function associated with the direct beam is stated in their Equation (3): , , 0 0 . Based on their stated convention and also their Equation (10), we may conclude that Equation (9)  Specifically, the integration should be over μ and not m¢.
In Equation (12), there is a typographical error: + S n should instead be -S n . The logical flow of the text into Equations (13) and (14) of Toon et al. (1989) makes this clear.
A.2. Meador & Weaver (1980) Our notations for μ, m¢, f and f¢ are the same as those of Meador & Weaver (1980). Equation(1) of Meador & Weaver (1980) concisely states the functional dependence of their scattering phase function: m f m f ¢ ¢ n ( ) P , ; , . It also states the functional dependence of the scattering phase function associated with the direct beam: m f m f n ( ) P , ; , 0 0 . Equation(7) of Meador & Weaver (1980) follows logically from their Equation(3). However, their Equation(9) suffers from the same logical inconsistency as what we pointed out for Toon et al. (1989), which is that the integration should be over μ and not m¢. This is clearly seen in the transition from Equations (7) to (9). It is not an issue of a mere change of notation for the integration variable, because it needs to correspond to the correct variable in the scattering phase function. Equation It is likely that this error in Meador & Weaver (1980) propagated into Toon et al. (1989).