Testing Approximate Infrared Scattering Radiative-transfer Methods for Hot Jupiter Atmospheres

The calculation of internal atmospheric (longwave) fluxes is a key component of any model of exoplanet atmospheres that requires radiative-transfer (RT) calculations. For atmospheres containing a strong scattering component such as cloud particles, most 1D multiple-scattering RT methods typically involve numerically expensive matrix inversions. This computational bottleneck is exacerbated when multitudes of RT calculations are required, such as in general circulation models (GCMs) and retrieval methods. In an effort to increase the speed of RT calculations without sacrificing too much accuracy, we investigate the applicability of approximate longwave scattering methods developed for the Earth science community to hot Jupiter atmospheres. We test the absorption approximation and variational iteration method (VIM) applied to typical cloudy hot Jupiter scenarios, using 64-stream DISORT calculations as reference solutions. We find the four-stream VIM variant is a highly promising method to explore for use in hot Jupiter GCM and retrieval modeling, and it shows excellent speed characteristics, with typical errors ∼1% for outgoing fluxes and within ∼50%, but with larger errors in the test case of a deep cloud layer, for vertical heating rates. Other methods explored in this study were found to typically produce similar error characteristics in vertical heating rates.


INTRODUCTION
The calculation of thermal internal atmospheric fluxes, traditionally called 'longwave' radiation, is a key part of investigating the emission properties of exoplanet atmospheres.This type of calculation is critical to modelling and fitting secondary eclipse (dayside) spectral data of tidally locked exoplanets (e.g.Crouzet et al. 2014;Kreidberg et al. 2014;Coulombe et al. 2023) as well as full orbit thermal spectral phase curves of exoplanets (e.g.Stevenson et al. 2014;Arcangeli et al. 2019;Mikal-Evans et al. 2022).
With the advent of JWST, exoplaneteers now have access to an expansive wavelength range, unparalleled sensitivity and continuous staring time for these types of observations than previous keystone observatories.However, this comes with an increase in the challenge of appropriately modelling this data for both forward and retrieval efforts.The wavelength dependent opacity and scattering properties of potential clouds and hazes in these atmospheres can have a large impact on an objects spectra and interpretation of observational data.
Several recent studies have shown the importance of incorporating longwave scattering when modelling cloud properties.Kitzmann et al. (2013) compare a 24-stream DISORT calculation with a two-stream method to investigate the longwave greenhouse effect of CO 2 clouds on exoplanet atmospheres.They find large deviations between the DISORT and two-stream methods alters the strength of the longwave greenhouse effect, impacting the climate of potential worlds containing CO 2 clouds.Taylor et al. (2021) provide analytical expectations for the impact of a scattering component on the outgoing flux of hot Jupiter atmospheres.They show that degeneracies and biases in retrieval results can occur when the longwave scattering is not taken into account and propose a method to retrieve optical constants of the scattering component.Taylor & Parmentier (2023) explore the possible presence of infrared scattering clouds on the dayside of HD 209458b and WASP-43b through retrieval modelling, finding that HD 209458b shows no evidence of a thermal scattering component, while the evidence for WASP-43b is dependent on the specific data reduction used.
Typically, 1D column approaches such as radiativeconvective-equilibrium (RCE) modelling are numerically expedient enough to include more sophisticated RT frameworks (e.g.Drummond et al. 2016;Mollière et al. 2017;Gandhi et al. 2019;Malik et al. 2019).However, 1D retrieval modelling of exoplanet emission spectra data requires many (sometimes millions) of RT iterations and computational efficiency is an important consideration when developing retrieval models (e.g. for emision spectra; Line et al. 2014;Waldmann et al. 2015;Gandhi & Madhusudhan 2018;Kitzmann et al. 2020;Cubillos et al. 2022).With JWST data, the increased spectral resolution as well as spectral wavelength range has substantially increased the RT model burden on retrieval efforts.With the observation of data rich phase curves of exoplanets, consistent 2D retrievals from thermal phase curve data has also been approached (e.g.Irwin et al. 2020;Feng et al. 2020;Changeat et al. 2021;Dobbs-Dixon & Blecic 2022), further increasing the number of RT calculations required to be performed.
3D General Circulation Models (GCMs) are also a class of models that require fast RT schemes.Typically, a 1D RT model is applied for every vertical column of the GCM, which for usual hot Jupiter GCM grid resolutions is thousands of columns (e.g.6144 for the C32 cubedsphere resolution commonly used by SPARC/MITgcm (Showman et al. 2009)) requiring calculation.For an integration time of 100 simulated Earth days, the total number of RT calculations can be in the 10s of million range depending on the radiative timestep used.
In this brief study, we investigate the use of longwave scattering approximation methods in use in the Earth atmospheric sciences and apply them to typical hot Jupiter conditions.In Section 2, we briefly review the RT equation for longwave radiation propagation.Section 3 provides the approximate solution of RT equation, the absorption approximation (AA) and extended absorption approximation (EAA).Section 4 shows the variational iteration method (VIM).In Section 5, we test each method for a cloudy HD 189733b scenario, comparing the outgoing longwave radiation (OLR) flux and vertical heating rates to 64 stream DISORT reference calculations.Section 6 contains the discussion and Section 7 the conclusions of our study.

RADIATIVE-TRANSFER EQUATION
At the core of the longwave radiative-transfer equation in plane-parallel geometry is an expression of energy conservation through a medium by thermal and scattering processes.The radiative-transfer equation for the thermal component in a plane-parallel atmosphere is expressed as (e.g.Li 2002) where µ = cos θ is the cosine angle of propagation of the beam, τ 0 the optical depth (used as the vertical coordinate system), I [erg s −1 cm −2 sr −1 ] the intensity of the radiation, ω 0 the single scattering albedo, B(T) [erg s −1 cm −2 sr −1 ] the Planck function at temperature T [K] and P(µ, µ ′ ) the scattering phase function, which gives the probability of radiation being scattered from the µ ′ direction towards the µ beam direction.
The term with the Planck function contains the contribution of medium's thermal energy to the intensity of the beam.This is assumed to be an isotropic source, so does not contain any directional dependent component.The second term contains the contribution of scattered light from other directions.This is governed by the properties of the scattering phase function P(µ, µ ′ ), which can take many forms, but usually assumed to either be isotropic or equal to the Henyey-Greenstein analytic phase function (Henyey & Greenstein 1941).This is integrated across all µ ′ directions, providing the net contribution of scattered light towards a given µ beam direction.In the following sections, we attempt to keep the expressions generally consistent with the Heng et al. (2014) and associated papers notation, adding the concepts for each approximation following Li & Fu (2000) and Zhang et al. (2017).

ABSORPTION APPROXIMATION
In the absorption approximation (AA), opacity sources are assumed to only contribute their absorption opacity to the atmospheric transmission, neglecting their scattering component to the total extinction opacity.In this approach, the scattering source function is ignored and incorporated into an alteration to the extinction opacity of the atmosphere, this reduces the RT equation to (e.g.Zhang et al. 2017) where ϵ = 1 -ω 0 is co-albedo.
To take into account non-isothermal layers an exponential functional form for B(τ 0 ) is used (e.g.Fu & Liou 1993;Zhang et al. 2017) where β = -[ln(B 2 /B 1 )]/τ , τ = τ 2 -τ 1 is the optical depth of the layer.Following Heng et al. (2014), here the index 1 denotes the upper level (located at lower optical depth) and 2 the lower level (located at higher optical depth).A commonly used functional form is the so-called linear in τ approximation given by (e.g.Toon et al. 1989) Li ( 2002) discuss the differences between the exponential and linear forms.We define the regular transmission function, T , of an atmospheric layer as and the absorption only transmission function, T a , of an atmospheric layer as The solution of Eq. 2, assuming the exponential Planck function form (Eq. 3), is (e.g.Zhang et al. 2017) for the upward and downward directions respectively.The upper boundary condition at layer index 0 is typically unless some other thermal source is required to be taken into account above the computational boundary.For a gas giant atmosphere with given internal temperature, T int [K], the lower boundary condition is given by where n is the index of the lowest boundary level.This condition ensures energy conservation across the boundary and the net flux is equal to the internal thermal energy.Integration proceeds by first calculating the downward intensity from the upper boundary to the lower, then using the downward intensity at the lowest level through Eq. 10 to the perform the return upward intensity sweep.
The AA method is simple to implement as there is no coupling between the upward and downward direction from a scattering component.Overall, AA effectively replaces the extinction opacity used in the regular transmission function with the absorption opacity, meaning any method that doesn't compute a scattering source function can use AA in a simple manner through a change in the transmission function.
The upward and downward flux, F ↑,↓ in each layer is calculated using Gaussian quadrature through respectively.For the two-stream variant (2AA) a single Gaussian quadrature value and weight is used, traditionally µ 1 = 1/1.66,w 1 = 1 following Elsasser (1942).
The four-stream variant (4AA) uses a double Gaussian quadrature scheme with µ 1 = 0.21132487, w 1 = 0.5 and µ 2 = 0.78867513, w 2 = 0.5, which comes at the additional cost of performing four sweeps in the vertical column for the flux calculation.Lastly, the heating rate ∂T/∂t [K s −1 ] for a layer is given by where g [cm s −2 ] is the gravitational acceleration, c p [erg is the net radiative flux into the layer and ∂p = p 2 -p 1 is the pressure difference across the layer.The AA method is only valid for weak scattering and does not consider the directional component of the scattering, assuming all scattering is in the exact forward direction (i.e.g 0 = 1).δ-M+ scaling of the albedo and optical depth for AA can be applied following Wiscombe (1977) and Lin et al. (2018).We denote this scaled version as δ-2AA/δ-4AA (Zhang et al. 2017).

Extended Absorption Approximation
To somewhat take into account the directional component in the AA formulation, we explore a slight extension to the co-albedo, ϵ ′ , following the hemispheric twostream solution (e.g.Pierrehumbert 2010;Heng et al. 2014;Heng & Kitzmann 2017) where g 0 ≥ 0, is assumed for a forward or isotropic scattering component only.The AA method is returned when g 0 = 1, with the regular definition of the co-albedo ϵ ′ = ϵ.When g 0 = 0 for isotropic scattering we more clearly see the effect of the altered co-albedo, this results in ϵ ′ = √ ϵ, increasing the co-albedo, in effect reducing the single scattering albedo.This increase attempts to account for the non-preference in the scattering direction through reducing the transmission function.Overall, this altered co-albedo attempts to relax the g 0 = 1 assumption in the AA approach in a simple manner, we name this alteration the 'Extended Absorption Approximation' (EAA) or when δ-M+ scaled δ-nEAA for n streams.In Appendix A we explore specifically the differences between AA and EAA for our proposed tests.Overall, we find EAA to be a general improvement over AA.
4. VARIATIONAL ITERATION METHOD Zhang et al. (2017) propose a variational iteration method (VIM) (e.g.He 1999) solution to the longwave RT equation with approximate scattering.In VIM, a differential equation is split into a linear and non-linear term.An iterative scheme can then be used to improve the solution (the so called correction functional) by successive orders through applying a Lagrange multiplier to the previous solution.Importantly, the VIM solution does not consider true multiple scattering across layers but only the scattering properties within a layer and their onward affect on adjacent layers.
Here we apply the four-stream variant of VIM outlined in Zhang et al. (2017), however, we note the theory generalises to n-streams.For brevity, we do not repeat the extensive derivation performed in Zhang et al. (2017), only stating the main result of the derivation.The fourstream longwave radiative-transfer equation is approximated as (Zhang et al. 2017) Here, i and j denote the positive cosine angles (i.e.µ i = µ i ) and -i and -j the negative cosine angles (i.e.µ −i = -µ i ), with the coefficients b, c, d and ψ given by The scattering coefficients S ↑,j , S ↓,j are given by and where In practice, VIM uses the (E)AA method to provide the zeroth order estimate of the intensity (here denoted by I 0 ↑,↓ ) in the upward and downward directions.This zeroth order solution is then used to estimate the scattering coefficients S ↑,↓ , with the final RT solution performed including the estimated scattering source function.VIM therefore attempts to calculate a first order scattering 'correction' to the zeroth order (E)AA method, following the iterative methodology of VIM.In addition, in comparison to the Toon et al. (1989) two-stream method, the Henyey-Greenstein phase function is used for the scattering component rather than a hemispheric-mean approximation.The main benefit of the VIM method is that when no scattering component is present in a layer, the scattering source function need not be calculated, substantially improving computational efficiency.Again, the δ-M+ scaled versions of VIM are denoted δ-nVIM.

TESTING THE APPROXIMATIONS
In this section, we briefly describe the other RT methods used in this study to benchmark against the AA and VIM approximations.We then perform tests on the computational speed, OLR and vertical heating rates for each method.Table 1.Relative speed scaled to δ-2EAA of each of the longwave algorithms for the grey atmosphere speed test.

Other RT methods
The Toon et al. (1989) (hereafter Toon89) method applies a two-stream approximation to calculate the scattering source function assuming a hemispheric mean scattering phase function.Then, the scattering source function is added to thermal component in the source function technique to find the intensity.Any number of streams can be used for the source function technique part of the calculation, but the scattering component is always assumed as the hemispheric mean two-stream solution.Here, we use four-stream version for the source function technique part of the calculation for a fairer comparison to the other four-stream approximate methods used in this study.
We note Marley et al. (2021), which also used Toon89, have applied 12-streams for the source functions part of the method in their brown dwarf RCE modelling.Toon89 has also been applied to numerous exoplanet GCM models (e.g.Showman et al. 2009;Roman et al. 2021;Wolf et al. 2022; Lee 2023) and retrieval/RCE models (e.g.Line et al. 2013;Burningham et al. 2017;Batalha et al. 2019) with various number of source function streams.However, the two-stream derived hemispheric mean scattering source function remains the same for all these studies.
A specific version of the discrete ordinance method code DISORT (Stamnes et al. 1988;Stamnes et al. 2000) optimised for two-streams (DISORT-twostr) is publicly available1 (Kylling et al. 1995).Some brown dwarf and exoplanet GCM studies have used this code as the RT driver (e.g.Tan & Showman 2021), mostly with application coupling to the MITgcm (Komacek et al. 2022).We use the DISORT code with 64 streams as our reference solutions to compare methods.

Grey speed tests
For a fair test of the pure RT component of each method, we time the longwave algorithm for each method for a grey opacity atmospheric model.This is ideal to test the direct speed properties as it greatly simplifies the RT problem and only one pass of the longwave subroutine is required per iteration.We use the two-stream AA as the baseline, adding the four-stream VIM, four-stream Toon89 method and the two-stream DISORT version to the speed test.We use a Guillot (2010) analytical radiative-equilibrium grey opacity profile, with first test without and scattering component and the second test with a constant single scattering albedo of 0.5 and asymmetry factor of 0.5.Table 1 shows the results of this speed test.The main result of this test is that the speed of the δ-4VIM method is on par with the Toon89 method when every layer contains a scattering component, but speed is much improved in δ-4VIM method the when no scattering components are present.

Outgoing longwave radiation fluxes
In this section, we test the accuracy of each scheme in producing OLR fluxes under cloud free and cloudy scenarios.We follow a similar methodology to Amundsen et al. (2014), where a representative HD 209458b hot Jupiter T-p profile from the Heng et al. (2011) polynomial fit to the Iro et al. (2005) simulations was used  2015) analytical picket-fence scheme as our test temperature-pressure (T-p) profile.Figure 1 shows this T-p profile.A deep boundary pressure of 1000 bar is chosen provide a large optical depth region to avoid undue boundary condition influences on the upper atmosphere.We use 54 layers, emulating a typical GCM vertical resolution.We use the same hybrid-sigma grid pressure spacing as used in the hot Jupiter Exo-FMS simulations (e.g. Lee et al. 2021).Strong evidence of cloud coverage on HD 189733b is seen in contemporary observational data (e.g.Sing et al. 2016), making the HD 189733b hot Jupiter parameter regime ideal for testing the longwave schemes.
For the gas phase opacity, we use a 34 band scheme the same as Amundsen et al. (2014) with pre-mixed correlated-k tables (Amundsen et al. 2017), assuming mixing ratios at chemical equilibrium.These tables use a 4+4 k-coefficient scheme the same as Kataria et al. (2013) and Marley et al. (2021) and have been previously used in hot Jupiter simulations using Exo-FMS (e.g. Lee et al. 2021).Pre-mixed tables such as these are in common use in GCMs of exoplanet atmospheres (e.g.Showman et al. 2009;Schneider et al. 2022) For the cloud opacities, we apply Mie theory using a Fortran conversion of LX-MIE code from Kitzmann & Heng (2018).We convert a given cloud mass mixing ratio, q c [kg/kg], to a total cloud particle number density, assuming a mean cloud particle size and standard deviation, following a log-normal distribution (e.g.Komacek et al. 2022).This results in realistic cloud opacities, sin-gle scattering albedo and asymmetry parameters across all the bands.
We investigate three cloud scenarios, a deep cloud case, an upper extended cloud case and a combination deep and extended case.For the deep cloud case, we parameterise a cloud base at 10 bar and cloud top at 0.01 bar, using a mean particle size of 2µm and standard deviation of σ = 1.The q c is given a maximum value at 10 bar of 10 −4 kg/kg, with a power law drop off with gradient = 0.5 to the cloud top pressure level.This mimics a compact deep cloud layer scenario, such as those seen in Ackerman & Marley (2001) models of sub-stellar objects (e.g.Gao et al. 2018), when mixing is slow and sedimentation efficient.
For the extended cloud case, the cloud base is at 0.01 bar and cloud top at 10 −6 bar, using a mean particle size of 0.1µm and standard deviation of σ = 0.1.We use a constant q c of 10 −4 kg/kg.This intends to represent a well mixed and lofted upper atmosphere extended cloud, commonly seen for models that contain strong mixing or those seen in cloud studies using GCM model output (e.g.Gao et al. 2018;Powell et al. 2018;Lee 2023).
For both cloud layers, we assume a MgSiO 3 cloud composition, taking the optical constants of amorphous MgSiO 3 from the Kitzmann & Heng (2018) collection.Figure 1 shows the q c structure of the two cloud parameterisations.After the calculation of the cloud opacity, several wavelength bands exhibit high single scattering albedos coupled with high asymmetry factors (e.g.ω 0 > 0.8, g 0 > 0.8) inside both cloud regions, which is an ideal case for testing the scattering properties of each method.For illustration, Figure 2 shows several wave-length bands and their cloud opacity, single scattering albedo and asymmetry parameter with pressure.
In Table 2, we present the total OLR error for each method compared to the reference 64-stream DISORT calculations.Here we see the EAA method varies strongly depending on the cloud component, with the deep clouds more accurately represented than the extended cloud, or even cloud free profile.VIM performs similarly across all tests, with an error of around ∼1%. Toon89 fares the best, with OLR errors of only ∼0.3% compared to the reference DISORT calculations, consistent across all tests.The DISORT-twostr methods produces the largest errors at around ∼4%.

Heating rates
In this section, we test the accuracy of each scheme in producing vertical heating fluxes under cloud free and cloudy scenarios.We use the same T-p, opacity and cloud component set-up as the previous section.Figure 3 presents the results of the tests, showing both the absolute value of the heating rates produced by each method and the error relative to the 64-stream DISORT reference solution.For the cloud free case, all methods produce very reliable results, with DISORT-twostr producing the largest errors of about 50% in deeper regions.
For the deep cloud case we start to see where the main differences between the methods arise, where all methods produce significant errors near the base of cloud layer where the opacity of the cloud is strongest.Interestingly, the Toon89 and DISORT-twostr methods produce heating locations similar to DISORT while the EAA and VIM methods results in a slightly lower pressure heating spike due to the cloud layer.The large errors in the AA and VIM method at these pressures reflect this difference in heating level.This suggests that multiple-scattering may be playing an important role here, allowing a heating component below the cloud base not captured by AA or VIM.However, the multiplescattering methods still show significant error at the cloud base, with up to 50% errors compared to the DIS-ORT reference solution.Above the deep cloud base, errors are highly comparable between each method, suggesting this effect does not persist through the cloud layer.
For the extended cloud case, the Toon89 and DISORT-twostr compare much better than EAA and VIM to the reference solution, with EAA and VIM producing a greater heating component at the cloud base than the Toon89 and DISORT-twostr methods.When comparing the results of both cloud layers, both errors in each scenario are generally added together, showing similar characteristics for each test separately.

DISCUSSION
Our OLR tests showed that the most consistent methods were Toon89 at ∼0.3% error and VIM at ∼1% error across all tests.However, the computational time to accuracy trade-off must be taken into account for GCM and retrieval modelling efforts.From the OLR tests, VIM is promising as in the best case scenario VIM is approximately 3x faster than Toon89 while at worst (a profile with all layers containing a scattering component) is on-par with the Toon89 computational time.We therefore suggest that EAA and VIM show strong promise for improving the speed of emission spectra retrieval modelling where retrieving vertical cloud properties is required.This is possibly an improvement for 2D retrieval efforts where hemispheres of cloud free and cloudy profiles are mixed.EAA and VIM may be most useful for parameterised cloud profiles or grey clouds that have a specified single scattering albedo that is not too large (we suggest ω 0 ≲ 0.8).Usefully, both EAA and VIM can be used to retrieve an asymmetry parameter for a parameterised cloud layer.
Our vertical heating rate tests shows the strongest promise for approximate methods, as the errors are very comparable to the Toon89 and DISORT-twostr methods.This suggests EAA and VIM to be very suitable for GCM modelling where RT efficiency is paramount.The computational efficiency of VIM will also benefit from the vertical and hemispheric inhomogeneous nature of clouds on hot Jupiter.
In Rooney et al. (2024), the Toon89 method was compared to a spherical harmonic (SH) method with source function technique.They found that the Toon89 method produced spectral OLR differences of ∼20% or more at certain wavelengths compared to their DISORT reference calculations for their cloudy brown dwarf examples.However, here we report a much closer agreement in OLR < 1% between the Toon89 and DISORT methods.We suggest a few reasons for this discrepancy • The larger gravity used in the Rooney et al. (2024) brown dwarf models, making the OLR more sensitive to strong scattering components deeper in the atmosphere.
• The consistently large single scattering albedo and asymmetry factors produced by their cloud model across a wide infrared wavelength regime.
We suspect these two factors, as well as the assumption of a hemispheric mean phase function in the Toon89 method compared to the HG function in the SH method, are responsible for the much larger errors seen in Rooney et al. (2024).Owing to this, the EAA and VIM methods will need to be separately tested for brown dwarf conditions similar to Rooney et al. (2024), though we speculate that since a lot of the errors in the Toon89 method stem from the use of the hemispheric-mean phase approximation (Rooney et al. 2024) that VIM may be more accurate than Toon89 in some cases, but probably not as accurate as the SH method itself.
The VIM method is also known to be less accurate for smaller optical depth layers with large single scattering albedos (Zhang et al. 2017), such as high altitude Earth clouds, this is primarily due to the zeroth order solution using AA which also performs worse at these conditions (Zhang et al. 2017).This can be seen in Figure 3, where the EAA and VIM method produce a larger error for the low pressure extended cloud case compared the Toon89 and DISORT-twostr methods.We therefore suggest this be taken into account when considering EAA or VIM for heating rate purposes, for example, this error may be acceptable for GCM modelling but not adequate for RCE modelling, where an inaccurate equilibrium profile would more greatly affect the interpretation of results.
Due to its application in Earth atmospheric modelling, the VIM method may be applied to Mars CO 2 clouds and their affect on the Mars climate (e.g.Forget & Pierrehumbert 1997;Kitzmann 2016;Heng & Kitzmann 2017).However, as an approximate method, VIM is likely to not be as accurate as the DISORT simulations performed in Kitzmann (2016).As stated in Zhang et al. (2017), low optical depths with large albedo is a parameter regime where VIM can be in significant error, which is the case for the Mars CO 2 cloud problem.

SUMMARY AND CONCLUSIONS
In this brief study, we investigated two approximations for calculating internal atmospheric (longwave) radiation with a scattering component in a hot Jupiter context, the (extended) absorption approximation (EAA) and the variational iteration method (VIM).
For retrieval modelling, EAA and VIM may be most useful for parameterised cloud profiles or grey clouds that have a specified single scattering albedo that is not too large.In the JWST era, where the number of wavelength points required for calculation is steadily increasing, and as retrieving specific cloud properties becomes more important, convergence runtimes of retrieval models are also on the rise.EAA and VIM offer a simple way to speed up the forward modelling process without too much sacrifice in accuracy.
For gas giant exoplanet GCMs, VIM offers a quick and simple method to take into account cloud particle scattering in the longwave, with comparable heating rate errors to other well used models in the field.We therefore suggest VIM is a highly promising method for use in exoplanet GCMs.
Due to the generalised nature of the RT methods used here, our findings are likely applicable to atmospheres across the exoplanet spectrum, from ultra hot Jupiters to warm Neptunes and rocky planets.However, we caution high gravity, strong scattering cloudy brown dwarf conditions such as those in Rooney et al. (2024) will need to be tested further before AAE and VIM can be recommended for brown dwarf modelling efforts.E.K.H. Lee is supported by the SNSF Ambizione Fellowship grant (#193448).
Software: DISORT and DISORT-twostr: http: //www.rtatmocn.com/disort/.Exo-FMS column emulation codes: https://github.com/ELeeAstro.Fortran 90 version of LX-MIE: https://github.com/ELeeAstro/gCMCRT/blob/main/src optools V2/lxmie mod.f90In Figure 4 we show the errors between the AA and EAA methods in vertical heating rates compared to the 64 stream DISORT reference values.These results generally show an improvement of EAA over AA, in particular for the deep cloud scenario.This also shows the superiority of the four-stream variants over the two-stream, the foursteam variants show a much increased accuracy over the two-steam, in particular at the low pressure regions.We can therefore recommend using the four-stream EAA variant for modelling calculations and as the zeroth order estimate for the VIM method.

Figure 1 .
Figure 1.HD 189733b global averaged T-p profile (orange solid line) produced from the picket-fence scheme and condensed mass fraction qc (blue dashed line) profiles for the parameterised deep and extended cloud layers used in the RT tests.

Figure 2 .
Figure 2. Left: cloud extinction opacity, κ ext,cld , with pressure for the combined upper and lower cloud layer case for several wavelength bands.Right: single scattering albedo, ω0 (solid lines), and asymmetry parameter, g0 (dashed lines), with pressure for several wavelength bands for the combined cloud case. to test the radiative-transfer schemes.Here, we use a global averaged HD 189733b profile derived from the Parmentier et al. (2015) analytical picket-fence scheme as our test temperature-pressure (T-p) profile.Figure 1 shows this T-p profile.A deep boundary pressure of 1000 bar is chosen provide a large optical depth region to avoid undue boundary condition influences on the upper atmosphere.We use 54 layers, emulating a typical GCM vertical resolution.We use the same hybrid-sigma grid pressure spacing as used in the hot Jupiter Exo-FMS simulations (e.g. Lee et al. 2021).Strong evidence of cloud coverage on HD 189733b is seen in contemporary observational data (e.g.Sing et al. 2016), making the HD 189733b hot Jupiter parameter regime ideal for testing the longwave schemes.For the gas phase opacity, we use a 34 band scheme the same asAmundsen et al. (2014) with pre-mixed correlated-k tables(Amundsen et al. 2017), assuming mixing ratios at chemical equilibrium.These tables use a 4+4 k-coefficient scheme the same asKataria et al. (2013) andMarley et al. (2021) and have been previously used in hot Jupiter simulations using Exo-FMS (e.g.Lee et al. 2021).Pre-mixed tables such as these are in common use in GCMs of exoplanet atmospheres (e.g.Showman et al. 2009;Schneider et al. 2022) For the cloud opacities, we apply Mie theory using a Fortran conversion of LX-MIE code fromKitzmann & Heng (2018).We convert a given cloud mass mixing ratio, q c [kg/kg], to a total cloud particle number density, assuming a mean cloud particle size and standard deviation, following a log-normal distribution (e.g.Komacek et al. 2022).This results in realistic cloud opacities, sin-

Figure 3 .
Figure 3. Absolute value of the vertical heating rates (left) and percent error in the heating rates (right) for each of the test cases.Each method is shown as a coloured solid line.

Figure 4 .
Figure 4. Percent error in the heating rates for the 2 and 4 stream AA and EAA methods.Each method is shown as a coloured solid line.

Table 2 .
Relative percentage error on the total OLR for each method and scenario compared to the reference 64 stream DISORT results.