# Spherical Harmonics for the 1D Radiative Transfer Equation I: Reflected
Light

Caoimhe M. Rooney NASA Ames Research Center, Moffett Field, CA, 94035 USA
Natasha E. Batalha NASA Ames Research Center, Moffett Field, CA, 94035 USA
Mark S. Marley Department of Planetary Sciences, Lunar and Planetary
Laboratory, University of Arizona, Tucson AZ 85721

###### Abstract

A significant challenge in radiative transfer theory for atmospheres of
exoplanets and brown dwarfs is the derivation of computationally efficient
methods that have adequate fidelity to more precise, numerically demanding
solutions. In this work, we extend the capability of the first open-source
radiative transfer model for computing the reflected light of exoplanets at
any phase geometry, PICASO: Planetary Intensity Code for Atmospheric
Spectroscopy Observations. Until now, PICASO has implemented two-stream
approaches to the solving the radiative transfer equation for reflected light,
in particular following the derivations of Toon et al. (1989) (Toon89). In
order to improve the model accuracy, we have considered higher-order
approximations of the phase functions, namely, we have increased the order of
approximation from 2 to 4, using spherical harmonics. The spherical harmonics
approximation decouples spatial and directional dependencies by expanding the
intensity and phase function into a series of spherical harmonics, or Legendre
polynomials, allowing for analytical solutions for low-order approximations to
optimize computational efficiency. We rigorously derive the spherical
harmonics method for reflected light and benchmark the 4-term method (SH4)
against Toon89 and two independent and higher-fidelity methods (CDISORT &
doubling-method). On average, the SH4 method provides an order of magnitude
increase in accuracy, compared to Toon89. Lastly, we implement SH4 within
PICASO and observe only modest increase in computational time, compared to
two-stream methods (20% increase).

Radiative transfer (1335) — Radiative transfer equation (1336)

††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

##  1 Introduction

A serious challenge in atmospheric studies is the derivation of
computationally efficient methods to solve the radiative transfer equation in
a scattering and absorbing medium. Exact solutions typically do not exist;
thus, we rely on approximate methods to estimate solutions (Stephens &
Preisendorfer, 1984; Thomas & Stamnes, 2002; Chandrasekhar, 1960; Liou, 2002).
In recent decades, there has been particular focus on the derivation of simple
and effective approximate methods that have sufficient fidelity to more exact,
numerically intricate solutions. This problem is further complicated by the
vast range of parameter values relevant for atmospheric models that make
spectral inference for exoplanets and Brown Dwarf problems computationally
expensive (e.g. Madhusudhan & Seager, 2009; Line et al., 2012; Barstow et al.,
2020).

The most popular approximate methods for solving the radiative transfer
equation are the (1) discrete-ordinates method (Chandrasekhar, 1960; Stamnes
et al., 1988, 2000), (2) Monte-Carlo method (Modest, 2013; Iwabuchi, 2006) and
(3) spherical harmonics method (Modest, 1989, 2013; Olfe, 1967; van
Wijngaarden & Happer, 2022).

The discrete-ordinates method (DOM) is arguably the most well-studied and
widely used approach when both rapidity and accuracy are required. The problem
to be solved is the transfer of monochromatic radiation in a scattering,
absorbing and emitting plane-parallel atmosphere with a given function for
bidirectional reflectivity at the lower boundary. The general approach is to
discretize the solid angle by a finite number $N$ of directions or “streams”,
along which the radiative intensities are tracked. DISORT (Stamnes et al.,
1988, 2000) is an example of a discrete ordinate algorithm for radiative
transfer in media that is assumed to be non-isothermal, vertically
inhomogeneous, but horizontally homogeneous. DISORT solves the radiative
transfer equation within a single layer without boundary conditions. The
boundary conditions, along with continuity conditions at interfaces between
layers, are used to establish a system of linear equations, which are solved
numerically to obtain the intensity. It is capable of simulating thermal
emission, absorption, and scattering for arbitrary phase functions across the
electromagnetic spectrum from UV to radio wavelengths. However, the
convergence of DOM has been shown to depreciate for optically thick media
(Modest, 2013; Fiveland & Jessee, 1996; Lewis & Miller, 1984). This is due to
the strongly coupled directional equations, which make the DOM method very
computationally expensive when calculating accurate solutions (Modest, 2013;
Ravishankar, 2009). However, there exist a number of acceleration schemes to
improve the convergence rate of DOM (Fiveland & Jessee, 1996; Lewis & Miller,
1984).

Monte-Carlo methods operate by tracking emitted photons throughout the media.
As a renowned accurate method to model radiation, it is often used for bench-
marking, however, it is computationally slow, making it unsuitable for
practical applications (Iwabuchi, 2006; Mayer, 2009). Furthermore, due to the
significant noise that arises as a consequence of the stochastic nature of
Monte-Carlo, the method cannot be easily coupled with other deterministic
solvers.

The spherical harmonics (SH) approximation, denoted $P_{N}$, decouples spatial
and directional dependencies by expanding the intensity and phase function
into a series of spherical harmonics, or Legendre polynomials. This approach
results in fewer equations than the discrete-ordinates method and potentially
obtains higher levels of accuracy at comparable computational expense. This
method has been widely applied to study radiative transfer in planetary and
stellar atmospheres (van Wijngaarden & Happer, 2022; Li & Ramaswamy, 1996;
Zhang & Li, 2013). However, higher order $P_{N}$ are mathematically complex
and increasingly difficult to implement as $N$ increases (Ge et al., 2015; van
Wijngaarden & Happer, 2022). This method is appropriate for scattering media,
however, low-order spherical harmonics can produce inaccuracies for optically
thin atmospheres, and are invalid for non-participating media111Media is
described as non-participating if it does not absorb, emit or scatter
radiation. (Ravishankar, 2009).

Both the discrete-ordinates and spherical harmonics methods allow for
increasing levels of complexity by modifying the order of the approximation.
For DOM, this involves increasing the number of streams, whereas for SH this
involves increasing the order of the Legendre expansion. The simplest and
widely used choice within the discrete-ordinate method is to limit the number
of ordinates, or streams $N$, to two. In other words, this retains only
hemispherical asymmetry information. We refer to the $N=2$ DOM method as the
“two-stream” approximation.

The two-stream problem is formulated as a pair of ordinary differential
equations which are solved within an atmosphere partitioned into homogeneous
layers. Within these layers, the differential equations have constant
coefficients and hence can be solved analytically to produce simple
exponential expressions for fluxes. There exist a number of variations of the
two-stream approach, including the two-stream source function technique (Toon
et al., 1989; Batalha et al., 2019) and the analytical formulation described
in Heng et al. (2014, 2018). The equivalent, lowest order approximation in
spherical harmonics is the $P_{1}$ or Eddington approximation (Meador &
Weaver, 1980; Irvine, 1975; Chandrasekhar, 1960; Mihalas & Mihalas, 2013),
where we expand the intensity and phase function to only two terms.

This low-order angular variation is generally justified for problems where
scattering dominates or is the principle interest; in particular, for slowly
varying phase functions with angle or for large optical depths when intensity
has been smoothed by multiple scattering (Cuzzi et al., 1982). However, there
are significant limitations associated with restricting angular variation to
such low-orders, and these methods are often inappropriate for use with highly
asymmetric phase functions typical of particulate scattering by water and
silicate clouds (Modest, 2013). For example, Wiscombe & Joseph (1977) studied
the accuracy of albedos and absorptivity computed using this Eddington
approximation, for a homogeneous layer with varying values for optical depth,
single-scattering albedo and zenith angle of the incident beam. When compared
to the doubling method, which was assumed by Wiscombe & Joseph (1977) to give
“exact” solutions to the problem, the errors tended to worsen when going from
the Rayleigh to the Henyey-Greenstein to the Mie scattering cases. Errors
generally fell within the 10%–33% range for Mie scattering. This decrease in
accuracy is due to the increase in the forward scattering peak of the phase
function. Such errors in the Mie case are unacceptably large for many
practical applications.

To account for inaccuracies encountered for asymmetric phase functions, the
implementation of the delta ($\delta$)-function adjustment, which does not
increase computational time, was introduced (Joseph et al., 1976). The
$\delta-$adjustment provides a third term closure through the second moment of
the phase function expansion Liou et al. (1988). This adjustment results in an
improvement in the accuracy of radiative flux calculations by accounting for
the strong forward scattering of large particles. Wiscombe & Joseph (1977)
analysed the impact of the $\delta$-adjustment on the accuracy of albedos and
absorptivity mentioned above, and found that the $\delta$-Eddington
approximation reduced the above relative error for Mie scattering of 10%–33%
to less than 10%. Such improvement implies that the forward peak of the Mie
phase function can be better captured using $\delta$-adjusted approximations.

However, King (1986) found that relative errors of 15%–20% could still be
produced for various two-stream approximations. Liou et al. (1988) conducted a
similar analysis for the discrete-ordinate $\delta$-two-stream model. The
$\delta$-two-stream performed comparatively to the spherical harmonics
$\delta$-Eddington approximation studied by King (1986). For conservative
scattering, both approximations display low accuracy for reflection values, on
the order of 10%–30% for incident solar angle $\mu_{0}<0.5$ and $\mu_{0}>0.9$
with $\tau<1$. For transmission, errors exceeded 10% for $\mu_{0}<0.2$.

Further improvement in angular accuracy, whilst maintaining analytical
capacity, can be obtained by increasing the number of streams $N$ or
polynomials to four. Considering $N=4$ in both DOM and SH results in an
improved representation of the phase function and subsequently greater
accuracy of the angle-integrated parameters (Liou, 1973, 1974; Li & Ramaswamy,
1996). The $\delta$-function adjustment to account for the forward diffraction
peak can also be generalised to higher order techniques (Wiscombe, 1977; Cuzzi
et al., 1982; Liou et al., 1988; Shibata & Uchiyama, 1992), including the
4-stream method, and offers significant accuracy improvement. Liou et al.
(1988) conducted the same relative accuracy analysis on the $\delta$-four-
stream approximation (DOM) as discussed above for $\delta$-two-stream and
$\delta$-Eddington. In contrast to the $>$10% errors observed for the
$\delta$-two-stream approximations, the reflection and transmission values
calculated using the $\delta$-four-stream approximation were within 5%
accuracy, except for a number of small regions within the parameter space
where the relative accuracies remained below 10%. Generally, these higher
inaccuracies occurred in optically thin regions with small solar zenith angles
$\mu_{0}<0.3$, where the reflection and transmission values themselves were
very small, thus the absolute errors were extremely small ($<1$%). The authors
concluded that the $\delta$-four stream approximation applied to the
computation of solar fluxes covering the entire solar spectrum can achieve
averaged accuracy within 5%, and thus may be sufficient for studying the flux
distribution for solar radiative transfer in cloudy atmospheres (Liou et al.,
1988; Liou, 1973).

Despite the evidence in the literature suggesting that four-stream methods
offer significant improvements in model accuracy, the conceptual simplicity
and ease of implementation of two-stream methods is attractive for use in
radiative transfer models. A widely used model is that of Toon et al. (1989)
who derived a general two-stream solution for the upward and downward fluxes
within a single homogeneous layer. This single-layer solution is extended to
multiple homogeneous layers through the solving of a matrix system of
equations, with flux continuity conditions enforced at the interfaces between
layers. Following this, the authors introduce the two-stream source function
technique, which involves solving the radiative transfer equation (1) with the
source function written in terms of the two-stream intensity; thus this
approximation of the source function can be inferred from the flux solution
obtained in the first step. With the source function specified, the radiative
transfer equation (1) can be easily solved to give the azimuthally averaged
intensity $I(\tau,\mu)$ at the top ($\tau=0$) and bottom
($\tau=\partial\tau_{n}$) of a layer and thereby the top and bottom of the
atmosphere. The Toon et al. (1989) methodology has been implemented in Python
code PICASO (Batalha et al., 2022). PICASO, however, does not have the ability
to move to higher order streams, and there hasn’t been an extension of Toon et
al. (1989) to four terms.

Here we present a novel approach using spherical harmonics. Similar to the
Toon et al. (1989) methodology, we derive and solve a system of equations for
the layer-wise upwards and downwards fluxes, however we choose a
$\delta$-adjusted four-term spherical harmonics ($P_{3}$) approximation in
place of the two-stream approach. We then extend the source-function technique
for the $P_{3}$ approach to derive the azimuthally averaged intensity emerging
from the top of a vertically inhomogeneous atmosphere. This is then used in an
example calculation to compute the full geometric albedo spectrum for some
cloudy model planets. The spherical harmonics method that we implement is
similar to the analysis conducted by Zhang & Li (2013), however, we extend the
method to multi-layer atmospheres via the source-function technique, rather
than the doubling-adding method studied by Zhang & Li (2013).

We study these with spherical harmonics and compare our results to other
models in the literature. Specifically the applications in atmospheric
modeling often require accurate results in reflection, transmission and
absorption for a wide variety of optical depths and solar zenith angles.
Therefore, this work is outlined as follows: we begin by outlining the
spherical harmonics method for solving the radiative transfer equation within
a single layer in Section 2. We describe both the $P_{1}$ and $P_{3}$
solutions before extending our analysis to multiple layers in Section 3. In
Section 4, we describe our extension to the Toon et al. (1989) source-function
technique. We then conduct a series of model comparisons in Section 5. The
first of these analyses is the comparison of reflection and transmission
values obtained via spherical harmonics to those obtained via the discrete-
ordinate method (DOM) and the doubling method. This comparison follows the
work of Liou (1973) who produced the DOM and doubling data that we use. Next,
we compare various model outputs to that of DISORT, in particular, the
azimuthally averaged intensity varying with both asymmetry and optical depth,
as well as the layer-wise vertical fluxes. Finally, we conclude our model
comparisons with an analysis of the geometric albedo produced by the spherical
harmonics method compared to PICASO’s original two stream Toon et al. (1989)
method.

We intend for this work to be a comprehensive explanation of the spherical
harmonics technique for 1D radiative transfer of reflected light that can
easily be cross-referenced with the numerical method implemented within
PICASO. For this reason, we step through the derivation of the model and
include the key mathematical expressions needed to understand the methodology,
and hence, the algorithm. We have also included 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
mathematical expressions in question.

##  2 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:

$\mu\frac{\partial
I}{\partial\tau}(\tau,\mu)=I(\tau,\mu)-\frac{w_{0}}{2}\int_{-1}^{1}I(\tau,\mu^{\prime})\mathcal{P}(\mu,\mu^{\prime})\mathrm{d}\mu^{\prime}-\frac{w_{0}}{4\pi}F_{\odot}e^{-\frac{\tau}{\mu_{0}}}\mathcal{P}(\mu,-\mu_{0}),$
(1)

where location within the atmosphere is specified by $\tau\in[0,\tau_{N}]$,
(where $\tau_{N}$ is the cumulative optical depth), $I$ is the azimuthally
averaged intensity and $w_{0}$ is the single scattering albedo. The incoming
solar flux is $F_{\odot}$, and the direction of incident solar and outgoing
scattered radiation is defined by the cosine of the zenith angles, denoted
$\mu_{0}$ and $\mu$ respectively. Finally, $\mathcal{P}(\mu,\mu^{\prime})$ is
the azimuthally averaged scattering phase function.

Both the phase function and intensity can be expanded in terms of Legendre
polynomials, up to given order $L$:

$\displaystyle\mathcal{P}(\mu,\mu^{\prime})$
$\displaystyle=\sum_{l=0}^{L}\chi_{l}P_{l}(\mu)P_{l}(\mu^{\prime}),$ (2)
$\displaystyle I(\tau,\mu)$
$\displaystyle=\sum_{l=0}^{L}(2l+1)I_{l}(\tau)P_{l}(\mu),$ (3)

where the coefficients $\chi_{l}$ of the phase function expansion can be
determined from the orthogonal property of Legendre polynomials (Liou, 2002):

$\chi_{l}=\frac{2l+1}{2}\int_{-1}^{1}\mathcal{P}(\cos\Theta)P_{l}(\cos\Theta)\mathrm{d}\cos\Theta.$
(4)

Substituting (2) and (3) into (1) and using both the orthogonality property
and recursion relation of Legendre polynomials, we obtain

$\sum_{l=0}^{L}\left[(l+1)\frac{\mathrm{d}I_{l+1}}{\mathrm{d}\tau}+l\frac{\mathrm{d}I_{l-1}}{\mathrm{d}\tau}\right]P_{l}(\mu)=\sum_{l=0}^{L}[a_{l}I_{l}(\tau)-b_{l}e^{-\frac{\tau}{\mu_{0}}}]P_{l}(\mu),$
(5)

where

$\displaystyle a_{l}=(2l+1)-w_{0}\chi_{l},$ (6) $\displaystyle
b_{l}=\frac{w_{0}\chi_{l}F_{\odot}}{4\pi}P_{l}(-\mu_{0}),$ (7)

for $l=0,\cdots,L$ .

###  2.1 $P_{1}$ single-layer solution

For clarity and to demonstrate the spherical harmonics methodology, we first
consider an atmosphere consisting of a single horizontally homogeneous layer
where the optical properties of which are characterized by its single
scattering albedo $w_{0}$, asymmetry parameter $g_{0}$ and optical thickness
$\tau_{N}$. We will extend this analysis to multiple layers in Section 3. The
two-stream spherical harmonics problem is denoted $P_{1}$, where $L=1$
represents the highest Legendre polynomial in the expansion. We can formulate
(5) as a matrix system given by

$\frac{\mathrm{d}}{\mathrm{d}\tau}\begin{pmatrix}I_{0}\\\
I_{1}\end{pmatrix}=\begin{pmatrix}0&a_{1}\\\
a_{0}&0\end{pmatrix}\begin{pmatrix}I_{0}\\\
I_{1}\end{pmatrix}-\begin{pmatrix}b_{1}\\\
b_{0}\end{pmatrix}e^{-\frac{\tau}{\mu_{0}}},$ (8)

which has solution

$\begin{pmatrix}I_{0}\\\
I_{1}\end{pmatrix}=\begin{pmatrix}e^{-\lambda\tau}&e^{\lambda\tau}\\\
-qe^{-\lambda\tau}&qe^{\lambda\tau}\end{pmatrix}\begin{pmatrix}X_{0}\\\
X_{1}\end{pmatrix}+\begin{pmatrix}\eta_{0}\\\
\eta_{1}\end{pmatrix}e^{-\frac{\tau}{\mu_{0}}},$ (9)

where $\lambda=\sqrt{a_{0}a_{1}}$, $q=\lambda/a_{1}$ and

$\displaystyle\eta_{0}$
$\displaystyle=\frac{1}{1/\mu_{0}^{2}-\lambda^{2}}\left(b_{1}/\mu_{0}-a_{1}b_{0}\right),$
(10) $\displaystyle\eta_{1}$
$\displaystyle=\frac{1}{1/\mu_{0}^{2}-\lambda^{2}}\left(b_{0}/\mu_{0}-a_{0}b_{1}\right).$
(11)

Coefficients $X_{0}$ and $X_{1}$ are to be determined from boundary
conditions. The details of this solution are outlined in Appendix B.

In order to solve for radiative-convective equilibrium temperature structure,
the net upward and downward fluxes are required at every layer in the
atmosphere. Since the ultimate goal of this work is to present the spherical
harmonics method for both incident optical and emitted thermal radiation, we
must calculate such fluxes. As the transformation to fluxes is non-trivial,
especially for non-expert readers, we present the explicit equations.

We can write (9) in terms of fluxes where $F^{\pm}(\tau)=2\pi\int_{0}^{\pm
1}\mu I(\tau,\mu)\,\mathrm{d}\mu$. Recalling the expansion of intensity $I$ in
Legendre polynomials (3), the upward and downward fluxes can be written as

$\displaystyle F^{\pm}(\tau)$ $\displaystyle=2\pi\int_{0}^{\pm
1}\mu[I_{0}(\tau)+3\mu I_{1}(\tau)]\mathrm{d}\mu,$ (12)
$\displaystyle=2\pi\left[\frac{1}{2}I_{0}(\tau)\pm I_{1}(\tau)\right].$ (13)

Therefore, (9) formulated in terms of fluxes is given as

$\begin{pmatrix}F^{-}\\\
F^{+}\end{pmatrix}=\begin{pmatrix}Q^{+}e^{-\lambda\tau}&Q^{-}e^{\lambda\tau}\\\
Q^{-}e^{-\lambda\tau}&Q^{+}e^{\lambda\tau}\end{pmatrix}\begin{pmatrix}X_{0}\\\
X_{1}\end{pmatrix}+\begin{pmatrix}Z^{-}\\\ Z^{+}\end{pmatrix},$ (14)

where $Q^{\pm}=\pi(1\pm 2q)$ and $Z^{\pm}(\tau)=\pi(\eta_{0}\pm 2\eta_{1})$.
For the single-layer problem, we enforce the following boundary conditions
(Toon et al., 1989) on (14)

$\displaystyle F^{-}(0)$ $\displaystyle=0,$ (15) $\displaystyle
F^{+}(\tau_{N})$
$\displaystyle=A_{S}[F^{-}(\tau_{N})+\mu_{0}F_{\odot}e^{-\frac{\tau}{\mu_{0}}}],$
(16)

where $A_{S}$ is the surface reflectivity. 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 the sum of the reflected downward diffuse flux
and the reflection of the unattenuated portion of the direct beam.

###  2.2 $P_{3}$ single-layer

The same principles applied to the $P_{1}$ problem can be extended to study
higher-order techniques. We increase the order of approximation to $L=3$, more
specifically, we consider the $P_{3}$ problem or the four-term spherical
harmonics technique within a single layer. In a similar method to the $P_{1}$
problem in Section 2.1, we can write (5) as

$\frac{\mathrm{d}}{\mathrm{d}\tau}\begin{pmatrix}I_{0}\\\ I_{1}\\\ I_{2}\\\
I_{3}\end{pmatrix}=\begin{pmatrix}0&a_{1}&0&-\frac{2a_{3}}{3}\\\
a_{0}&0&0&0\\\ 0&0&0&\frac{a_{3}}{3}\\\
-\frac{2a_{0}}{3}&0&\frac{a_{2}}{3}&0\\\ \end{pmatrix}\begin{pmatrix}I_{0}\\\
I_{1}\\\ I_{2}\\\ I_{3}\end{pmatrix}-\begin{pmatrix}b_{1}-\frac{2b_{3}}{3}\\\
b_{0}\\\ \frac{b_{3}}{3}\\\
\frac{b_{2}}{3}-\frac{2b_{0}}{3}\end{pmatrix}e^{-\frac{\tau}{\mu_{0}}},$ (17)

which has solution

$\begin{pmatrix}I_{0}\\\ I_{1}\\\ I_{2}\\\
I_{3}\end{pmatrix}=\begin{pmatrix}e^{-\lambda_{1}\tau}&e^{\lambda_{1}\tau}&e^{-\lambda_{2}\tau}&e^{\lambda_{2}\tau}\\\
R_{1}e^{-\lambda_{1}\tau}&-R_{1}e^{\lambda_{1}\tau}&R_{2}e^{-\lambda_{2}\tau}&-R_{2}e^{\lambda_{2}\tau}\\\
Q_{1}e^{-\lambda_{1}\tau}&Q_{1}e^{\lambda_{1}\tau}&Q_{2}e^{-\lambda_{2}\tau}&Q_{2}e^{\lambda_{2}\tau}\\\
S_{1}e^{-\lambda_{1}\tau}&-S_{1}e^{\lambda_{1}\tau}&S_{2}e^{-\lambda_{2}\tau}&-S_{2}e^{\lambda_{2}\tau}\end{pmatrix}\begin{pmatrix}X_{0}\\\
X_{1}\\\ X_{2}\\\ X_{3}\end{pmatrix}+\begin{pmatrix}\eta_{0}\\\ \eta_{1}\\\
\eta_{2}\\\ \eta_{3}\end{pmatrix}e^{-\frac{\tau}{\mu_{0}}}$ (18)

where

$\displaystyle\lambda_{1,2}=\sqrt{\frac{1}{2}(\beta\pm\sqrt{\beta^{2}-4\gamma})},\qquad\beta=a_{0}a_{1}+\frac{1}{9}a_{2}a_{3}+\frac{4}{9}a_{0}a_{3},\qquad\gamma=\frac{1}{9}a_{0}a_{1}a_{2}a_{3},$
(19) $\displaystyle R_{1,2}=-\frac{a_{0}}{\lambda_{1,2}},\qquad
Q_{1,2}=\frac{1}{2}\left(\frac{a_{0}a_{1}}{\lambda_{1,2}^{2}}-1\right),\qquad
S_{1,2}=-\frac{3}{2a_{3}}\left(\frac{a_{0}a_{1}}{\lambda_{1,2}}-\lambda_{1,2}\right),$
(20)

and $\eta_{l}=\Delta_{l}/\Delta$ for $\Delta=9f(1/\mu_{0})$, $f(x)=x^{4}-\beta
x^{2}+\gamma$ and

$\displaystyle\Delta_{0}$
$\displaystyle=(a_{1}b_{0}-b_{1}/\mu_{0})(a_{2}a_{3}-9/\mu_{0}^{2})+2(a_{3}b_{2}-2a_{3}b_{0}-3b_{3}/\mu_{0})/\mu_{0}^{2},$
(21) $\displaystyle\Delta_{1}$
$\displaystyle=(a_{0}b_{1}-b_{0}/\mu_{0})(a_{2}a_{3}-9/\mu_{0}^{2})-2a_{0}(a_{3}b_{2}-3b_{3}/\mu_{0})/\mu_{0},$
(22) $\displaystyle\Delta_{2}$
$\displaystyle=(a_{3}b_{2}-3b_{3}/\mu_{0})(a_{0}a_{1}-1/\mu_{0}^{2})-2a_{3}(a_{0}b_{1}-b_{0}/\mu_{0})/\mu_{0},$
(23) $\displaystyle\Delta_{3}$
$\displaystyle=(a_{2}b_{3}-3b_{2}/\mu_{0})(a_{0}a_{1}-1/\mu_{0}^{2})+2(3a_{0}b_{1}-2a_{0}b_{3}-3b_{0}/\mu_{0})/\mu_{0}^{2}.$
(24)

We wish to consider the problem in terms of fluxes, namely
$F^{\pm}(\tau)=2\pi\int_{0}^{\pm 1}I(\tau,\mu)P_{1}(\mu)\,\mathrm{d}\mu$ and
$f^{\pm}(\tau)=2\pi\int_{0}^{\pm 1}I(\tau,\mu)P_{3}(\mu)\,\mathrm{d}\mu$.
Considering the Legendre expansion of intensity (3), these fluxes are given as

$\displaystyle F^{\pm}(\tau)$ $\displaystyle=2\pi\int_{0}^{\pm
1}\mu\left[I_{0}(\tau)+3\mu
I_{1}(\tau)+\frac{5}{2}(3\mu^{2}-1)I_{2}(\tau)+\frac{7}{2}(5\mu^{3}-3\mu)I_{3}(\tau)\right]\mathrm{d}\mu,$
(25) $\displaystyle=2\pi\left(I_{0}(\tau)\pm
I_{1}(\tau)+\frac{5}{8}I_{2}(\tau)\right),$ (26) $\displaystyle f^{\pm}(\tau)$
$\displaystyle=2\pi\int_{0}^{\pm
1}\frac{1}{2}\left(5\mu^{3}-3\mu\right)\left[I_{0}(\tau)+3\mu
I_{1}(\tau)+\frac{5}{2}(3\mu^{2}-1)I_{2}(\tau)+\frac{7}{2}(5\mu^{3}-3\mu)I_{3}(\tau)\right]\mathrm{d}\mu,$
(27)
$\displaystyle=2\pi\left(-\frac{1}{8}I_{0}(\tau)+\frac{5}{8}I_{2}(\tau)\pm
I_{3}(\tau)\right).$ (28)

System (18) can therefore be rewritten in terms of fluxes as

$\begin{pmatrix}F^{-}\\\ f^{-}\\\ F^{+}\\\
f^{-}\end{pmatrix}=\begin{pmatrix}p_{1}^{-}e^{-\lambda_{1}\tau}&p_{1}^{+}e^{\lambda_{1}\tau}&p_{2}^{-}e^{-\lambda_{2}\tau}&p_{2}^{+}e^{\lambda_{2}\tau}\\\
q_{1}^{-}e^{-\lambda_{1}\tau}&q_{1}^{+}e^{\lambda_{1}\tau}&q_{2}^{-}e^{-\lambda_{2}\tau}&q_{2}^{+}e^{\lambda_{2}\tau}\\\
p_{1}^{+}e^{-\lambda_{1}\tau}&p_{1}^{-}e^{\lambda_{1}\tau}&p_{2}^{+}e^{-\lambda_{2}\tau}&p_{2}^{-}e^{\lambda_{2}\tau}\\\
q_{1}^{+}e^{-\lambda_{1}\tau}&q_{1}^{-}e^{\lambda_{1}\tau}&q_{2}^{+}e^{-\lambda_{2}\tau}&q_{2}^{-}e^{\lambda_{2}\tau}\\\
\end{pmatrix}\begin{pmatrix}X_{0}\\\ X_{1}\\\ X_{2}\\\
X_{3}\end{pmatrix}+\begin{pmatrix}Z_{1}^{-}\\\ Z_{2}^{-}\\\ Z_{1}^{+}\\\
Z_{2}^{+}\end{pmatrix}$ (29)

where $p_{1,2}^{\pm}=\pi(1\pm 2R_{1,2}+\frac{5}{4}Q_{1,2})$,
$q_{1,2}^{\pm}=\pi(-\frac{1}{4}+\frac{5}{4}Q_{1,2}\pm 2S_{1,2})$,
$Z_{1}^{\pm}(\tau)=\pi(\eta_{0}\pm 2\eta_{1}+\frac{5}{4}\eta_{2})$,
$Z_{2}^{\pm}(\tau)=\pi(-\frac{1}{4}\eta_{0}+\frac{5}{4}\eta_{2}\pm
2\eta_{3})$.

For the single-layer problem, we enforce the following boundary conditions on
(29),

$\displaystyle F^{-}(0)$ $\displaystyle=0,$ (30) $\displaystyle f^{-}(0)$
$\displaystyle=0,$ (31) $\displaystyle F^{+}(\tau_{N})$
$\displaystyle=A_{S}[F^{-}(\tau_{N})+\mu_{0}F_{\odot}e^{-\frac{\tau}{\mu_{0}}}],$
(32) $\displaystyle f^{+}(\tau_{N})$
$\displaystyle=A_{S}[f^{-}(\tau_{N})-\frac{1}{4}\mu_{0}F_{\odot}e^{-\frac{\tau}{\mu_{0}}}],$
(33)

where $A_{S}$ is the surface reflectivity. The derivation of the bottom
boundary condition (33) for $f^{+}$ is outlined in Appendix G.

##  3 Extension to multiple layers

The procedures outlined in Sections 2.1 and 2.2 can be extended to multiple
layers by solving the flux systems (14) for $P_{1}$ and (29) for $P_{3}$
within each layer and enforcing continuity boundary conditions between layers
along with the boundary conditions given above for the top and bottom of the
atmosphere.

We divide the atmosphere into $N$ homogeneous layers of optical depth
$\partial\tau_{n}=\tau_{n}-\tau_{n-1}$ for $n=1,\dots,N$. To solve the
radiative transfer equation (1) in the $n^{\text{th}}$ layer we rescale the
optical depth as

$\hat{\tau}=\tau-\tau_{n-1},\qquad\hat{\tau}\in[0,\partial\tau_{n}].$ (34)

Dropping the hats, we continue with the solutions within layer $n$ for
$\tau\in[0,\partial\tau_{n}]$.

###  3.1 $P_{1}$ multiple layers

For $N$ layers, we must solve system (8) within every homogeneous layer. The
solution at the top of layer $n$ is given by

$\begin{pmatrix}I_{0,n}\\\
I_{1,n}\end{pmatrix}=\begin{pmatrix}e^{-\lambda_{n}\tau}&e^{\lambda_{n}\tau}\\\
-q_{n}e^{-\lambda_{n}\tau}&q_{n}e^{\lambda_{n}\tau}\end{pmatrix}\begin{pmatrix}X_{0,n}\\\
X_{1,n}\end{pmatrix}-\begin{pmatrix}\eta_{0,n}\\\
\eta_{1,n}\end{pmatrix}e^{-\frac{1}{\mu_{0}}(\tau+\tau_{n-1})},$ (35)

for $\tau\in[0,\partial\tau_{n}]$, and where we extend the single-layer
definitions for $\eta$ given by (11) to obtain

$\displaystyle\eta_{0,n}$
$\displaystyle=\frac{1}{1/\mu_{0}^{2}-\lambda_{n}^{2}}\left(b_{1,n}/\mu_{0}-a_{1,n}b_{0,n}\right),$
(36) $\displaystyle\eta_{1,n}$
$\displaystyle=\frac{1}{1/\mu_{0}^{2}-\lambda_{n}^{2}}\left(b_{0,n}/\mu_{0}-a_{0,n}b_{1,n}\right),$
(37)

and $a_{l,n}=(2l+1)-w_{0,n}\chi_{l,n}$,
$b_{l,n}=w_{0,n}\chi_{l,n}F_{\odot}P_{l}(-\mu_{0})/4\pi$,
$\lambda_{n}=\sqrt{a_{0,n}a_{1,n}}$, and $q_{n}=\lambda/a_{1,n}$ , similarly
extended from single-layer definitions (7). In terms of fluxes,

$\begin{pmatrix}F_{n}^{-}\\\
F_{n}^{+}\end{pmatrix}=\begin{pmatrix}Q_{n}^{+}e^{-\lambda_{n}\tau}&Q_{n}^{-}e^{\lambda_{n}\tau}\\\
Q_{n}^{-}e^{-\lambda_{n}\tau}&Q_{n}^{+}e^{\lambda_{n}\tau}\end{pmatrix}\begin{pmatrix}X_{0,n}\\\
X_{1,n}\end{pmatrix}+\begin{pmatrix}Z_{n}^{-}\\\
Z_{n}^{+}\end{pmatrix}e^{-\frac{1}{\mu_{0}}(\tau+\tau_{n-1})},$ (38)

where $Q_{n}^{\pm}=\pi(1\pm 2q_{n})$, $Z_{n}^{\pm}(\tau)=\pi(\eta_{0,n}\pm
2\eta_{1,n})$ , and with boundary conditions

$\displaystyle F_{1}^{-}(0)$ $\displaystyle=0,$ (39) $\displaystyle
F_{n}^{-}(\partial\tau_{n})$ $\displaystyle=F_{n+1}^{-}(0),$ (40)
$\displaystyle F_{n}^{+}(\partial\tau_{n})$ $\displaystyle=F_{n+1}^{+}(0),$
(41) $\displaystyle F_{N}^{+}(\tau_{N})$
$\displaystyle=A_{S}[F_{N}^{-}(\tau_{N})+\mu_{0}F_{\odot}e^{-\frac{\tau_{N}}{\mu_{0}}}].$
(42)

The 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).

###  3.2 $P_{3}$ multiple layers

Consider $L=3$. In a similar method to the $P_{1}$ problem, the solution to
the radiative transfer equation within layer $n$ is given by

$\begin{pmatrix}I_{0,n}\\\ I_{1,n}\\\ I_{2,n}\\\
I_{3,n}\end{pmatrix}=\begin{pmatrix}e^{-\lambda_{1,n}\tau}&e^{\lambda_{1,n}\tau}&e^{-\lambda_{2,n}\tau}&e^{\lambda_{2,n}\tau}\\\
R_{1,n}e^{-\lambda_{1,n}\tau}&-R_{1,n}e^{\lambda_{1,n}\tau}&R_{2,n}e^{-\lambda_{2,n}\tau}&-R_{2,n}e^{\lambda_{2,n}\tau}\\\
Q_{1,n}e^{-\lambda_{1,n}\tau}&Q_{1,n}e^{\lambda_{1,n}\tau}&Q_{2,n}e^{-\lambda_{2,n}\tau}&Q_{2,n}e^{\lambda_{2,n}\tau}\\\
S_{1,n}e^{-\lambda_{1,n}\tau}&-S_{1,n}e^{\lambda_{1,n}\tau}&S_{2,n}e^{-\lambda_{2,n}\tau}&-S_{2,n}e^{\lambda_{2,n}\tau}\end{pmatrix}\begin{pmatrix}X_{0,n}\\\
X_{1,n}\\\ X_{2,n}\\\ X_{3,n}\end{pmatrix}-\begin{pmatrix}\eta_{0,n}\\\
\eta_{1,n}\\\ \eta_{2,n}\\\
\eta_{3,n}\end{pmatrix}e^{-\frac{1}{\mu_{0}}(\tau+\tau_{n-1})},$ (43)

for $\tau\in[0,\partial\tau_{n}]$, where

$\displaystyle\lambda_{1,2,n}=\sqrt{\frac{1}{2}(\beta_{n}\pm\sqrt{\beta_{n}^{2}-4\gamma_{n}})},\qquad\beta_{n}=a_{0,n}a_{1,n}+\frac{1}{9}a_{2,n}a_{3,n}+\frac{4}{9}a_{0,n}a_{3,n},\qquad\gamma_{n}=\frac{1}{9}a_{0,n}a_{1,n}a_{2,n}a_{3,n},$
(44)

and

$\displaystyle R_{1,2,n}=-\frac{a_{0,n}}{\lambda_{1,2,n}},\qquad
Q_{1,2,n}=\frac{1}{2}\left(\frac{a_{0,n}a_{1,n}}{\lambda_{1,2,n}^{2}}-1\right),\qquad
S_{1,2,n}=-\frac{3}{2a_{3,n}}\left(\frac{a_{0,n}a_{1,n}}{\lambda_{1,2,n}}-\lambda_{1,2,n}\right).$
(45)

We also define $\eta_{i,n}=\Delta_{i,n}/\Delta$ , for $\Delta=9f(1/\mu_{0})$,
$f(x)=x^{4}-\beta_{n}x^{2}+\gamma_{n}$ and

$\displaystyle\Delta_{0,n}$
$\displaystyle=(a_{1,n}b_{0,n}-b_{1,n}/\mu_{0})(a_{2,n}a_{3,n}-9/\mu_{0}^{2})+2(a_{3,n}b_{2,n}-2a_{3,n}b_{0,n}-3b_{3,n}/\mu_{0})/\mu_{0}^{2},$
(46) $\displaystyle\Delta_{1,n}$
$\displaystyle=(a_{0,n}b_{1,n}-b_{0,n}/\mu_{0})(a_{2,n}a_{3,n}-9/\mu_{0}^{2})-2a_{0,n}(a_{3,n}b_{2,n}-3b_{3,n}/\mu_{0})/\mu_{0},$
(47) $\displaystyle\Delta_{2,n}$
$\displaystyle=(a_{3,n}b_{2,n}-3b_{3,n}/\mu_{0})(a_{0,n}a_{1,n}-1/\mu_{0}^{2})-2a_{3,n}(a_{0,n}b_{1,n}-b_{0,n}/\mu_{0})/\mu_{0},$
(48) $\displaystyle\Delta_{3,n}$
$\displaystyle=(a_{2,n}b_{3,n}-3b_{2,n}/\mu_{0})(a_{0,n}a_{1,n}-1/\mu_{0}^{2})+2(3a_{0,n}b_{1,n}-2a_{0,n}b_{3,n}-3b_{0,n}/\mu_{0})/\mu_{0}^{2}.$
(49)

This problem can be written in terms of fluxes as

$\begin{pmatrix}F_{n}^{-}\\\ f_{n}^{-}\\\ F_{n}^{+}\\\
f_{n}^{-}\end{pmatrix}=\begin{pmatrix}p_{1,n}^{-}e^{-\lambda_{1,n}\tau}&p_{1,n}^{+}e^{\lambda_{1,n}\tau}&p_{2,n}^{-}e^{-\lambda_{2,n}\tau}&p_{2,n}^{+}e^{\lambda_{2,n}\tau}\\\
q_{1,n}^{-}e^{-\lambda_{1,n}\tau}&q_{1,n}^{+}e^{\lambda_{1,n}\tau}&q_{2,n}^{-}e^{-\lambda_{2,n}\tau}&q_{2,n}^{+}e^{\lambda_{2,n}\tau}\\\
p_{1,n}^{+}e^{-\lambda_{1,n}\tau}&p_{1,n}^{-}e^{\lambda_{1,n}\tau}&p_{2,n}^{+}e^{-\lambda_{2,n}\tau}&p_{2,n}^{-}e^{\lambda_{2,n}\tau}\\\
q_{1,n}^{+}e^{-\lambda_{1,n}\tau}&q_{1,n}^{-}e^{\lambda_{1,n}\tau}&q_{2,n}^{+}e^{-\lambda_{2,n}\tau}&q_{2,n}^{-}e^{\lambda_{2,n}\tau}\\\
\end{pmatrix}\begin{pmatrix}X_{0,n}\\\ X_{1,n}\\\ X_{2,n}\\\
X_{3,n}\end{pmatrix}+\begin{pmatrix}Z_{1,n}^{-}\\\ Z_{2,n}^{-}\\\
Z_{1,n}^{+}\\\
Z_{2,n}^{+}\end{pmatrix}e^{-\frac{1}{\mu_{0}}(\tau+\tau_{n-1})},$ (50)

where $p_{1,2,n}^{\pm}=\pi(1\pm 2R_{1,2,n}+\frac{5}{4}Q_{1,2,n})$,
$q_{1,2,n}^{\pm}=\pi(-\frac{1}{4}+\frac{5}{4}Q_{1,2,n}\pm 2S_{1,2,n})$ ,
$Z_{1,n}^{\pm}=\pi(\eta_{0,n}\pm 2\eta_{1,n}+\frac{5}{4}\eta_{2,n})$,
$Z_{2^{,}n}\pm=\pi(-\frac{1}{4}\eta_{0,n}+\frac{5}{4}\eta_{2,n}\pm
2\eta_{3,n})$ , and with boundary conditions

$\displaystyle F_{1}^{-}(0)$ $\displaystyle=0,$ $\displaystyle\quad\qquad
f_{1}^{-}(0)=0,$ (51) $\displaystyle F_{n}^{-}(\partial\tau_{n})$
$\displaystyle=F_{n+1}^{-}(0),$ $\displaystyle\qquad
f_{n}^{-}(\partial\tau_{n})=f_{n+1}^{-}(0),$ (52) $\displaystyle
F_{n}^{+}(\partial\tau_{n})$ $\displaystyle=F_{n+1}^{+}(0),$
$\displaystyle\qquad f_{n}^{+}(\partial\tau_{n})=f_{n+1}^{+}(0),$ (53)
$\displaystyle F_{N}^{+}(\tau_{N})$
$\displaystyle=A_{S}[F^{-}_{N}(\tau_{N})+\mu_{0}F_{\odot}e^{-\frac{\tau_{N}}{\mu_{0}}}],$
$\displaystyle\,\,\qquad
f_{N}^{+}(\tau_{N})=A_{S}[f_{N}^{-}(\tau_{N})-\frac{1}{4}\mu_{0}F_{\odot}e^{-\frac{\tau}{\mu_{0}}}]$
(54)

for $n={1,2,\cdots,N-1}$ . The derivation of the bottom boundary condition for
$f_{N}^{+}(\tau_{N})$ is outlined in Appendix G. As for the $P_{1}$ case, the
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).

##  4 Source function technique

We follow the methodology outlined in Toon et al. (1989) and 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 $\mu$ at the top of the
$n^{\text{th}}$ layer ($\tau=0$) as

$I_{n}(0,\mu)=I_{n}(\partial\tau_{n},\mu)e^{-\frac{\partial\tau_{n}}{\mu}}+\frac{1}{\mu}\int_{0}^{\partial\tau_{n}}S_{vt}e^{-\frac{\tau}{\mu}}\mathrm{d}\tau,$
(55)

for

$S_{vt}=\frac{w_{0,n}}{2}\int_{-1}^{1}I_{t}(\tau,\mu^{\prime})\mathcal{P}(\mu,\mu^{\prime})\mathrm{d}\mu^{\prime}+S_{n}(\tau,\mu),$
(56)

where

$S_{n}(\tau,\mu)=\frac{w_{0}}{4\pi}F_{\odot}e^{-\frac{\tau+\tau_{n-1}}{\mu_{0}}}\mathcal{P}(\mu,-\mu_{0}).$
(57)

Note that we have used $I_{t}$ to denote the approximated intensity (3) in
place of the true intensity in the source term (56). Therefore we can rewrite
(56) as

$S_{vt}=w_{0,n}\sum_{l=0}^{L}\chi_{l}I_{l}(\tau)P_{l}(\mu)+S_{n}(\tau,\mu).$
(58)

Let us consider the integral term in (55). Using (58), this can be written as

$\int_{0}^{\partial\tau_{n}}S_{vt}e^{-\frac{\tau}{\mu}}\mathrm{d}\tau=w_{0}\sum_{l=0}^{L}\chi_{l}P_{l}(\mu)\int_{0}^{\partial\tau_{n}}I_{l}(\tau)e^{-\frac{\tau}{\mu}}\mathrm{d}\tau+\int_{0}^{\partial\tau_{n}}S_{n}(\tau,\mu)e^{-\frac{\tau}{\mu}}\mathrm{d}\tau.$
(59)

We can calculate the second term on the right-hand side of (59) to be

$\int_{0}^{\partial\tau_{n}}S_{n}(\tau,\mu)e^{-\frac{\tau}{\mu}}\mathrm{d}\tau=\frac{w_{0}}{4\pi}F_{\odot}\mathcal{P}(\mu,-\mu_{0})\frac{\mu\mu_{0}}{\mu+\mu_{0}}\left(1-e^{-\partial\tau_{n}\frac{\mu+\mu_{0}}{\mu\mu_{0}}}\right)e^{-\frac{\tau_{n}}{\mu_{0}}}.$
(60)

Next, let us write
$A_{n,\text{int}}=\int_{0}^{\partial\tau_{n}}I_{l}(\tau)e^{-\frac{\tau}{\mu}}\mathrm{d}\tau$.
Therefore we can calculate $A_{n,\text{int}}$ for the $P_{1}$ problem :

$A_{n,\text{int}}=\begin{pmatrix}A_{0,n}&A_{1,n}\\\ -QA_{0,n}&QA_{1,n}\\\
\end{pmatrix}\begin{pmatrix}X_{0,n}\\\
X_{1,n}\end{pmatrix}+\begin{pmatrix}N_{0,n}\\\
N_{1,n}\end{pmatrix}e^{-\frac{\tau_{n}}{\mu_{0}}}$ (61)

where

$\displaystyle A_{0,n}$
$\displaystyle=\frac{1}{1/\mu+\lambda_{n}}\left(1-e^{-\partial\tau_{n}\left(\frac{1}{\mu}+\lambda_{n}\right)}\right),$
(62) $\displaystyle A_{1,n}$
$\displaystyle=\frac{1}{1/\mu-\lambda_{n}}\left(1-e^{-\partial\tau_{n}\left(\frac{1}{\mu}-\lambda_{n}\right)}\right),$
(63)

and

$N_{i,n}=\eta_{i,n}\frac{\mu\mu_{0}}{\mu+\mu_{0}}\left(1-e^{-\partial\tau_{n}\frac{\mu+\mu_{0}}{\mu\mu_{0}}}\right)$
(64)

where $n_{i,n}=\Delta_{i,n}/\Delta$ are defined by equations (46)–(49).

Similarly for the $P_{3}$ problem :

$A_{n,\text{int}}=\begin{pmatrix}A_{00,n}&A_{01,n}&A_{02,n}&A_{03,n}\\\
A_{10,n}&A_{11,n}&A_{12,n}&A_{13,n}\\\ A_{20,n}&A_{21,n}&A_{22,n}&A_{23,n}\\\
A_{30,n}&A_{31,n}&A_{32,n}&A_{33,n}\\\ \end{pmatrix}\begin{pmatrix}X_{0,n}\\\
X_{1,n}\\\ X_{2,n}\\\ X_{3,n}\end{pmatrix}+\begin{pmatrix}N_{0,n}\\\
N_{1,n}\\\ N_{2,n}\\\ N_{3,n}\end{pmatrix}e^{-\frac{\tau_{n}}{\mu_{0}}},$ (65)

where

$\displaystyle A_{00,n}$ $\displaystyle=\alpha_{1,n},\quad\qquad
A_{01,n}=\beta_{1,n},\,\,\,\quad\qquad A_{02,n}=\alpha_{2,n},\quad\qquad
A_{03,n}=\beta_{2,n},$ (66) $\displaystyle A_{10,n}$
$\displaystyle=R_{1,n}\alpha_{1,n},\quad A_{11,n}=-R_{1,n}\beta_{1,n},\quad
A_{12,n}=R_{2,n}\alpha_{2,n},\quad A_{13,n}=-R_{2,n}\beta_{2,n},$ (67)
$\displaystyle A_{20,n}$ $\displaystyle=Q_{1,n}\alpha_{1,n},\quad
A_{21,n}=Q_{1,n}\beta_{1,n},\,\,\,\,\quad
A_{22,n}=Q_{2}\alpha_{2,n},\,\,\,\,\quad A_{23,n}=Q_{2,n}\beta_{2,n},$ (68)
$\displaystyle A_{30,n}$ $\displaystyle=S_{1,n}\alpha_{1,n},\,\quad
A_{31,n}=-S_{1,n}\beta_{1,n},\,\quad
A_{32,n}=S_{2}\alpha_{2,n},\,\,\,\,\,\quad A_{33,n}=-S_{2,n}\beta_{2,n},$ (69)

and

$\displaystyle\alpha_{i,n}=\frac{1}{1/\mu+\lambda_{i,n}}\left(1-e^{-\partial\tau_{n}\left(\frac{1}{\mu}+\lambda_{i,n}\right)}\right),\qquad\beta_{i,n}=\frac{1}{1/\mu-\lambda_{i,n}}\left(1-e^{-\partial\tau_{n}\left(\frac{1}{\mu}-\lambda_{i,n}\right)}\right).$
(71)

Coefficients $R_{i,n},Q_{i,n},S_{i,n}$ are defined using equations (45), and
$N_{i,n}$ using (64).

##  5 Analysis

###  5.1 Liou Comparison

|  |  | Reflection ($w_{0}=1$) | Transmisssion ($w_{0}=1$)
---|---|---|---|---
|  |  |  |  | $\mu_{0}$ | $\mu_{0}$
$\tau$ |  | Method |  | 0.1 |  | 0.5 |  | 0.9 |  |  | 0.1 |  | 0.5 |  | 0.9 |
0.25 |  |  | Toon89 |  |  | 0.41133 |  | 0.07635 |  | -0.01294 |  |  | 0.58867 |  | 0.92365 |  | 1.01294 |
|  | SH2 |  |  | 0.41742 |  | 0.09176 |  | 0.00423 |  |  | 0.58258 |  | 0.90824 |  | 0.99577 |
|  | SH4 |  |  | 0.40745 |  | 0.06746 |  | 0.02331 |  |  | 0.59255 |  | 0.93254 |  | 0.97669 |
|  | Doubling |  |  | 0.41610 |  | 0.07179 |  | 0.02250 |  |  | 0.58390 |  | 0.92821 |  | 0.97751 |
1 |  |  | Toon89 |  |  | 0.51778 |  | 0.07635 |  | 0.02389 |  |  | 0.48222 |  | 0.77441 |  | 0.97610 |
|  | SH2 |  |  | 0.51577 |  | 0.09176 |  | 0.05904 |  |  | 0.48423 |  | 0.75109 |  | 0.94096 |
|  | SH4 |  |  | 0.55781 |  | 0.06746 |  | 0.09758 |  |  | 0.44218 |  | 0.76739 |  | 0.90242 |
|  | Doubling |  |  | 0.58148 |  | 0.07179 |  | 0.09672 |  |  | 0.41852 |  | 0.75952 |  | 0.90328 |
4 |  |  | Toon89 |  |  | 0.68564 |  | 0.07635 |  | 0.31611 |  |  | 0.31436 |  | 0.50001 |  | 0.68388 |
|  | SH2 |  |  | 0.67143 |  | 0.09176 |  | 0.32974 |  |  | 0.32857 |  | 0.50002 |  | 0.67025 |
|  | SH4 |  |  | 0.72698 |  | 0.06746 |  | 0.34904 |  |  | 0.27301 |  | 0.48146 |  | 0.65095 |
|  | Doubling |  |  | 0.73254 |  | 0.07179 |  | 0.34823 |  |  | 0.26746 |  | 0.48069 |  | 0.65178 |
16 |  |  | Toon89 |  |  | 0.86859 |  | 0.07635 |  | 0.71338 |  |  | 0.13140 |  | 0.20899 |  | 0.28659 |
|  | SH2 |  |  | 0.85624 |  | 0.09176 |  | 0.70623 |  |  | 0.14374 |  | 0.21874 |  | 0.29373 |
|  | SH4 |  |  | 0.87891 |  | 0.06746 |  | 0.70853 |  |  | 0.12107 |  | 0.21351 |  | 0.29143 |
|  | Doubling |  |  | 0.88103 |  | 0.07179 |  | 0.70722 |  |  | 0.11897 |  | 0.21342 |  | 0.29279 |
|  |  | Reflection ($w_{0}=0.8$) | Transmisssion ($w_{0}=0.8$)
|  |  |  |  | $\mu_{0}$ | $\mu_{0}$
$\tau$ |  | Method |  | 0.1 |  | 0.5 |  | 0.9 |  |  | 0.1 |  | 0.5 |  | 0.9 |
0.25 |  |  | Toon89 |  |  | 0.31802 |  | 0.05739 |  | -0.01125 |  |  | 0.46566 |  | 0.84979 |  | 0.95403 |
|  | SH2 |  |  | 0.31987 |  | 0.06802 |  | 0.00111 |  |  | 0.45934 |  | 0.83725 |  | 0.94049 |
|  | SH4 |  |  | 0.30146 |  | 0.04662 |  | 0.01676 |  |  | 0.45385 |  | 0.85230 |  | 0.92677 |
|  | Doubling |  |  | 0.28961 |  | 0.04855 |  | 0.01547 |  |  | 0.43017 |  | 0.84756 |  | 0.92669 |
1 |  |  | Toon89 |  |  | 0.37519 |  | 0.05739 |  | -0.00064 |  |  | 0.29023 |  | 0.55267 |  | 0.76333 |
|  | SH2 |  |  | 0.36348 |  | 0.06802 |  | 0.01710 |  |  | 0.28732 |  | 0.53231 |  | 0.73607 |
|  | SH4 |  |  | 0.36158 |  | 0.04662 |  | 0.05275 |  |  | 0.22624 |  | 0.52033 |  | 0.72010 |
|  | Doubling |  |  | 0.35487 |  | 0.04855 |  | 0.04929 |  |  | 0.20556 |  | 0.51606 |  | 0.71772 |
4 |  |  | Toon89 |  |  | 0.40411 |  | 0.05739 |  | 0.05152 |  |  | 0.06605 |  | 0.12281 |  | 0.20828 |
|  | SH2 |  |  | 0.38034 |  | 0.06802 |  | 0.05024 |  |  | 0.06584 |  | 0.11729 |  | 0.19690 |
|  | SH4 |  |  | 0.38125 |  | 0.04662 |  | 0.09327 |  |  | 0.04441 |  | 0.10603 |  | 0.22080 |
|  | Doubling |  |  | 0.37148 |  | 0.04855 |  | 0.08925 |  |  | 0.04539 |  | 0.10718 |  | 0.21953 |
16 |  |  | Toon89 |  |  | 0.40571 |  | 0.05739 |  | 0.05636 |  |  | 0.00018 |  | 0.00034 |  | 0.00060 |
|  | SH2 |  |  | 0.38128 |  | 0.06802 |  | 0.05297 |  |  | 0.00018 |  | 0.00033 |  | 0.00057 |
|  | SH4 |  |  | 0.38204 |  | 0.04662 |  | 0.09700 |  |  | 0.00027 |  | 0.00063 |  | 0.00138 |
|  | Doubling |  |  | 0.37229 |  | 0.04855 |  | 0.09297 |  |  | 0.00027 |  | 0.00062 |  | 0.00139 |
Table 1: Comparison of reflection and transmission as computed by Toon89,
2-term spherical harmonics (SH2), 4-term spherical harmonics (SH4) and by the
doubling method for conservative and non-conservative scattering. The doubling
method results are taken from Liou (1973). Figure 1: To illustrate model
performance, we calculated the percentage difference between each model
(Toon89, 2 and 4-term spherical harmonics [SH2, SH4], 2 and 4-stream discrete
ordinates [DOM2, DOM4]) and the doubling method for reflection and
transmission. The raw data for Toon89 and spherical harmonics is given in
Table 1 and the DOM and doubling method results are taken from Liou (1973),
where the results for the doubling method are considered to be the “true”
solution.

To quantitatively justify the effective performance of the spherical harmonics
method for radiative transfer, we compare it with established methods in the
literature. The physical quantities of interest for radiative balance studies
in planetary atmospheres are the vertical distributions of fluxes. To study
the capabilities of the spherical harmonics method in flux calculation for
solar radiation, we repeat the analysis of Liou (1973), who numerically
studied the discrete-ordinate method (DOM) for radiative transfer by comparing
the reflection and transmission values obtained via DOM with those obtained
using the doubling method.

The doubling method, first introduced by Van de Hulst (1963), is widely
considered one of the most accurate tools for multiple-scattering calculations
(Liu & Weng, 2006; Zhang & Li, 2013; Chou, 1992; Li & Ramaswamy, 1996; Ayash
et al., 2008). In this approach, it is assumed that the reflection and
transmission properties of a single, very thin, homogeneous layer are known.
The reflection and transmission properties of a layer twice as thick is then
found by considering two adjacent layers with properties described above and
summing the contributions from each layer (van de Hulst, 1980). The reflection
and transmission values for an arbitrarily thick layer is thus calculated by
repeatedly doubling until arriving at the desired thickness. Derived from the
work of Stokes (1862), the method computes the radiative properties of the
medium instead of solving for the radiance field explicitly (Evans & Stephens,
1991). This allows for the easy computation of radiance exiting the atmosphere
for many boundary conditions after the solution has been found. Although a
simple concept compared to the discrete ordinates and spherical harmonics
methods, the doubling method is very computational expensive due to the large
iteration loop over the layers, therefore it is not optimal for retrievals or
data assimilation. Instead, the doubling method is often used as an accurate
and detailed benchmark for comparison with fast, approximate radiative
transfer models.

We calculated the reflection and transmission values obtained by using the
Toon89 implementation in PICASO (Batalha et al., 2022) and 2-term and 4-term
spherical harmonics (SH2 & SH4). Following the analysis of Liou (1973), we
consider single-layer atmospheres of optical thickness $\tau_{N}=0.25,1,4$ and
$16$, where the scattering phase function is given by the Henyey-Greenstein
phase function (99) with asymmetry factor $g_{0}=0.75$. For single scattering
albedos we choose values that represent the potential values for visible and
near infrared of the solar spectrum for clouds and aerosols, $w_{0}=1$ and
$w_{0}=0.8$. Additionally, we choose to benchmark a range of cosine solar
zenith angles $\mu_{0}=0.1,0.5$ and 0.9. We note that for comparison with the
DOM and doubling methods discussed by Liou (1973), we are restricted to the
same parameter values considered by the authors. The reflection $r$, diffuse
transmission $t_{\text{dif}}$ and direct transmission $t_{\text{dir}}$ for
solar radiation are given by

$\displaystyle r$ $\displaystyle=\frac{F^{+}(0)}{\pi\mu_{0}F_{0}},$ (72)
$\displaystyle t_{\text{dif}}$
$\displaystyle=\frac{F^{-}(\tau_{N})}{\pi\mu_{0}F_{0}},$ (73) $\displaystyle
t_{\text{dir}}$ $\displaystyle=\exp(\tau_{N}/\mu_{0}),$ (74)

The total transmission $t$ is equal to $t_{\text{dif}}+t_{\text{dir}}$. Table
1 compares the reflection and total transmission values for each model to
those obtained using the doubling method which are taken from Table 1 in Liou
(1973).

In Figure 1 we plot the percentage difference between each model (Toon89, 2
and 4-term spherical harmonics [SH2, SH4], 2 and 4-stream discrete ordinates
[DOM2, DOM4]) and the doubling method. The values for DOM2, DOM4 and the
doubling method are taken from Liou (1973), where the results for the doubling
method are considered to be the “true” solution. We use orange lines for
4-stream models (SH4 and DOM4), pink for 2-stream models (SH2 and DOM2) with
the exception of Toon89, which is plotted in purple. The solid lines represent
the results produced by the authors of this work (SH and Toon89) whereas the
dashed lines represent those taken from (Liou, 1973) (DOM).

Figure 1 shows that the results produced by Toon89 and DOM2 are visibly
identical ($<0.01$%). This is not surprising since Toon89 is built upon the
Toon et al. (1989) 2-stream methodology which utilizes the quadrature
approximation when deriving the solution coefficients. The quadrature
approximation is analogous to the 2-stream discrete-ordinates approach to
solving the radiative transfer equation.

In most cases considered, SH4 agrees more closely with the accurate doubling
method than SH2, DOM2 and Toon89. This is expected due to the higher-order
approximation of the phase function. The 4-term advantage is most evident for
the $\mu_{0}=0.9$ tests (Figures 1(c), 1(f), 1(i) and 1(l)), where the
percentage difference between the 2-stream methods and the doubling method is
approximately two orders of magnitude greater than that of the 4-stream
methods.

However, in a few cases, particularly for small cosine zenith angle
$\mu_{0}=0.1$, the 2-stream methods incur a lower percentage difference with
the doubling method than their 4-stream counterparts for optically thin
layers. For example, both the reflection and transmission values produced by
SH2 in the $w_{0}=1,\mu_{0}=0.1$ case (Figures 1(a)&(g)), agree more closely
than SH4 with the doubling benchmark for optical depth less than around 0.6.
It is important to note, however, that in this case, the 4-stream DOM method
exhibits similar behaviour to SH4 and also suffers from worse agreement with
the doubling method than its 2-stream counterpart, indicating that this
behaviour in the optical thin region is not due to the choice of modelling
method (SH versus DOM) but rather that 4-stream techniques are still limited
in the improvement that can be achieved over 2-stream techniques in particular
cases.

In general, in the regions where SH2 exceeds the accuracy of the SH4, such as
in the case discussed above, both models remain within a 5% difference with
the doubling method. Therefore, even though SH4 might not be the “optimal”
model in these regions, it is still exhibiting excellent agreement with the
benchmark model. Notably, when considering the full range of parameters
considered in this analysis, the SH4 values for reflectivity and transmission
are consistently within 10% of the doubling method, whereas SH2 can reach
differences of close to 100%. This can be clearly realised in Figure 1(f).

Furthermore, SH4 performs comparatively to, if not better than, DOM4. We
notice that in cases such as Figures 1(c), for $w_{0}=1,\mu_{0}=0.9$ the
maximum percentage difference of both SH4 and DOM4 with the doubling method
occurs at the optically thin region at the top of the atmosphere. However, SH4
observes a percentage difference of approximately 4% versus a 40% difference
for DOM4. Similar discrepancies are observed for other cases throughout the
analysis. Liou (1973) concluded that the 4-stream DOM method may be of
adequate accuracy for studies of the flux distribution in the transfer of
solar irradiance through cloudy atmospheres. We thus extend this conclusion to
the 4-term spherical harmonics method.

###  5.2 CDISORT Comparison

The discrete ordinate solver, DISORT (Stamnes et al., 1988, 2000) is a
versatile, well-tested and one of the most widely used one-dimensional
radiative transfer solvers. The numerical capabilities of DISORT extend to
$N$-stream discrete ordinates approximations, where $N$ is arbitrary and
considerably greater than 4 (typically DISORT is run for 32 streams). DISORT
was originally written in FORTRAN (Stamnes et al., 2000). T. Dowling rewrote
the code in C, which we will refer to as CDISORT (Mayer & Kylling, 2005; Buras
et al., 2011).

To further study the efficacy of the spherical harmonics method of solving the
azimuthally-averaged radiative transfer equation (1), we benchmark our model
outputs against that of CDISORT.

We consider a test atmosphere that is comprised of 30 layers with cumulative
optical depth $\tau=10^{-2}-10^{2}$, with constant values for the single
scattering albedo $w_{0}$ and asymmetry parameter $g_{0}$, and use the Henyey-
Greenstein phase function (99) to describe scattering. We consider incident
solar radiation $F_{\odot}=1$ with a cosine solar zenith angle of $\mu_{0}$.

#### Azimuthally-averaged intensity with asymmetry

Figure 2: Azimuthally-averaged intensity emergent from the top of the
atmosphere with asymmetry $g_{0}$ for a variety of zenith angles $\mu,\mu_{0}$
and single scattering albedo $w_{0}=0.4,0.9$. We compare CDISORT, Toon89 and
SH4 and notice that the agreement worsens for decreasing cosine zenith angle.

In Figure 2 we plot the azimuthally-averaged intensity emergent from the top
of the atmosphere $I(0,\mu)$ against the asymmetry parameter $g_{0}$ for
32-stream CDISORT, Toon89 and 4-term spherical harmonics (SH4). This quantity
is important as it is used to calculate the geometric albedo of the planet.
Each panel of Figure 2 represents a different cosine zenith angle, namely
$\mu=0.9,0.7,0.4,0.2$. For each case, we have chosen to set $\mu$ to be equal
to the cosine solar zenith angle $\mu_{0}$. We also consider two single
scattering albedos $w_{0}=0.4$ and $w_{0}=0.9$ in each plot, similar to our
doubling method benchmarking exercise against Liou (1973). In each plot, we
notice that SH4 more closely matches 32-stream CDISORT than Toon89 does,
especially towards high asymmetry. However, overall the agreement between all
models worsens as $g_{0}$ increases. This is as expected because low-order
approximations are inadequate representations of highly asymmetric phase
functions. We also notice that the digression of the solutions for high
asymmetry is more pronounced for smaller values of $\mu$ and $\mu_{0}$. This
is similar to the behaviour observed in the Liou (1973) analysis in Section
5.1.

#### Azimuthally-averaged intensity with optical depth

Figure 3: Azimuthally-averaged intensity emergent from the top of the
atmosphere with optical depth $\tau$ for a variety of zenith angles
$\mu=\mu_{0}=0.2,0.4,0.7,0.9$ and asymmetry $g_{0}=0.4,0.9$. The single
scattering albedo is $w_{0}=0.9$. We compare CDISORT, Toon89 and SH4 and
notice that the agreement worsens for decreasing zenith angle and increasing
$g_{0}$.

In Figure 3 we study how the azimuthally-averaged intensity $I(\tau,\mu)$
varies with optical depth $\tau$ for each of the models. For this analysis, we
set $w_{0}=0.9$ and consider two values of asymmetry, namely $g_{0}=0.4$ and
$g_{0}=0.9$. We notice that the agreement between all three models is greater
for low asymmetry ($g_{0}=0.4$) – the models are more disparate for the high
asymmetry value of $g_{0}=0.9$. Similar to the analysis of Figure 2, the
agreement between the solutions worsens as $\mu=\mu_{0}$ decreases. This
effect is consistent with what was observed in the Liou analysis in Section
5.1. However, despite its limitations in high asymmetry and low cosine zenith
angles, SH4 does offer a significant improvement in agreement with CDISORT
over Toon89 for all cases.

#### Layer-wise vertical fluxes

Figure 4: Fluxes with optical depth $\tau$ for cosine solar zenith angle (a)
$\mu_{0}=0.9$ and (b) $\mu_{0}=0.2$, single scattering albedo $w_{0}=0.5$ and
asymmetry $g_{0}=0.0,0.9$. We compare CDISORT with Toon89 and SH4 and notice
that SH4 agrees marginally better with CDISORT than Toon89.

Figure 4 depicts the layer-wise downward, upward and net fluxes calculated by
each model, where Figure 4(a) is for incoming angle $\mu_{0}=0.9$ and 4(b) is
for $\mu_{0}=0.2$. We consider a constant single scattering albedo of
$w_{0}=0.5$ and include results for two extreme asymmetry values in each plot,
namely $g_{0}=0.0$ and $g_{0}=0.9$. We again notice an improvement in the
agreement with CDISORT of SH4 from Toon89. In particular, we observe negative
upward flux values for Toon89 in Figure 4(a). Negative fluxes have been
reported for delta-Eddington calculations in the literature and are discussed
in detail by (Wiscombe, 1977). The authors explain that, although unphysical,
a negative flux does is not particularly more problematic than an inaccurate
positive flux; the importance lies in how great the error is. SH4 does not
exhibit this negative-flux behaviour.

For each of the 30-layer, single-wavelength simulations studied in this
section, the Toon89 framework in PICASO took on average 4e-5 seconds to run,
while SH2 and SH4 took 2.5e-4 seconds and 3e-4 seconds respectively, using a
2021 Macbook Pro with the Apple M1 chip. All of the PICASO packages leverage
numba’s CPU just-in-time compiling. (Lam et al., 2015). We note that the
Toon89 framework has been extensively used and sufficiently optimized, whereas
the SH framework within PICASO is still in its infancy. With adequate
computational optimization and code refinement, we expect SH2 to operate at
similar speeds to Toon89, and for SH4 to exhibit a comparable increase in
speed. The equivalent 30-layer, single-wavelength simulations ran by 32-stream
CDISORT took on average 2.4e-1 seconds – three orders of magnitude slower than
the current implementation of SH4. Therefore, even without the appropriate
speed optimization and code refinement, SH4 offers significantly faster
solutions than the higher-order simulations of CDISORT, whilst allowing
notable improvement in model agreement over the two-stream Toon89.

###  5.3 Geometric Albedo Spectrum

Figure 5: Geometric albedo for a Jupiter atmosphere with (a) Jupiter clouds
and (b) single parameter box cloud profile. We compare the spectra produced by
Toon89 and SH4, where we consider two different scattering options for the
latter; one for which the single scattering phase function is taken to be two-
term Henyey-Greenstein (TTHG) while multi-scattering is one-term Henyey
Greenstein, and the other for which all scattering is assumed to be TTHG. The
former of these two SH4 options is most similar to the Toon89 analysis.

Finally, we compare the ultimate geometric albedo produced by Toon89 and SH4
for a Jupiter-like system (1$\times$Solar metallicty, 25 ms-2 gravity, 5 AU
from a sun-like star) with two different cloud parameterizations: 1) Jupiter-
like water clouds with a sedimentation efficiency (Ackerman & Marley, 2001) of
f${}_{\rm sed}=1$, and 2) a single-parameter box cloud profile described by a
single asymmetry and single scattering albedo. The former case was also used
in Batalha et al. (2019) as the benchmark reflected light cloudy case. For the
latter, we specify a cloud layer located between $P=1$ and $P=0.1$ bars with
asymmetry $g_{0}=0.9$, single scattering $w_{0}=0.8$ and layer optical depth
$\tau=0.5$. These cases allow us to sample two different scattering regimes.
The Jupiter-like water cloud represents the scenario where the cloud is
dominant source of scattering throughout the 0.3-1$\mu$m. The single-parameter
box cloud represents the scenario where Rayleigh is the dominant source of
scattering at 0.5$<\mu$m, and clouds are the dominant source 0.5$>\mu$m. This
is particularly important to test because of the different methodologies for
including Rayleigh scattering.

To model these two cases, we test two different scattering approaches for SH4
against the default scattering approach used in Toon89. In the Toon89
framework of PICASO, multiple-scattering for clouds is always modelled using a
one-term Henyey-Greenstein (OTHG) phase function (99), whereas single-
scattering is represented using a weighted sum of a two-term Henyey-Greenstein
(TTHG) phase function (103) for clouds, and Rayleigh, as described in Appendix
F.1. SH4 has more flexibility: namely, we can prescribe either OTHG or TTHG
for both single and multiple scattering for clouds, and we can include
Rayleigh scattering in either scattering regimes. For this comparison, we
replicate the Toon89 framework for single scattering, but consider two cases
for multiple scattering: (a) OTHG (with Rayleigh) and (b) TTHG (with
Rayleigh). The geometric albedo for each parameterization is plotted in Figure
5. The different modelling approaches produce evident differences in the
geometric albedo.

For the Jupiter-like-cloud case [Figure 5(a)], the SH4 model with an OTHG
multiple-scattering function agrees more closely with Toon89 than the TTHG
case, where the TTHG case is significantly brighter than the other two models.
The purpose of the TTHG framework is to better capture the back scattering
peak, which contributes to the overall brightness (Batalha et al., 2019).
Therefore, one major drawback of the standard PICASO methodology is its
inability to account for the back scattering peak for multiply scattered
photons. Figure 5(a) shows how SH4 can better account for this.

For the single parameter box cloud profile [Figure 5(b)], both SH4 cases
produce different spectra than Toon89. In the regime dominated by clouds
($>0.58\mu$m), the SH4 TTHG is still the brightest case, similar to the case
shown in Figure 5(a). This suggests that, similar to the Jupiter-like case,
the differences in this regime are due to the additional brightness from the
inclusion of the TTHG treatment. In the region dominated by Rayleigh
($<0.58\mu$m) the behaviour flips and Toon89 is brighter than both SH4 OTHG
and TTHG. Because the SH4 OTHG and TTHG spectra are similar in this regime,
the spectral differences between SH4 and Toon89 are due to the differences in
how the Rayleigh scattering is incorporated. PICASO post-processes Rayleigh
into the multiple- and single-scattering phase function after the Toon89
methodology. SH2/4 incorporate Rayleigh directly in the phase functions
expansion (see Appendix F.1).

##  6 Conclusion

In this work, we presented a spherical harmonics approach to solving the
radiative transfer equation for reflected light, which has been implemented in
modelling software PICASO Batalha et al. (2022). This new methodology allows
us to approximate scattering phase functions to higher orders than the
original, two-stream implementation in PICASO, which we denoted Toon89 to
reflect its heritage from Toon et al. (1989). The primary purpose of this work
was to rigorously outline the derivation of the model, explaining the steps of
the spherical harmonics analysis and explicitly stating the matrix systems
being solved by the model. To demonstrate the efficacy of the higher-order
model, we compared our results to two independent models from Liou (1973) and
Stamnes et al. (2000). These analyses illustrated the expected superiority of
higher-order approximations in radiative transfer calculations, whilst also
elucidating the extent of the improvement when moving from 2 to 4 term
expansions. We studied the vertical distribution of fluxes by comparing the
reflection and transmission values obtained via 2 and 4-term discrete
ordinates methods (DOM), as well as the more accurate doubling method with
those produced by 2 and 4-term spherical harmonics in Section 5.1. This
investigation highlighted that the choice of model (DOM or SH) is negligible
compared to the order of approximation (2 or 4), where the 4-term
approximations significantly out-performed the 2-term approximations. However,
in the majority of the cases studied in this analysis, the SH approach was
marginally better than the DOM approach. Extending the conclusion of Liou
(1973) regarding the comparison of low-order approximations with the doubling
method, we conclude that the 4-term spherical harmonics method may be of
adequate accuracy for studied of the flux distribution in the transfer of
solar irradiance through cloudy atmospheres. When comparing the 2-stream
Toon89 and four-term SH4 to 32-stream CDISORT in Section 5.2, we noticed a
significant enhancement of model accuracy when moving from 2 to 4 term
approximations; we achieve notable improvement with minimal cost to
computational expense.

Finally, we considered the impact of SH4 and the different scattering
capabilities on the geometric albedo for a Jupiter atmosphere with two types
of clouds: (a) Jupiter clouds and (b) single parameter box cloud profile. This
analysis illustrated that the choice of multiple-scattering function can
significantly affect the reflected light spectra, thus the user-specification
of scattering behaviour renders the SH4 model useful for retrievals and fine-
tuning of atmospheric models. A major drawback of the standard PICASO
methodology is its inability to account for the back scattering peak for
multiply scattered photons, a phenomenon that SH4 can better capture. This
analysis further highlighted the impact on spectra depending on how the
Rayleigh scattering is incorporated.

The spherical harmonics technique outlined in this paper is also appropriate
for modelling thermal emission. This will be the focus of a future
publication.

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.

##  Appendix A Spherical harmonics methodology

The azimuthally-averaged radiative transfer equation for reflected light is
given by (1). By substituting (2)–(3) into (1), we obtain

$\begin{split}\sum_{l=0}^{L}(2l+1)\mu\frac{\mathrm{d}I_{l}}{\mathrm{d}\tau}P_{l}(\mu)&=\sum_{l=0}^{L}(2l+1)I_{l}(\tau)P_{l}(\mu)\\\
&-\frac{w_{0}}{2}\int_{-1}^{1}\sum_{l=0}^{L}(2l+1)I_{l}(\tau)P_{l}(\mu^{\prime})\sum_{k=0}^{L}\chi_{k}P_{k}(\mu)P_{k}(\mu^{\prime})\mathrm{d}\mu^{\prime}-\sum_{l=0}^{L}b_{l}e^{-\frac{\tau}{\mu_{0}}}P_{l}(\mu),\end{split}$
(75)

where

$b_{l}=\frac{w_{0}\chi_{l}F_{\odot}}{4\pi}\chi_{l}P_{l}(-\mu_{0}),$ (76)

Using the recurrence relation and orthogonality principle of Legendre
polynomials, i.e.

$\displaystyle(2l+1)\mu P_{l}(\mu)$
$\displaystyle=lP_{l-1}(\mu)+(l+1)P_{l+1}(\mu),$ (77)
$\displaystyle\int_{-1}^{1}P_{l}(\mu)P_{k}(\mu)\mathrm{d}\mu$
$\displaystyle=\frac{2}{2l+1}\delta_{l,k},$ (78)

we obtain

$\sum_{l=0}^{L}\left[lP_{l-1}(\mu)+(l+1)P_{l+1}(\mu)\right]\frac{\mathrm{d}I_{l}}{\mathrm{d}\tau}=\sum_{l=0}^{L}[(2l+1-w_{0}\chi_{l})I_{l}(\tau)-b_{l}e^{-\frac{\tau}{\mu_{0}}}]P_{l}(\mu),$
(79)

which can be rewritten as

$\sum_{l=0}^{L}\left[(l+1)\frac{\mathrm{d}I_{l+1}}{\mathrm{d}\tau}+l\frac{\mathrm{d}I_{l-1}}{\mathrm{d}\tau}\right]P_{l}(\mu)=\sum_{l=0}^{L}[(2l+1-w_{0}\chi_{l})I_{l}(\tau)-b_{l}e^{-\frac{\tau}{\mu_{0}}}]P_{l}(\mu),$
(80)

##  Appendix B $P_{1}$ single-layer solution

The $P_{1}$ single-layer problem (8) is linear in intensity, so is
traditionally solved by combining the “homogeneous” and “particular”
solutions, where the former must satisfy (8) with no source term whereas the
latter is a solution to (8) with the source term included. We first consider
the homogeneous problem, namely

$\frac{\mathrm{d}}{\mathrm{d}\tau}\begin{pmatrix}I_{0}\\\
I_{1}\end{pmatrix}=\begin{pmatrix}0&a_{1}\\\
a_{0}&0\end{pmatrix}\begin{pmatrix}I_{0}\\\ I_{1}\end{pmatrix}.$ (81)

We assume that this homogeneous solution is of the form
$I_{l}=G_{l}e^{\lambda\tau}$ for $l=0,1$. Substituting this general form into
(81) yields

$\begin{pmatrix}-\lambda&a_{1}\\\
a_{0}&-\lambda\end{pmatrix}\begin{pmatrix}G_{0}\\\ G_{1}\end{pmatrix}=0.$ (82)

To ensure that (82) has a solution, the value of the determinant must be zero,
namely

$\displaystyle\lambda^{2}-a_{0}a_{1}=0\Rightarrow\lambda=\pm\sqrt{a_{0}a_{1}}.$
(83)

The coefficients of $G_{l}$ are not independent, rather one can define

$G_{1}=\frac{\lambda}{a_{1}}G_{0}.$ (84)

For the particular solution, we assume the form
$I_{l}=\eta_{l}e^{-\tau/\mu_{0}}$ for $l=0,1$. Substituting the particular
solution into (8), we can obtain the coefficients of $\zeta_{l}$ by solving

$\begin{pmatrix}1/\mu_{0}&a_{1}\\\
a_{0}&1\mu_{0}\end{pmatrix}\begin{pmatrix}\eta_{0}\\\
\eta_{1}\end{pmatrix}=\begin{pmatrix}b_{1}\\\ b_{0}\end{pmatrix}.$ (85)

By inverting the matrix in (85), we obtain the particular solution

$\displaystyle\eta_{0}=\frac{1}{1/\mu_{0}^{2}-a_{0}a_{1}}\left(\frac{b_{1}}{\mu_{0}}-a_{1}b_{0}\right),$
(86)
$\displaystyle\eta_{1}=\frac{1}{1/\mu_{0}^{2}-a_{0}a_{1}}\left(\frac{b_{0}}{\mu_{0}}-a_{0}b_{1}\right).$

Thus, we can write the full solution to (8):

$\begin{pmatrix}I_{0}\\\
I_{1}\end{pmatrix}=\begin{pmatrix}e^{-\lambda\tau}&e^{\lambda\tau}\\\
-qe^{-\lambda\tau}&qe^{\lambda\tau}\end{pmatrix}\begin{pmatrix}X_{0}\\\
X_{1}\end{pmatrix}+\begin{pmatrix}\eta_{0}\\\
\eta_{1}\end{pmatrix}e^{-\frac{\tau}{\mu_{0}}},$ (87)

where $q=\lambda/a_{1}$ and $\lambda=\sqrt{a_{0}a_{1}}$.

##  Appendix C $P_{3}$ single-layer solution

To solve the $P_{3}$ single-layer problem (17), we first consider the
homogeneous problem, namely

$\frac{\mathrm{d}}{\mathrm{d}\tau}\begin{pmatrix}I_{0}\\\ I_{1}\\\ I_{2}\\\
I_{3}\end{pmatrix}=\begin{pmatrix}0&a_{1}&0&-\frac{2a_{3}}{3}\\\
a_{0}&0&0&0\\\ 0&0&0&\frac{a_{3}}{3}\\\
-\frac{2a_{0}}{3}&0&\frac{a_{2}}{3}&0\\\ \end{pmatrix}\begin{pmatrix}I_{0}\\\
I_{1}\\\ I_{2}\\\ I_{3}\end{pmatrix},$ (88)

As with the $P_{1}$ case, we assume the solution takes the form
$I_{l}=G_{l}e^{\lambda\tau}$ for $l=0,1,2,3$. Substituting this general form
into (88) yields

$\begin{pmatrix}-\lambda&a_{1}&0&-\frac{2a_{3}}{3}\\\ a_{0}&-\lambda&0&0\\\
0&0&-\lambda&\frac{a_{3}}{3}\\\
-\frac{2a_{0}}{3}&0&\frac{a_{2}}{3}&-\lambda\end{pmatrix}\begin{pmatrix}G_{0}\\\
G_{1}\\\ G_{2}\\\ G_{3}\end{pmatrix}=0.$ (89)

To ensure that (89) has a solution, the value of the determinant must be zero,
namely

$\displaystyle\lambda^{4}-\beta\lambda^{2}+\gamma=0\Rightarrow\lambda_{1,2}=\sqrt{\frac{\beta\pm\sqrt{\beta^{2}-4\gamma}}{2}},$
(90)

where

$\displaystyle\beta=a_{0}a_{1}+\frac{1}{9}a_{2}a_{3}+\frac{4}{9}a_{0}a_{3},\qquad\gamma=\frac{1}{9}a_{0}a_{1}a_{2}a_{3}.$
(91)

The coefficients of $G_{l}$ are not independent, rather by considering (88)
they are related by

$\displaystyle\lambda G_{0}$ $\displaystyle=a_{1}G_{1}-\frac{2a_{3}}{3}G_{3},$
(92) $\displaystyle\lambda G_{1}$ $\displaystyle=a_{0}G_{0},$
$\displaystyle\lambda G_{2}$ $\displaystyle=\frac{a_{3}}{3}G_{3},$
$\displaystyle\lambda G_{3}$
$\displaystyle=-\frac{2a_{0}}{3}G_{0}+\frac{a_{2}}{3}G_{2},$

and $G_{l}$ for $l=1,2,3$ can be written in terms of $G_{0}$:

$\displaystyle G_{1}$ $\displaystyle=\frac{a_{0}}{\lambda}G_{0},$ (93)
$\displaystyle G_{2}$
$\displaystyle=\frac{1}{2}\left(\frac{a_{0}a_{1}}{\lambda^{2}}-1\right)G_{0},$
$\displaystyle G_{3}$
$\displaystyle=\frac{3}{2a_{3}}\left(\frac{a_{0}a_{1}}{\lambda}-\lambda\right)G_{0}.$

For the particular solution, we assume the form
$I_{l}=\eta_{l}e^{-\tau/\mu_{0}}$ for $l=0,1,2,3$. Substituting the particular
solution into (8), we can obtain the coefficients of $\eta_{l}$ by solving

$\begin{pmatrix}1/\mu_{0}&a_{1}&0&-\frac{2a_{3}}{3}\\\ a_{0}&1/\mu_{0}&0&0\\\
0&0&1/\mu_{0}&\frac{a_{3}}{3}\\\
-\frac{2a_{0}}{3}&0&\frac{a_{2}}{3}&1/\mu_{0}\end{pmatrix}\begin{pmatrix}\eta_{0}\\\
\eta_{1}\\\ \eta_{2}\\\
\eta_{3}\end{pmatrix}=\begin{pmatrix}-\frac{2b_{3}}{3}+b_{1}\\\ b_{0}\\\
\frac{b_{3}}{3}\\\ -\frac{2b_{0}}{3}+\frac{b_{2}}{3}\end{pmatrix}.$ (94)

By inverting the matrix in (94), we obtain $\eta_{l}=\Delta_{l}/\Delta$ where

$\displaystyle\Delta_{0}$
$\displaystyle=(a_{1}b_{0}-b_{1}/\mu_{0})(a_{2}a_{3}-9/\mu_{0}^{2})+2(a_{3}b_{2}-2a_{3}b_{0}-3b_{3}/\mu_{0})/\mu_{0}^{2},$
(95) $\displaystyle\Delta_{1}$
$\displaystyle=(a_{0}b_{1}-b_{0}/\mu_{0})(a_{2}a_{3}-9/\mu_{0}^{2})-2a_{0}(a_{3}b_{2}-3b_{3}/\mu_{0})/\mu_{0},$
$\displaystyle\Delta_{2}$
$\displaystyle=(a_{3}b_{2}-3b_{3}/\mu_{0})(a_{0}a_{1}-1/\mu_{0}^{2})-2a_{3}(a_{0}b_{1}-b_{0}/\mu_{0})/\mu_{0},$
$\displaystyle\Delta_{3}$
$\displaystyle=(a_{2}b_{3}-3b_{2}/\mu_{0})(a_{0}a_{1}-1/\mu_{0}^{2})+2(3a_{0}b_{1}-2a_{0}b_{3}-3b_{0}/\mu_{0})/\mu_{0}^{2},$

and

$\displaystyle\Delta=9f(1/\mu_{0})\qquad\text{where}\qquad f(x)=x^{4}-\beta
x^{2}+\gamma.$ (96)

Thus, we can write the full solution to (17):

$\begin{pmatrix}I_{0}\\\ I_{1}\\\ I_{2}\\\
I_{3}\end{pmatrix}=\begin{pmatrix}e^{-\lambda_{1}\tau}&e^{\lambda_{1}\tau}&e^{-\lambda_{2}\tau}&e^{\lambda_{2}\tau}\\\
R_{1}e^{-\lambda_{1}\tau}&-R_{1}e^{\lambda_{1}\tau}&R_{2}e^{-\lambda_{2}\tau}&-R_{2}e^{\lambda_{2}\tau}\\\
Q_{1}e^{-\lambda_{1}\tau}&Q_{1}e^{\lambda_{1}\tau}&Q_{2}e^{-\lambda_{2}\tau}&S_{2}e^{\lambda_{2}\tau}\\\
S_{1}e^{-\lambda_{1}\tau}&-S_{1}e^{\lambda_{1}\tau}&S_{2}e^{-\lambda_{2}\tau}&-S_{2}e^{\lambda_{2}\tau}\end{pmatrix}\begin{pmatrix}X_{0}\\\
X_{1}\\\ X_{2}\\\ X_{3}\end{pmatrix}+\begin{pmatrix}\eta_{0}\\\ \eta_{1}\\\
\eta_{2}\\\ \eta_{3}\end{pmatrix}e^{-\frac{\tau}{\mu_{0}}},$ (97)

where

$\displaystyle R_{1,2}$ $\displaystyle=-\frac{a_{0}}{\lambda_{1,2}},\qquad
Q_{1,2}=\frac{1}{2}\left(\frac{a_{0}a_{1}}{\lambda_{1,2}^{2}}-1\right),\qquad
S_{1,2}=-\frac{3}{2a_{3}}\left(\frac{a_{0}a_{1}}{\lambda_{1,2}}-\lambda_{1,2}\right).$
(98)

##  Appendix D Henyey-Greeinstein phase function

Many of the analyses discussed in this work use the Henyey-Greenstein phase
function (Henyey & Greenstein, 1941):

$\mathcal{P}_{\text{HG}}(\cos\Theta)=\frac{1-g_{0}^{2}}{(1+g_{0}^{2}-2g_{0}\cos\Theta)^{3/2}},$
(99)

where the scattering angle $\Theta$ is defined as

$\displaystyle\cos\Theta$
$\displaystyle=\mu\mu^{\prime}-\sqrt{1-\mu^{2}}\sqrt{1-\mu^{\prime
2}}\cos(\phi-\phi^{\prime}).$ (100)

The direction of incident and outgoing scattered radiation is defined by the
cosine of the zenith angles, denoted $\mu^{\prime}$ and $\mu$ respectively,
and the azimuth angles, $\phi^{\prime}$ and $\phi$. The asymmetry parameter is
denoted $g_{0}$. The azimuthally-averaged Henyey-Greenstein function can be
expressed in terms of Legendre polynomials (Liou, 2002):

$\mathcal{P}(\mu,\mu^{\prime})=\sum_{l=0}^{N}\chi_{l}P_{l}(\mu)P_{l}(\mu^{\prime}),$
(101)

for moments

$\chi_{l}=(2l+1)g_{0}^{l}.$ (102)

The PICASO (Batalha et al., 2022) methodology also utilises a two-term Henyey-
Greenstein phase function in an attempt to capture back-scattering. The two-
term Henyey-Greenstein (TTHG) phase function is given by

$\mathcal{P}_{\text{TTHG}}(\cos\Theta)=\alpha\frac{1-g_{1}^{2}}{(1+g_{1}^{2}-2g_{1}\cos\Theta)^{3/2}}+(1-\alpha)\frac{1-g_{2}^{2}}{(1+g_{2}^{2}-2g_{2}\cos\Theta)^{3/2}},$
(103)

where we have two asymmetry factors, $g_{1}$ for forward scattering and
$g_{2}$ for backward scattering, and new parameter $\alpha$ to determine the
fraction of forward to back scattering:

$\alpha=C_{1}+C_{2}g_{2}^{C_{3}}.$ (104)

Users can specify $C_{1}$, $C_{2}$, and $C_{3}$, however, PICASO defaults
these to $C_{1}=1$, $C_{2}=-1$, and $C_{3}=2$. The moments for TTHG are given
by

$\chi_{l}=(2l+1)\left[\alpha g_{1}^{l}+(1-\alpha)g_{2}^{l}\right].$ (105)

##  Appendix E $\delta-$M approximation

Low order Legendre expansions are not capable of accurately representing
highly forward scattering phase functions. High scattering asymmetries are
produced by Mie scattering particles with sizes larger than typical optical
wavelengths, thus we must consider this limitation of the Legendre
methodology. To combat this, we apply the $\delta-$M technique (Joseph et al.,
1976; Wiscombe, 1977), which involves approximating the phase function by a
Dirac delta function forward scatter peak and an M-term expansion of the phase
function. Following the derivations of Joseph et al. (1976); Wiscombe (1977);
Cuzzi et al. (1982), we outline the $\delta-$M methodology for our problem.

Consider the radiative transfer equation (1) without the source function:

$\mu\frac{\partial
I}{\partial\tau}(\tau,\mu)=I(\tau,\mu)-\frac{w_{0}}{2}\int_{-1}^{1}I(\tau,\mu^{\prime})\mathcal{P}(\mu,\mu^{\prime})\mathrm{d}\mu^{\prime}.$
(106)

The azimuthally-averaged phase function $\mathcal{P}(\mu,\mu^{\prime})$ can be
expanded in Legendre polynomials $P_{l}$ (Chandrasekhar, 1960)

$\displaystyle\mathcal{P}(\mu,\mu^{\prime})$
$\displaystyle=\mathcal{P}(\cos\Theta)=\sum_{l=0}^{M-1}\chi_{l}P_{l}(\mu)P_{l}(\mu^{\prime}),$
(107) $\displaystyle\chi_{l}$
$\displaystyle=\frac{2l+1}{2}\int_{-1}^{1}\mathcal{P}(\cos\Theta)P_{l}(\cos\Theta)\,\mathrm{d}\cos\Theta,$
(108)

Approximating the phase function by a Dirac delta function forward scatter
peak and an M-term expansion, we obtain

$\displaystyle\mathcal{P}_{\delta}(\mu,\mu^{\prime})$
$\displaystyle=2f\delta(\mu-\mu^{\prime})+(1-f)\sum_{l=0}^{M-1}\chi_{l}^{*}P_{l}(\mu)P_{l}(\mu^{\prime}),$
(109)
$\displaystyle=2f\delta(\mu-\mu^{\prime})+(1-f)\mathcal{P}^{*}(\mu,\mu^{\prime}),$
(110)

where $f$ is the fractional scattering into the forward peak. If we substitute
$\mathcal{P}_{\delta}(\mu,\mu^{\prime})$ (110) in place of the azimuthally-
averaged phase function $\mathcal{P}(\mu,\mu^{\prime})$ in (106) we obtain

$\mu\frac{\partial
I}{\partial\tau}(\tau,\mu)=(1-w_{0}f)I(\tau,\mu)-\frac{(1-f)w_{0}}{2}\int_{-1}^{1}I(\tau,\mu^{\prime})\mathcal{P}^{*}(\mu,\mu^{\prime})\mathrm{d}\mu^{\prime},$
(111)

which we rewrite as

$\mu\frac{\partial
I}{\partial\tau^{*}}(\tau^{*},\mu)=I(\tau^{*},\mu)-\frac{w_{0}^{*}}{2}\int_{-1}^{1}I(\tau^{*},\mu^{\prime})\mathcal{P}^{*}(\mu,\mu^{\prime})\mathrm{d}\mu^{\prime},$
(112)

where

$\tau^{*}=(1-w_{0}f)\tau,\qquad w_{0}^{*}=\frac{(1-f)w_{0}}{1-w_{0}f}.$ (113)

Note that we require the moments of $\mathcal{P}_{\delta}(\mu,\mu^{\prime})$
be identical to that of the original phase function
$\mathcal{P}(\mu,\mu^{\prime})$ which allows us to determine the coefficients
$\chi_{l}^{*}$ of $\mathcal{P}^{*}(\mu,\mu^{\prime})$ to be

$\displaystyle\chi_{l}^{*}=\frac{\chi_{l}-(2l+1)f}{1-f}.$ (114)

The fractional scattering coefficient $f$ is evaluated by ensuring that the
$M^{\text{th}}$-order coefficient $\chi^{*}_{M}$ in the new phase function
$\mathcal{P}^{*}(\mu,\mu)$ is equal to zero, thus

$f=\frac{\chi_{M}}{2M+1}.$ (115)

In PICASO (Batalha et al., 2022), we primarily use either one-term Henyey-
Greenstein (OTHG) phase functions or two-term Henyey-Greenstein (TTHG) phase
function. For OTHG given by (99), the fractional scattering coefficient is

$f=g_{0}^{M},$ (116)

for asymmetry $g_{0}$ .

##  Appendix F Scattering capabilities

The spherical harmonics method implemented within the PICASO framework
(Batalha et al., 2022) has user-specific scattering functionality; namely, the
user can specify the type of phase function for single and multiple scattering
separately. The three types of phase functions currently offered by SH4 are
one-term Henyey-Greenstein (OTHG), two-term Henyey-Greenstein (TTHG) and
Rayleigh. The choice of phase function impacts the coefficients $\chi_{l}$ in
the Legendre approximation, which we recall is given by

$\mathcal{P}(\mu,\mu^{\prime})=\sum_{l=0}^{L}\chi_{l}P_{l}(\mu)P_{l}(\mu^{\prime}).$
(117)

For cloud-scattering captured by a OTHG phase function (99) with asymmetry
$g_{0}$, the coefficients are given by

$\chi_{l}^{\text{cld}}=(2l+1)g_{0}^{l},$ (118)

whereas for a TTHG phase function (103), we have

$\chi_{l}^{\text{cld}}=(2l+1)\left[\alpha
g_{1}^{l}+(1-\alpha)g_{2}^{l}\right].$ (119)

We note that for OTHG and TTHG we apply the $\delta-$M approximation described
in Appendix E, thus the coefficients (118) and (119) must be adjusted using
(114) and (115).

###  F.1 Rayleigh scattering

We include Rayleigh scattering by constructing a phase function that is
comprised of a weighted sum of the cloud-scattering and Rayleigh properties
(Batalha et al., 2019)

$\mathcal{P}(\mu,\mu^{\prime})=\frac{\tau_{\text{cld}}}{\tau_{\text{scat}}}\mathcal{P}_{\text{cld}}(\mu,\mu^{\prime})+\frac{\tau_{\text{ray}}}{\tau_{\text{scat}}}\mathcal{P}_{\text{ray}}(\mu,\mu^{\prime}),$
(120)

where $\tau_{\text{cld}}$ and $\tau_{\text{ray}}$ denote the opacity
contributed from the clouds and Rayleigh scattering respectively, and
$\tau_{\text{scat}}$ is the total scattering opacity. The Rayleigh scattering
phase function is taken to be

$\mathcal{P}_{\text{ray}}(\cos\Theta)=\frac{3}{4}(1+\cos^{2}\Theta),$ (121)

and the cloud-scattering phase function
$\mathcal{P}_{\text{cld}}(\mu,\mu^{\prime})$ will either be the OTHG (99) or
TTHG phase function (103). We note that before constructing this weighted sum,
the cloud-scattering phase function has already been adjusted using the
$\delta-$M approximation discussed in Appendix E. Let us denote

$f_{\tau_{\text{cld}}}=\frac{\tau_{\text{cld}}}{\tau_{\text{scat}}},\qquad
f_{\tau_{\text{ray}}}=\frac{\tau_{\text{ray}}}{\tau_{\text{scat}}}.$ (122)

Thus, to conduct the spherical harmonics analysis outlined in this paper, we
must be able to approximate the phase function in terms of Legendre
polynomials; that is, we must be able to write the phase function in the form

$\mathcal{P}(\mu,\mu^{\prime})=\sum_{l=0}^{L}\chi_{l}P_{l}(\mu)P_{l}(\mu^{\prime}).$
(123)

The coefficients $\chi_{l}$ for the weighted phase function (120) are

$\chi_{l}=f_{\tau_{\text{cld}}}\chi_{l}^{\text{cld}}+f_{\tau_{\text{ray}}}\chi_{l}^{\text{ray}}.$
(124)

Since the cloud-scattering phase function has already been adjusted using the
$\delta-$M approximation, the Legendre coefficients $\chi_{l}^{\text{cld}}$
are given by (114). The Rayleigh phase function (121) can be written in terms
of Legendre polynomials as

$\mathcal{P}_{\text{ray}}(\cos\Theta)=1+\frac{1}{2}P_{2}(\mu)P_{2}(\mu^{\prime}),$
(125)

where $P_{2}$ is the third Legendre polynomial. Therefore, the Legendre
coefficients $\chi_{l}^{\text{ray}}$ for the Rayleigh phase function are

$\chi_{l}^{\text{ray}}=\begin{cases}1,&\quad l=0,\\\ 0.5,&\quad l=2,\\\
0,&\quad\text{otherwise}.\end{cases}$ (126)

We notice that for 2-term expansions the Rayleigh contribution to the
scattering function is neglected, as it only appears in the third term ($l=2$)
of the phase function expansion. This is an issue in the Toon89 framework
within PICASO, but is managed by including the Rayleigh contribution by
forcing a third moment equal to $0.5f_{\tau_{\text{ray}}}$, so that when
Rayleigh dominates the total opacity, the correct scattering moments are
included. This is discussed in more detail in Batalha et al. (2019). No such
alternative measures are needed for 4-term spherical harmonics as the critical
$l=2$ contribution appears naturally in the phase function expansion.

##  Appendix G Derivation of boundary condition

To derive the boundary condition for upward flux $F^{+}(\tau_{N})$ at a
Lambertian surface, we must determine the total amount of flux being reflected
off of the surface. This will be the sum of the reflected diffuse downward
flux $F^{-}(\tau_{N})$ and the reflected direct flux from the solar beam. The
intensity of the solar beam at the surface is

$I_{direct}(\tau_{N},\mu)=\frac{F_{\odot}}{\pi}\mu_{0}e^{\frac{-\tau_{N}}{\mu_{0}}},$
(127)

from which we can derive the direct flux to be

$F^{-}_{direct}(\tau_{N})=2\pi\int_{0}^{-1}I_{direct}(\tau_{N},\mu)\mu\mathrm{d}\mu=F_{\odot}\mu_{0}e^{\frac{-\tau_{N}}{\mu_{0}}}.$
(128)

Thus, for a Lambertian surface with surface reflectivity $A_{S}$, the upward
flux $F^{+}(\tau_{N})$ is given by the reflected portion of the downward
diffuse and direct fluxes:

$F^{+}(\tau_{N})=A_{S}[F^{-}(\tau_{N})+\mu_{0}F_{\odot}e^{-\frac{\tau}{\mu_{0}}}].$
(129)

We can derive the equivalent boundary condition for $f^{+}(\tau_{N})$ by first
calculating $f^{-}_{direct}(\tau_{N})$ using

$f^{-}_{direct}(\tau_{N})=2\pi\int_{0}^{-1}I_{direct}(\tau_{N},\mu)\frac{1}{2}(5\mu^{3}-3\mu)\mathrm{d}\mu=-\frac{1}{4}F_{\odot}\mu_{0}e^{\frac{-\tau_{N}}{\mu_{0}}}.$
(130)

Similarly, for a Lambertian surface with surface reflectivity $A_{S}$, the
upward $f^{+}(\tau_{N})$ is given by

$f^{+}(\tau_{N})=A_{S}[f^{-}(\tau_{N})-\frac{1}{4}\mu_{0}F_{\odot}e^{-\frac{\tau}{\mu_{0}}}].$
(131)

##  Appendix H Modelling recommendations

In Section 5.3 we investigated the geometric albedo produced by the Toon89 and
SH4 frameworks within PICASO. Here, we aggregate our modelling recommendations
in a single table. In the current version of PICASO, these represent the
radiative transfer defaults. All of these toggles are controlled through the
PICASO approx picaso.justdoit.iputs.approx function. A tutorial on the use of
this function is available as Jupyter notebook.

  * •

Single Scattering

    * –

Use default specification for direct scattering:

      * *

Toon89: single_phase=‘TTHG_ray’.

      * *

SH4: psingle_form=‘TTHG’ with psingle_rayleigh=‘on’, and w_single_form=‘TTHG’
with w_single_rayleigh=‘on’.

    * –

Specify the functional form of the fraction, f, of forward and backward
scattering according to the problem being addressed.

  * •

Multiple Scattering

    * –

For planet cases with some degree of asymmetric clouds scatterers:

      * *

Toon89: always use N=2 Legendre polynomial expansion with the
$\delta-$Eddington approximation.

      * *

SH4: choose either w_multi_form=‘TTHG’ if you expect the back-scattering
fraction to be non-negligible, or ‘OTHG’ if you wish to primarily capture
forward scattering. Both phase functions should be adjusted with the
$\delta-$Eddington approximation.

    * –

Additionally for SH4, include Rayleigh effects in multiple scattering by
setting w_multi_rayleigh=‘on’.

## References

  * Ackerman & Marley (2001) Ackerman, A. S., & Marley, M. S. 2001, The Astrophysical Journal, 556, 872
  * Ayash et al. (2008) Ayash, T., Gong, S., & Jia, C. Q. 2008, Journal of the atmospheric sciences, 65, 2448
  * Barstow et al. (2020) Barstow, J. K., Changeat, Q., Garland, R., et al. 2020, MNRAS, 493, 4884, doi: 10.1093/mnras/staa548
  * Batalha et al. (2020) Batalha, N., caoimherooney11, & sagnickm. 2020, natashabatalha/virga: Initial Release, v0.0, Zenodo, doi: 10.5281/zenodo.3759888
  * Batalha et al. (2022) Batalha, N., Rooney, C., Blanch, N. R., & MacDonald, R. 2022, natashabatalha/picaso: Release 2.3, v2.3.0, Zenodo, doi: 10.5281/zenodo.6419943
  * Batalha et al. (2019) Batalha, N. E., Marley, M. S., Lewis, N. K., & Fortney, J. J. 2019, The Astrophysical Journal, 878, 70
  * Bokeh Development Team (2014) Bokeh Development Team. 2014, Bokeh: Python library for interactive visualization.  http://www.bokeh.pydata.org
  * Buras et al. (2011) Buras, R., Dowling, T., & Emde, C. 2011, Journal of Quantitative Spectroscopy and Radiative Transfer, 112, 2028
  * Chandrasekhar (1960) Chandrasekhar, S. 1960, Chandrasekhar, S, 20, 2
  * Chou (1992) Chou, M.-D. 1992, Journal of Atmospheric Sciences, 49, 762
  * Cuzzi et al. (1982) Cuzzi, J. N., Ackerman, T. P., & Helmle, L. C. 1982, Journal of the Atmospheric Sciences, 39, 917
  * Evans & Stephens (1991) Evans, K. F., & Stephens, G. 1991, Journal of Quantitative Spectroscopy and Radiative Transfer, 46, 413
  * Fiveland & Jessee (1996) Fiveland, V., & Jessee, J. 1996, Journal of thermophysics and heat transfer, 10, 445
  * Ge et al. (2015) Ge, W., Marquez, R., Modest, M. F., & Roy, S. P. 2015, Journal of Heat Transfer, 137
  * Heng et al. (2018) Heng, K., Malik, M., & Kitzmann, D. 2018, The Astrophysical Journal Supplement Series, 237, 29
  * Heng et al. (2014) Heng, K., Mendonça, J. M., & Lee, J.-M. 2014, The Astrophysical Journal Supplement Series, 215, 4
  * Henyey & Greenstein (1941) Henyey, L. G., & Greenstein, J. L. 1941, The Astrophysical Journal, 93, 70
  * Irvine (1975) Irvine, W. M. 1975, Icarus, 25, 175
  * Iwabuchi (2006) Iwabuchi, H. 2006, Journal of the atmospheric sciences, 63, 2324
  * Joseph et al. (1976) Joseph, J. H., Wiscombe, W., & Weinman, J. 1976, Journal of Atmospheric Sciences, 33, 2452
  * King (1986) King, M. D. 1986, Journal of Atmospheric Sciences, 43, 784
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in ELPUB, 87–90
  * Lam et al. (2015) Lam, S. K., Pitrou, A., & Seibert, S. 2015, in Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC, LLVM ’15 (New York, NY, USA: ACM), 7:1–7:6, doi: 10.1145/2833157.2833162
  * Lewis & Miller (1984) Lewis, E. E., & Miller, W. F. 1984
  * Li & Ramaswamy (1996) Li, J., & Ramaswamy, V. 1996, Journal of the atmospheric sciences, 53, 1174
  * Line et al. (2012) Line, M. R., Zhang, X., Vasisht, G., et al. 2012, ApJ, 749, 93, doi: 10.1088/0004-637X/749/1/93
  * Liou (1973) Liou, K.-N. 1973, Journal of Atmospheric Sciences, 30, 1303
  * Liou (1974) Liou, K.-n. 1974, J. Atmos. Sci, 31, 1473
  * Liou (2002) Liou, K.-N. 2002, An introduction to atmospheric radiation (Elsevier)
  * Liou et al. (1988) Liou, K.-N., Fu, Q., & Ackerman, T. P. 1988, Journal of Atmospheric Sciences, 45, 1940
  * Liu & Weng (2006) Liu, Q., & Weng, F. 2006, Journal of the Atmospheric Sciences, 63, 3459
  * Madhusudhan & Seager (2009) Madhusudhan, N., & Seager, S. 2009, ApJ, 707, 24, doi: 10.1088/0004-637X/707/1/24
  * MATLAB (2010) MATLAB. 2010, version 7.10.0 (R2010a) (Natick, Massachusetts: The MathWorks Inc.)
  * Mayer (2009) Mayer, B. 2009, in EPJ web of Conferences, Vol. 1, EDP Sciences, 75–99
  * Mayer & Kylling (2005) Mayer, B., & Kylling, A. 2005, Atmospheric Chemistry and Physics, 5, 1855
  * McKinney (2010) McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference, ed. S. van der Walt & J. Millman, 51 – 56
  * Meador & Weaver (1980) Meador, W., & Weaver, W. 1980, Journal of Atmospheric Sciences, 37, 630
  * Mihalas & Mihalas (2013) Mihalas, D., & Mihalas, B. W. 2013, Foundations of radiation hydrodynamics (Courier Corporation)
  * Modest (1989) Modest, M. F. 1989, Journal of Thermophysics and Heat Transfer, 3, 283
  * Modest (2013) —. 2013, Radiative Heat Transfer (Academic Press)
  * Olfe (1967) Olfe, D. 1967, AIAA Journal, 5, 638
  * Pérez & Granger (2007) Pérez, F., & Granger, B. E. 2007, Computing in Science & Engineering, 9
  * Ravishankar (2009) Ravishankar, M. 2009, PhD thesis, The Ohio State University
  * Rooney et al. (2022) Rooney, C. M., Batalha, N. E., Gao, P., & Marley, M. S. 2022, The Astrophysical Journal, 925, 33
  * Shibata & Uchiyama (1992) Shibata, K., & Uchiyama, A. 1992, Journal of the Meteorological Society of Japan. Ser. II, 70, 1097
  * Stamnes et al. (1988) Stamnes, K., Tsay, S.-C., Wiscombe, W., & Jayaweera, K. 1988, Applied optics, 27, 2502
  * Stamnes et al. (2000) Stamnes, K., Tsay, S.-C., Wiscombe, W., & Laszlo, I. 2000
  * Stephens & Preisendorfer (1984) Stephens, G. L., & Preisendorfer, R. W. 1984, Journal of the atmospheric sciences, 41, 725
  * Stokes (1862) Stokes, G. G. 1862, Proceedings of the Royal Society of London, 545
  * Thomas & Stamnes (2002) Thomas, G. E., & Stamnes, K. 2002, Radiative transfer in the atmosphere and ocean (Cambridge University Press)
  * Toon et al. (1989) Toon, O. B., McKay, C., Ackerman, T., & Santhanam, K. 1989, Journal of Geophysical Research: Atmospheres, 94, 16287
  * Van de Hulst (1963) Van de Hulst, H. 1963, Rept, Inst. Space Studies, NASA, New York
  * van de Hulst (1980) van de Hulst, H. C. 1980, Multiple light scattering, Elsevier
  * van Wijngaarden & Happer (2022) van Wijngaarden, W. A., & Happer, W. 2022, 2n-Stream Radiative Transfer, arXiv, doi: 10.48550/ARXIV.2205.09713
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, doi: 10.1038/s41592-019-0686-2
  * Walt et al. (2011) Walt, S. v. d., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science & Engineering, 13, 22
  * Wiscombe (1977) Wiscombe, W. 1977, Journal of Atmospheric Sciences, 34, 1408
  * Wiscombe & Joseph (1977) Wiscombe, W., & Joseph, J. 1977, Icarus, 32, 362
  * Zhang & Li (2013) Zhang, F., & Li, J. 2013, Journal of the atmospheric sciences, 70, 3084