Spherical Harmonics for the 1D Radiative Transfer Equation. II. Thermal Emission

Approximate methods for radiative transfer equations that are fast, reliable, and accurate are essential for the understanding of atmospheres of exoplanets and brown dwarfs. The simplest and most popular choice is the “two-stream method,” which is often used to produce simple yet effective models for radiative transfer in scattering and absorbing media. Toon et al. (hereafter, Toon89) outlined a two-stream method for computing reflected light and thermal spectra that was later implemented in the open-source radiative transfer model PICASO. In Part I of this series, we developed an analytical spherical harmonics method for solving the radiative transfer equation for reflected solar radiation that was implemented in PICASO to increase the accuracy of the code by offering a higher-order approximation. This work is an extension of this spherical harmonics derivation, to study thermal emission spectroscopy. We highlight the model differences in the approach for thermal emission and benchmark the four-term method (SH4) against Toon89 and a high-stream discrete-ordinates method, CDISORT. By comparing the spectra produced by each model, we demonstrate that the SH4 method provides a significant increase in accuracy, compared to Toon89, which can be attributed to the increased order of approximation and to the choice of phase function. We also explore the trade-off between computational time and model accuracy. We find that our four-term method is twice as slow as our two-term method, but is up to five times more accurate, when compared with CDISORT. Therefore, SH4 provides excellent improvement in model accuracy with minimal sacrifice in numerical expense.


INTRODUCTION
Studying the atmospheres of planets and substellar objects relies on computationally efficient methods to solve the radiative transfer equation in scattering and absorbing media.However, exact solutions typically do not exist.Scientists rely on approximate, parameterized methods to estimate solutions, but inclusion of intricate microphysical detail can render these models computationally intractable for useful applications (Stephens & Preisendorfer 1984;Thomas & Stamnes 2002;Chandrasekhar 1960;Liou 2002).The goal of radiative transfer parameterization in numerical models for exoplanet and brown dwarf atmospheres is to provide computationally efficient yet accurate methods to calculate radiative fluxes and heating rates (Stephens 1984).The approach to solving the radiative transfer equation is frequently chosen through a balance between accuracy and computational efficiency.In many practical cases, there are significant uncertainties associated with defining characteristics of the atmosphere, such as the composition, scattering phase function and opacities.Such uncertainties often dominate over small model errors in the solution, and therefore, obtaining a computationally efficient solution often takes precedence.
The most popular approximate methods for solving the radiative transfer equation are the (1) discrete-ordinates method (Chandrasekhar 1960;Stamnes et al. 1988Stamnes et al. , 2000)), (2) Monte-Carlo method (Modest 2013;Iwabuchi 2006) and (3) spherical harmonics method (Modest 1989(Modest , 2013;;Olfe 1967;van Wijngaarden & Happer 2022).The general approach of the discrete-ordinates method (DOM) is to discretize the solid angle by a finite number L of directions or "streams", along which the radiative intensities are tracked.DISORT (Stamnes et al. 1988(Stamnes et al. , 2000) ) is an example of a discrete ordinate algorithm for radiative transfer that is capable of simulating thermal emission, absorption, and scattering for arbitrary phase functions across the electromagnetic spectrum.The convergence of DOM can depreciate for optically thick media (Modest 2013;Fiveland & Jessee 1996;Lewis & Miller 1984), however, there exist a number of acceleration schemes to improve the convergence rate of DOM (Fiveland & Jessee 1996;Lewis & Miller 1984).Monte-Carlo methods track emitted photons throughout the media, and although accepted to be a largely accurate method, it is computationally taxing which makes it unsuitable for some applications (Iwabuchi 2006;Mayer 2009).The spherical harmonics (SH) approximation, denoted P L−1 , operates by expanding the intensity and phase function into a series of L spherical harmonics, or Legendre polynomials.This decouples spatial and directional dependencies.This method involves fewer equations than DOM and is potentially more accurate with comparable computational expense, but higher order expansions are mathematically complex and increasingly difficult to implement as L increases (Ge et al. 2015;van Wijngaarden & Happer 2022).
Such models have been studied for reflected solar radiation (e.g., most recently in our Part 1. Rooney et al. 2023).For example, the two-stream discrete-ordinates (L = 2) and two-term spherical harmonics (P 1 ) techniques are often used to provide simple yet effective models for atmospheric radiative transfer and are widely considered some of the simplest and most prolific approximations (Meador & Weaver 1980;King 1986;Chandrasekhar 1960;Mihalas & Mihalas 2013;van Wijngaarden & Happer 2022;Li & Ramaswamy 1996;Zhang & Li 2013).Two-stream methods are most useful for obtaining angle-averaged quantities such as heating rates and albedos (Schuster 1905;Meador & Weaver 1980;Heng & Marley 2018).
These studies of reflected solar radiation have shown that even though the two-stream methods are computationally preferable, they are often unsuitable for certain physical conditions.For example, non-physical solutions to the twostream method are obtained for the case of a collimated incident beam (Meador & Weaver 1980).However, there exist model adjustments to correct for these limitations in the two-stream method.In particular, applying the delta (δ)-adjustment to the two-stream technique improves the accuracy of radiative flux calculations by taking into account strong forward scattering due to large particles (Joseph et al. 1976;Wiscombe 1977;King 1986;Liou et al. 1988).
Though these corrections often help to improve accuracy, it has also been shown that improved accuracy can be achieved by considering four-stream approximations (Liou 1974;Liou et al. 1988;Cuzzi et al. 1982;Rooney et al. 2023).Four-stream methods are also still able to leverage the δ-approximation to allow for strong forward scattering.However, an increase in the order of approximation comes with the penalty of an increase in computational expense.Four-stream approximations are significantly more efficacious in cases of non-isotropic scattering and can even be used in general circulation models (Liou et al. 1988;Heng & Marley 2018).Their performance has been explored for both homogeneous and inhomogeneous atmospheres in reflected solar radiation (Liou 1973;Liou et al. 1988;Fu 1991;Shibata & Uchiyama 1992).
As well as solar radiation, these approximations can also be applied to study scattering in the presence of thermal emission (Mihalas 1978;Toon et al. 1989;Fu et al. 1997), which is the focus of this work.Infrared scattering is essential to understand emission spectra of exoplanets and brown dwarfs (Taylor et al. 2021), due to the ubiquity of clouds in atmospheres (Marley et al. 2013;Gao & Powell 2021) and the defining role they play in sculpting thermal emission.Specifically, the thermal emission for some classes of extrasolar planets and brown dwarfs (e.g., the L dwarfs) arises, in some wavelengths, from within the scattering, absorbing cloud layers.Therefore, in those cases it is particularly important to treat the radiative transfer within the cloud as carefully as possible.Additionally, the radiative-equilibrium temperature structure of an atmosphere (e.g., Mukherjee et al. 2023) depends upon the difference between upwards and downwards incident and emergent fluxes.Overall, an accurate treatment within the scattering cloud decks is required to have confidence in computed thermal profiles.
Two-stream methods have been implemented for infrared thermal radiation, such as the well-known work of Toon et al. (1989), who derived a general two-stream solution for the upward and downward fluxes within a single homogeneous layer.By considering continuity of flux across a number of stacked, homogeneous layers, the single-layer solution is extended to a multi-layer atmosphere.The final solution is obtained through the two-stream source function technique, with the source function written in terms of the two-stream intensity.As a side note, related to the two-stream technique are popular analytical calculations of temperature-pressure profiles to further understand the thermal atmospheric structure of atmospheres (Hubeny et al. 2003;Hansen et al. 2008;Guillot 2010;Heng & Kopparla 2012;Heng & Showman 2015;Robinson & Catling 2012;Parmentier & Guillot 2014).
The Toon et al. (1989) methodology has in particular been utilized extensively for the study of planetary and substellar atmospheres (e.g.McKay et al. 1989;Marley et al. 1999;Burrows et al. 1997;Fortney et al. 2005;Marley et al. 2021) and this implementation is available in open-source Python code PICASO (Batalha et al. 2022).This approach, however, is currently limited to two-stream approximations.Despite its usefulness in radiative transfer calculations due to its simplicity and ease of implementation, Toon et al. (1989) reported that relative errors in the emissivity calculated by such approaches can be as much as 10% in optically thin cases.
In addition to the potential errors reported by Toon et al. (1989), it has also been shown that increasing the approximations to four-streams improves the accuracy of the infrared models, similar to the reflected solar case (Fu et al. 1997;Liou et al. 1988;Lin et al. 2013).Fu et al. (1997) applied the δ-two and four-stream discrete-ordinates method (Chandrasekhar 1960;Stamnes et al. 1988Stamnes et al. , 2000) ) to solve the infrared radiative transfer equation in a vertically inhomogeneous atmosphere.By comparing the approximations to the high-order δ-129-stream model, the authors found that the δ-two-stream scheme can produce acceptable results under most atmospheric conditions, but suffers from large errors for small optical depth.The δ-four-stream method yields high accuracy in radiative fluxes and heating rates under all atmospheric conditions considered, however, the authors acknowledge a significant increase in computational cost.Zhang et al. (2016) also investigated δ-two and four-stream discrete-ordinates for infrared radiative transfer, demonstrating an analytical approach.By comparing the methods to a δ-64-stream DOM method, the authors similarly conclude that the four-stream method outperforms two-streams, particularly for small optical depths, reporting relative errors as high as 15% for two-stream versus 2% for four-stream.
These studies motivated the development of an analytical spherical harmonics method for solving the radiative transfer equation.The first component for solar radiation was published recently in Rooney et al. (2023).In a similar vein to the Toon et al. (1989) methodology, Rooney et al. (2023) derived and solved a system of equations for the upwards and downwards fluxes at every layer of our atmosphere, with the critical difference of a δ-adjusted four-term spherical harmonics (P 3 ) approximation in place of Toon's two-stream approach.By applying the sourcefunction technique to calculate the azimuthally averaged intensity emerging from the top of a vertically inhomogeneous atmosphere, we can compute the spectrum for clear and cloudy planets or brown dwarfs.Though the spherical harmonics approaches for reflected light and thermal emission are largely identical, the primary differences lie in the source terms, boundary conditions and the related applications.
Therefore, the present work is an extension to the derivation in Rooney et al. (2023), namely, applying the spherical harmonics model to thermal emission spectroscopy.We aim to make this manuscript easily cross-referenced with the numerical method implemented within PICASO source code.Throughout the manuscript, we refer the reader to Rooney et al. (2023) for more intricate detail into the derivation of the spherical harmonics method for 1D radiative transfer, when necessary.Here, we include only the key mathematical expressions that define the thermal emission model.We have also included persistent hyperlinks that can be accessed by clicking the following icon: , that will redirect the reader to the relevant lines of code (stored on Github) corresponding to the relevant mathematical expression.
We outline this work as follows: in Sections 2 and 3, we briefly explain the derivation of the spherical harmonics (SH) method for thermal radiation.As aforementioned, the SH method for reflected light and thermal emission are largely identical, with the exception of the source term and boundary conditions.In this section, we focus on the differences incurred by considering the thermal source term and relevant boundary conditions.We consider both two and four-stream approximations for plane-parallel atmospheres of many layers, where we apply the source-function technique to handle the multi-layer aspect of the model.
In Section 4, we compare the two and four-term spherical harmonics models and the Toon et al. (1989) approach implemented in PICASO to a 16 and 32-stream discrete ordinates method, CDISORT to illustrate the accuracy gains by increasing the number of streams from two to four.We also explore the impact of this order increase on computational time, and discuss the timing-accuracy trade-off that might be considered when choosing a model.

SOLVING THE RADIATIVE TRANSFER EQUATION USING SPHERICAL HARMONICS
We wish to use the spherical harmonics technique to solve the azimuthally-averaged, one-dimensional radiative transfer equation: where the location within the atmosphere is specified by τ ∈ [0, τ N ], (where τ N is the cumulative optical depth), I is the azimuthally averaged intensity and w 0 is the single scattering albedo, B(T ) is the Planck function at temperature T , and P (µ, µ ) is the azimuthally averaged scattering phase function.
We note the similarities between the radiative transfer equation for thermal emission (1) and that for reflected light, outlined in Rooney et al. (2023).The difference lies in the final term on the right-hand side, the source term S(T ), defined as ( We emphasise that all other terms in the azimuthally-averaged, one-dimensional radiative transfer equation ( 1) are identical for reflected light and thermal emission.This allows us to largely follow the spherical harmonics model derivation outlined in Rooney et al. (2023) for reflected light, with a few modifications to allow for the different source term.We will highlight these differences throughout this work.However, we refer the reader to Rooney et al. (2023) for a more in-depth discussion of the general model derivation.
By expanding the phase function and intensity in terms of Legendre polynomials P l , up to given order L: where the coefficients χ l of the phase function expansion can be determined from the orthogonal property of Legendre polynomials (Liou 2002): we can substitute (3) and ( 4) into (1) and use both the orthogonality property and recursion relation of Legendre polynomials to obtain Here, δ 0l is the Dirac-delta function (δ 0l = 1 for l = 0, and 0 otherwise) and Here, a l is identical to that derived for reflected light (Rooney et al. 2023), whereas the expressions for b l differ.This is because the source term ( 2) is relevant only for the b l terms.Thus, any analysis involving only the a l terms and not the b l terms will be identical for reflected light and thermal emission.We assume that the Planck function with a single layer B(T ) can be represented as a Taylor series expansion (as done in Toon et al. 1989), namely where B 0 is the Planck function evaluated at τ = 0 (or the top of the layer) and B 1 is related to the Planck function at temperature T bot at the bottom of the layer τ N :

P 1 Multiple Layers
For clarity and to demonstrate the spherical harmonics methodology, Rooney et al. (2023) began with an atmosphere consisting of a single horizontally homogeneous layer, before extending the analysis to the more practical case of multiple layers.Here, we proceed immediately to the multiple-layer solution.
Let us first study the two-stream spherical harmonics problem, denoted P 1 , where L = 1 represents the highest Legendre polynomial in the expansion.Consider an atmosphere consisting of N horizontally homogeneous layers, where layer n is characterized by single scattering albedo w 0,n , asymmetry parameter g 0,n and optical thickness ∂τ n = τ n − τ n−1 for n = 1, . . ., N .
To solve the radiative transfer equation (1) in the n th layer we rescale the optical depth as Dropping the hats, we continue with the solutions within layer n for τ ∈ [0, ∂τ n ].
We can formulate (6) as a matrix system within layer n: Closely following the methodology outlined in Rooney et al. (2023), we arrive at the layer-wise solution: where is the multi-layer extension of ( 7), and Following Rooney et al. (2023), we can rewrite system (13) in terms of fluxes, where with boundary conditions and where A S is the surface reflectivity.The final boundary condition ( 21) is derived from Mihalas (1978) in Appendix A. These boundary conditions enforce that there is no incident diffuse flux at the top of the atmosphere, and the upward flux at the surface is either that from a (potentially) reflective surface or an estimate of the upwards flux emerging from an atmosphere that continues below the lowermost model grid point (e.g., a giant planet or brown dwarf atmosphere).The spherical harmonics flux problem is formulated in PICASO by representing the system in terms of banded matrices, and solved using the solve banded functionality of SciPy (Virtanen et al. 2020) .

SOURCE FUNCTION TECHNIQUE
Following the methodology of Toon et al. (1989), we apply the source function technique to calculate the emergent intensity from the top of the atmosphere.The radiative transfer equation ( 1) can be solved to yield the azimuthally integrated intensity at angle µ at the top of the n th layer (τ = 0) as for where Toon et al. (1989) showed that infrared intensities can be estimated with sufficient accuracy by using the two-stream approximation to define the source function in the equation of radiative transfer, therefore, we use I t , the solution to the P 1 /P 3 problem outlined in Section 2, in place of the true intensity in the source term (34).Therefore we can rewrite (34) as Let us consider the integral term in (33).Using (36), this can be written as We can calculate the second term on the right-hand side of (37) to be Next, let us write A n,int = ∂τn 0 I l (τ )e − τ µ dτ .This is calculated identically for reflected light, and is outlined in detail in Rooney et al. (2023), where the solution is given by where matrix A n is defined in Rooney et al. ( 2023) and X n are the coefficients we solve for in the P 1 and P 3 flux problems ( 16)-( 21) and ( 25)-( 32) respectively.Vector N n differs from that for reflected light.
For infrared sources, N n is defined as Substituting A n,int (39) and the integrated source term (38) back into (37), we can use (33) to calculate the azimuthally integrated intensity emerging from the top of the n th layer.By beginning at the bottom of the atmosphere (n = N ) and working our way up layer-by-layer, we can derive the azimuthally integrated intensity at the top of the atmosphere.This intensity is used to calculate the infrared flux to predict the observed atmospheric spectra.

ANALYSIS
To quantitatively analyze the performance of the spherical harmonics method for infrared radiative transfer, we compare our results with Toon89 and CDISORT, a version of the discrete ordinate solver, DISORT, written in C rather than FORTRAN (Stamnes et al. 1988(Stamnes et al. , 2000;;Mayer & Kylling 2005;Buras et al. 2011).CDISORT is a versatile, welltested and widely used radiative transfer software, with advanced numerical capabilities.CDISORT has the capacity to model L-stream discrete ordinates approximations, where L is arbitrary and considerably greater than 4 (we will study 16 and 32 stream calculations in this work).
One important note to consider before delving into comparisons of Toon89 and CDISORT is the different scattering phase functions.The Toon89 methodology utilizes the hemispheric mean phase function for infrared scattering (Toon et al. 1989).The hemispheric mean approach is derived by assuming that the phase function takes the value of 1+g 0 in the forward scattering hemisphere, and 1 − g 0 in the backward scattering hemisphere, where g 0 denotes the asymmetry parameter.Toon et al. (1989) chose this technique because for infrared wavelengths, it assumes the correct relationship between flux and intensity and produces the proper emissivity in the limiting case of dominant absorption (w 0 = 0) for a semi-infinite atmosphere.
On the other hand, CDISORT, SH2 and SH4 all utilize the Henyey-Greenstein phase function (Henyey & Greenstein 1941): where the scattering angle Θ is defined as for incoming and outgoing radiation angular directions (µ, φ) and (µ , φ ) respectively.We emphasize this difference in computational methods to foreshadow differences that arise between the methodologies.In what follows, we first compare two benchmark spectra in §4.1.Then, we isolate the dependence of each method's accuracy on scattering parameters (single scattering, asymmetry) in §4.2.Lastly, we baseline the timing of these methodologies in §4.3 to investigate the trade-off between computational time and accuracy.

Comparison of benchmark spectra
We consider two different benchmark atmospheres on which to conduct our analysis: (i) a brown dwarf with effective temperature T eff = 1200 K, gravity g = 200 m/s 2 , solar metallicity, solar C/O and forsterite, iron, corundum clouds, and (ii) planet similar to Jupiter with g =25m/s 2 , semi major axis=5 AU, orbiting a Sun-like star, with H 2 O and NH 3 clouds.The cloudy and non-cloudy infrared spectra, as predicted by the Toon et al. (1989) implementation in PICASO (Batalha et al. 2022), is denoted Toon89.
As shown in Figure 1 the cases chosen both have spectra that are largely affected by the presence of clouds.In both cases, the clouds act to prevent photon contribution from the deepest, hottest, layers.As a result, the pressure range probed by the cloudy models is limited to the upper, cooler layers, creating spectral features that appear muted, relative to the cloud-free counterpart.We choose these cases in order to test the accuracy of these methodologies in the scattering-dominated limit for typical cloudy objects.We note that the different methodologies agree in the case of no scattering, indicating that any spectral differences in the cloudy cases are a consequence of the approximations implemented to deal with scattering.
In Figure 2(b) we plot the infrared spectra for cloudy atmosphere (i), predicted by 16-stream CDISORT, Toon89, twoterm spherical harmonics (SH2) and four-term spherical harmonics (SH4).Note that we indicate whether the models utilise the Henyey-Greenstein (HG) or hemispheric mean (hem-mean) phase function in the figure legend.Figure 2(b) depicts the single scattering, asymmetry and optical depth profiles with pressure, averaged in the wavelength range 1-1.4µm, as indicated by the grey dashed lines on the spectra plot.We choose this wavelength window to average the scattering parameters as it is a region with significant differences in the spectra produced by Toon89 and the other models.The cloud profile is shaped by two cloud layers: one smaller cloud layer below 1 bar overlayed by a larger optical depth cloud layer above 1 bar.In the high optical depth region (τ >0.5), the associated optical properties range between 0.6-0.75 for the asymmetry and 0.8-0.94 for the single scattering.These values are typical for these condensate species and present a less forward scattering example, as compared with the "Jupiter-like" example.
There is an immediately noticeable difference between Toon89, CDISORT and SH2/4.Given the close agreement of SH2 with 16-stream CDISORT compared to Toon89, which is also a two-stream technique, we can isolate that it is   the choice of phase function that leads to this difference.Rooney et al. (2023) conducted an investigation into the accuracy gain when applying four-term spherical harmonics to predict the scattering of reflected light in atmospheres, and compared geometric albedo produced by SH2, SH4 and Toon89 with the doubling method, calculated by Liou (1973).The authors concluded that SH2 and Toon89 performed comparably, and that the choice between spherical harmonics or discrete-ordinates had little impact on the solution accuracy when compared to the doubling method.
The only difference between the SH2 method applied in this work and that applied in Rooney et al. (2023) is the thermal source and boundary conditions, which are identically applied to Toon89 in PICASO.However, the PICASO implementation of Toon89 for reflected light leverages a post-processed Henyey-Greenstein phase function for direct scattering, opposed to the hemispheric mean approach used in infrared (Batalha et al. 2022).
For a clearer understanding of how SH4 compares to SH2, we also plot the percentage difference between 16-stream CDISORT and Toon89, SH2 and SH4 in Figure 3.We notice that the largest deviance of SH2 from 16-stream CDISORT is around 8.5%, whereas SH4 is always within 2.5%.
We conduct the same analysis for the Jupiter-like profile in Figure 4, where the infrared spectra predicted by 16stream CDISORT, Toon89, SH2 and SH4 are plotted on the right, alongside the single scattering, asymmetry and optical depth profiles with pressure, averaged in the wavelength range 8.2-9µm.We notice immediately that the four models are in close agreement throughout the entire wavelength range, with the greatest differences evident between 8-9µm and for wavelengths greater than 13µm.We attribute this to the different cloud condensate optical properties of the Jupiter-like profile.The cloud profile in this case is also shaped by two cloud layers: one larger cloud deck around 1 bar and another smaller optical depth cloud layer around 0.02 bar.In the high optical depth region, the associated optical properties range between 0.8-0.94 for the asymmetry and 0.53-0.8for the single scattering.Since this atmosphere exhibits greater forward scattering but with lower values for the single scattering albedo, this suggests that the accuracy of SH4, SH2, and Toon89 methods are highly dependent on what the strength of scattering and asymmetry of the cloud is.
We again plot the percentage difference between 16-stream CDISORT and Toon89, SH2 and SH4 in Figure 5, and observe a maximum deviance around 11% for Toon89 (ignoring wavelengths less than 6.6µm where spectra values themselves are very small).We also notice that SH2 and SH4 exhibit effectively identical agreement with 16-stream CDISORT, with percentage differences staying below 3.5%.

Dependence of accuracy on scattering parameters
Since it is clear that there is a large accuracy dependence on single scattering and asymmetry, we compare the spectra produced by Toon89 and SH4 with 32-stream CDISORT for a range of single-scattering albedos and asymmetry parameters in a test atmosphere.We define the test atmosphere with the same pressure-temperature and optical depth profile as the T eff = 1270K case studied in Section 4.1, however we force the cloud single-scattering albedo w 0 and asymmetry parameter g 0 to take constant values for clarity.We sweep over a grid of parameter values and calculate  the resulting spectra for each of our three models.We consider w 0 in the range 0.1-1.0, and g 0 in the range 0.0-0.9.Note that we consider a finer grid for high (w 0 >0.9) single-scattering albedo.Taking an average of the spectra values over the wavelength range 1-10µm for each pair of w 0 and g 0 values, we calculate the percentage difference between 32-stream CDISORT and each of Toon89 and SH4.We plot the results as heatmaps in Figure 6, where Figure 6(a) depicts the percentage difference in infrared flux between Toon89 and 32-stream CDISORT, and Figure 6(b) displays that of SH4 and CDISORT.To further elucidate for which parameters SH4 and Toon89 better agree with CDISORT, we subtract the absolute percentage difference of SH4 with CDISORT from that of Toon89, and plot the result in Figure 6(c).In the red-colored regions, SH4 out-performs Toon89 when compared to CDISORT.The white regions represent the cases where SH4 and Toon89 have comparable agreement with CDISORT.The blue-colored regions represent when Toon89 is in closest agreement with CDISORT.We see from Figure 6 that Toon89 exhibits a maximum percentage difference of around 60% with 32-stream CDISORT for high asymmetry (g 0 > 0.8) and moderate single-scattering albedo (0.5 < w 0 < 0.8).The lowest errors occur for extreme values of single scattering, namely w 0 = 0.1 and w 0 > 0.99.We note the excellent agreement for w 0 = 0, which validates the justification of using the hemispheric mean phase function by Toon et al. (1989) to ensure correct emissivity in the w 0 = 0 limit.
Superior agreement is achieved by SH4 when compared with 32-stream CDISORT, with the maximum error of -6% occurring for single-scattering albedo exceeding 0.95.By comparing the two heatmaps in this region, we see that Toon89, using the hemispheric mean approximation, agrees more closely with 32-stream CDISORT than SH4, even though both SH4 and CDISORT use the Henyey-Greenstein phase function.This implies that for high single-scattering (w 0 = 1), the hemispheric-mean approximation marginally outperforms low-order spherical harmonic expansions for the Henyey-Greenstein phase function.

Timing-accuracy trade-off
Despite the improvements of moving to SH4, we still must consider the timing-accuracy trade-off.To elucidate this, we analyze the computational expense of SH2 and SH4 as the number of layers is increased from 40 to 140, alongside the maximum percentage difference of their thermal spectra with that of a 16-stream, 140-layer CDISORT model.We run this analysis on the T eff = 1270K atmosphere studied above, over a wavelength range of 0.7-2µm and plot our results in Figure 7.
As the focus of the present study is to assess the trade-off between computational expense and model accuracy, we compare only SH2 and SH4 to CDISORT to illustrate the increase in cost and improvement in agreement with higher-fidelity models when moving from two to four terms.The modelling approach of SH2 and SH4 is identical bar the number of terms, whereas the Toon89 methodology differs both in model choice (discrete-ordinates versus spherical harmonics) and phase function (hemispheric mean versus Henyey-Greenstein).In an attempt to attribute any differences in computational expense and agreement with CDISORT to only the number of layers chosen for the model, we compare SH2 with SH4.
From Figure 7(a) we see an increase in computational expense, t, for both SH2 and SH4 when the number of layers, N , is increased from 40 to 140 scaling approximately as t = O(N ).Overall, SH4 is approximately twice as expensive as SH2.However, in Figure 7(b) SH4 has a significant increase in model agreement with the CDISORT test case as we increase the number of layers from 40 to 140.For 140 layers, SH4 is within 2% of CDISORT versus 9.7% for SH2, illustrating that although twice as slow as SH2, SH4 is nearly five times more accurate when benchmarked against 16-stream CDISORT.
In cases where model accuracy is important and in the single scattering/asymmetry regions outlined in Section 4.2, SH4 is the obvious choice due to its superior agreement with higher-fidelity models over its lower-term counter-model SH2.However, in the instance when rapid solutions are required, the additional computational expense of SH4 might be undesirable and the efficient SH2 will be the model of choice.

CONCLUSION
Following from the analysis conducted by Rooney et al. (2023), we extended the spherical harmonics approach to solving the radiative transfer equation, implemented in modelling software PICASO Batalha et al. (2022), to thermal emission.In particular, we considered a four-term expansion of spherical harmonics, an increase from the original, two-stream implementation in PICASO, which we denoted Toon89 to reflect its heritage from Toon et al. (1989).The general spherical harmonics methodology for reflected light and thermal emission is largely the same, except for the source function, boundary conditions, and use cases.The main objective of this work was to build on the rigorous derivation of the model for reflected light studied by Rooney et al. (2023), and explain the differences in the model for thermal emission.Without re-deriving every equation in the model, we highlighted the differences incurred by considering a thermal source, and outlined the relevant matrix systems being solved by the model.
To explore the accuracy performance of the four-stream spherical harmonics model, we compared our results to CDISORT (Stamnes et al. 2000).When considered alongside two-term spherical harmonics and the two-stream Toon89 method, this analysis elucidated the increased efficacy of higher-order approximations in radiative transfer calculations for thermal emission, and also demonstrated the impact of the choice of phase function on the resulting spectra.We studied the thermal spectra obtained via the two-stream Toon89 implementation, two and four-term spherical harmonics and CDISORT in Section 4.1 for two different sample atmospheres.This investigation highlighted that the choice of phase function has a large impact on the resultant spectra.The use of hemispheric mean in Toon89 created spectra that were largely different (up to 60%) than those computed from SH2, which utilizes a Henyey-Greenstein phase function.Additionally, we find that accuracy of the order of approximation (two versus four term) is highly dependent on the single scattering albedo and asymmetry of the cloud profile.This motivated a deeper exploration of how the accuracy of the radiative transfer method depends on both values.
Computational time.
Figure 7. Analyzing how the (a) computational time and (b) maximum percentage difference of SH2 and SH4 with CDISORT changes as the number of layers is increased from 40 to 140.We use a 16-stream, 140-layer CDISORT model to benchmark the spherical harmonics against.We see an evident increase in computational expense with number of layers, where SH4 is twice as slow as SH2 for the 140-layer case, however, the maximum percentage difference with the benchmark decreases significantly with layers for SH4, where SH4 is within 2% of CDISORT versus 9.7% for SH2.
Therefore, we created a grid of models with a fixed atmosphere profile and varied asymmetry parameters and singlescattering albedos to study the Toon89 and SH4 models performance when compared to 32-stream CDISORT.We found that Toon89 performs particularly well for the limiting cases of single-scattering albedo, namely w 0 = 0 and w 0 = 1, however suffered from substantial errors of around 60% for high asymmetry and moderate single-scattering.SH4 experiences maximum error of around -6% for high single-scattering.
Finally, we analyzed the timing-accuracy trade-off for the spherical harmonics methods when increasing the number of model layers.By calculating the maximum percentage difference between the thermal spectra produced by SH2 and SH4 with 16-stream, 140-layer CDISORT, we discussed the sacrifice in computational speed for model agreement.This study elucidated that, although increasing model approximation order from two to four terms results in an increase in computational expense, the increase in accuracy when bench-marked against CDISORT is significant.The SH4 model took twice as long as SH2 to calculate the thermal spectra, but produced a result that was nearly five times more accurate when compared to 16-stream CDISORT, with a maximum percentage error of 2%.This analysis demonstrates that a sacrifice of computational expense might be acceptable when a significant increase in accuracy is required from the observational data accuracy, but may not be necessary if numerical efficiency is the priority.
In conclusion, we have demonstrated that increasing the order of approximation from two to four streams can produce significant improvement on model accuracy when compared with high-order CDISORT.The spherical harmonics analysis outlined in this paper is implemented in the PICASO framework, alongside the Toon89 methodology, and is available for download and use Batalha et al. (2022) (a) T eff = 1270K.(b)Jupiter-like.

Figure 1 .
Figure 1.Infrared spectra with and without clouds predicted by Toon89 in PICASO for two different atmospheres: (a) brown dwarf with effective temperature T =1200 K, gravity g =200 m/s 2 , solar metallicity, solar C/O and forsterite, iron, corundum clouds; (b) planet similar to Jupiter with g =25m/s 2 , semi major axis=5 AU, orbiting a Sun-like star, with H2O and NH3 clouds.

Figure 2 .
Figure2.Comparison between the infrared spectra predicted by 16-stream CDISORT, PICASO, 2-term spherical harmonics (SH2) and 4-term spherical harmonics (SH4) for a cloudy brown dwarf with effective temperature T =1200 K, gravity g =200 m/s 2 , solar metallicity, solar C/O and forsterite, iron, corundum clouds.We use HG and hem-mean to indicate that the model phase function is Henyey-Greenstein or hemispheric mean respectively.The scattering properties plotted in subfigure (a) correspond to the average values within the 1-1.4µm wavelength region, as marked by the grey dashed lines on the spectra plot (b).

Figure 3 .
Figure3.Percentage differences between the spectra produced by CDISORT and Toon89, SH2 and SH4 for the T eff = 1270K profile, where HG and hem-mean indicates that the model phase function is Henyey-Greenstein or hemispheric mean respectively.

Figure 4 .
Figure 4. Comparison between the infrared spectra predicted by 16-stream CDISORT, PICASO, 2-term spherical harmonics (SH2) and 4-term spherical harmonics (SH4) for a Jupiter-like profile.The scattering properties plotted in subfigure (a) correspond to the average values within the 8.2-9µm wavelength region, as marked by the grey dashed lines on the spectra plot (b).

Figure 5 .
Figure 5. Differences in spectra for cloudy Jupiter.

Figure 6 .
Figure 6.Heatmaps depicting the percentage difference in average flux produced by (a) Toon89 and (b) SH4 with 32-stream CDISORT. Figure (c) is produced by subtracting the absolute percentage differences of SH4 from that of Toon89 to elucidate for which parameters one method agree with CDISORT better than the other.Dark-red represents cases where SH4 outperforms Toon89, as compared to CDISORT.White (toward zero) represents cases where Toon89 and SH4 perform comparably.
Software: numba(Lam et al. 2015), pandas(McKinney 2010), bokeh (Bokeh Development Team 2014), NumPy(Virtanen et al. 2020),(Walt et al. 2011), IPython(Pérez & Granger 2007), Jupyter,(Kluyver et al. 2016), VIRGA(Batalha et al. 2020;Rooney et al. 2022), PICASO(Batalha et al. 2022), MATLAB (MATLAB 2010), A version of PICASO corresponding to these hyperlinks and the software used in this work is archived on Zenodo as v3.1 with DOI: 10.5281/zenodo.7765171 . The Jupyter notebook, which reproduces our results, can be found on Github as well .C.R.'s research was supported by an appointment to the NASA Postdoctoral Program at the NASA Ames Research Center, administered by Universities Space Research Association under contract with NASA.N.B. & C.R. both acknowledge support from the NASA Astrophysics Division.Additionally, N.B. acknowledges support from NASA'S Interdisciplinary Consortia for Astrobiology Research (NNH19ZDA001N-ICAR) under award number 19-ICAR19 2-0041.We thank Jeff Cuzzi and Sanford Davis for enlightening discussions about some of the finer points of radiative transfer in higher order approximations.Lastly we thank Arve Kylling for helpful discussions regarding CDISORT's radiative transfer methodology.