Clouds form on the hot Saturn JWST ERO target WASP-96b

WASP-96b is a hot Saturn exoplanet, with an equilibrium temperature well within the regime of thermodynamically expected extensive cloud formation. Prior observations with Hubble/WFC3, Spitzer/IRAC, and VLT/FORS2 have been combined into a single spectra for which retrievals suggest a cold but cloud-free atmosphere. Recently, the planet was observed with the James Webb Space Telescope (JWST) as part of the Early Release Observations (ERO). 1D profiles are extracted from the 3D GCM expeRT/MITgcm results and used as input for a kinetic, non-equilibrium model to study the formation of mineral cloud particles of mixed composition. The ARCiS retrieval framework is applied to the pre-JWST WASP-96b transit spectra to investigate the apparent contradiction between cloudy models and assumed cloud-free transit spectra. Clouds are predicted to be ubiquitous throughout the atmosphere of WASP-96b. Silicate materials contribute between 40% and 90%, hence, also metal oxides contribute with up to 40% in the low-pressure regimes that effect the spectra. We explore how to match these cloudy models with currently available atmospheric transit spectra. A reduced vertical mixing acts to settle clouds to deeper in the atmosphere, and an increased cloud particles porosity reduces the opacity of clouds in the near-IR and optical region. Both these effects allow for clearer molecular features to be observed, while still allowing clouds to be in the atmosphere. The atmosphere of WASP-96b is unlikely to be cloud free. Also retrievals of HST, Spitzer and VLT spectra show that multiple cloudy solutions reproduce the data. JWST observations will be affected by clouds, where within even the NIRISS wavelength range the cloud top pressure varies by an order of magnitude. The long wavelength end of NIRSpec and short end of MIRI may probe atmospheric asymmetries between the limbs of the terminator on WASP-96b.


Introduction
The hot Saturn WASP-96b was one of the first1 extrasolar planets for which the James Webb Space Telescope (JWST) observations were released to the community on the 12 th July 2022.The planet's primary transit was observed using the telescope's NIRISS/SOSS instrument as part of the JWST Early Release Observations (ERO) (Pontoppidan et al. 2022).It therefore presents itself as a desirable target to analyse from a modelling perspective, in preparation for newly reduced transit spectra expected to become available in the near future.WASP-96b has a mass of 0.48 ± 0.03 M Jup , a mean radius of 1.2 ± 0.06 R Jup , and orbits a G-type stars in 3.4 days at a distance of 1120 lyr from Earth (Hellier et al. 2014).The orbit appears non-eccentric and its global temperature is T glob ≈ 1300K (1285 ± 40, Hellier et al. 2014).Optical transmission spectra of WASP-96b have been observed pre-JWST by the FORS2 spectrograph on the ground-based Very Large Telescope (VLT) (Nikolov et al. 2018).The resolution of these VLT spectra allowed the Na doublet to be studied in detail; Na abundances were derived and a cloud-free atmosphere Fig. 1: Pressure-temperature structure of the 120 1D profiles extracted from the GCM results.Orange lines represent the dayside profiles, blue lines the nightside profiles, and the grey lines the terminator profiles.Black dashed and dot-dashed lines highlight the anti-stellar points respectively.Day-and night-side hemisphere averages shown in red and blue dashed lines.Inset shows the upper atmosphere GCM of the profiles.recent observations by Nikolov et al. (2022) suggest that WASP-96 has enhanced metallicity.
A recently analysed observation of brown dwarf VHS 1256-1257 b by Miles et al. (2022) using JWST's NIRSpec IFU and MIRI MRS modes, revealed solid-states silicate features in the longer wavelength part of the atmospheric spectrum, giving strong evidence for small silicate particles in the atmosphere.
We demonstrate in this paper that mineral cloud particles that are made of a mix of metal-oxide and silicate components are expected to form in the atmosphere of the hot Saturn WASP-96b.We present a cloud map for this planet which demonstrate its complete coverage with clouds.We identify WASP-96b as sitting right at the boundary between the weather/climate regime for cool planets and the less homogeneous transition planets that were introduced in Helling et al. (2022)..The absence of a strong morning/evening terminator asymmetry should improve the chances of observing spectral features of the ensemble of cloud particles that our microphysical models predict should populate the terminator regions.

Approach
We examine the cloud structure on the hot Saturn WASP-96b by adopting a hierarchical approach similar to works on the hot Jupiters HD 189733b and HD 209458b (Lee et al. 2015;Helling et al. 2016), and the ultra-hot Jupiters WASP-18b (Helling et al. 2019) and HAT-P-7b (Helling et al. 2019;Molaverdikhani et al. 2020): The first modelling step produces a cloud-free 3D GCM result for a model representing WASP-96b.These results are used as input for the second modelling step which is a kinetic cloud formation model consistently combined with equilibrium gas-chemistry calculations.We utilise 120 1D (T gas (z), p gas (z), v z (z))-profiles for WASP-96b similar to our previous works.T gas (z) is the local gas temperature [K], p gas (z) is the local gas pressure [bar], and v z (z) is the local vertical velocity component [cm s −1 ].
Such a hierarchical approach has the limitation of not explicitly taking into account the potential effect of horizontal winds on cloud formation.However, processes governing the formation of mineral clouds are determined by local thermodynamic properties which are the result of 3D dynamic atmosphere simulations.Cloud particle properties such as particle size or particle composition could be somewhat smeared out in longitude compared to the results shown here.The temperature structure may change if the cloud particle opacity is fully taken into account in the solution of the radiative transfer.This may change the precise location of the cloud in pressure space but not the principle result of clouds forming in WASP-96b.
3D atmosphere modelling: We utilise expeRT/MITgcm (Carone et al. 2020;Baeyens et al. 2021) to model WASP-96b.expeRT/MITgcm builds on the dynamical core of MITgcm and has been adapted to model tidally locked gas giants.Recent extensions in Schneider et al. (2022) include non-grey radiative transfer coupling.The model parameters used for the GCM, representative of the hot Saturn WASP-96b, are: R p = 8.58×10 9 cm, P rot = 3.4 days, log 10 (g [cm s −2 ]) = 2.9, and the substellar point irradiation temperature T irr = 1819 K (Eq.20 Schneider et al. 2022).The model is run for 1400 days, and we note that the deepest layers may not be fully converged.A constant mean molecular weight of µ = 2.3 is assumed for the 3D atmospheric modelling.This constant value is a reasonable assumption given that the thermodynamic structure of the atmosphere does not cause the gas composition to deviate from a H 2 -dominated gas.Additional details for the 3D GCM setup can be found in Table A.1.Nikolov et al. (2022) fit an isothermal black-body to the dayside of WASP-96b using Spitzer observations, finding a dayside brightness temperature of 1545 ± 90 K.The equilibrium temperature of the planet can be calculated using T eq = R /a p T ( f (1 − A B )) 1/4 , the final factor includes the heat redistribution factor ( f ) and the bond Albedo of the planet (A B ). Assuming full heat redistribution ( f = 1/4) and full absorption of incident radiation (A B = 0), one arrives at the equilibrium temperature of the planet: T eq = 1286 K.As input to the initialisation of the GCM, the substellar irradiation temperature is used, this is the temperature of a black-body re-radiating the energy received from the star, by the planet at the substellar point (T irr = 1819 K, see Table A .1).Comparing these numbers to the observed brightness temperature is difficult, the brightness temperature sits between the planet equilibrium temperature and the substellar irradiation temperature.This indicates some level of heat re-distribution around the planet, as well as the potential impact of cloud albedo.However, observations of similar temperature exoplanets have revealed low dayside albedos (e.g.Schwartz & Cowan (2015); Fraine et al. (2021); Brandeker et al. (2022)).

Kinetic cloud formation:
We apply the same set-up of our kinetic cloud formation model (nucleation, growth, evaporation, gravitational settling, element consumption and replenishment) and equilibrium gas-phase calculations as for our grid study in Helling et al. (2022).It is essential that the seed forming species are also considered as surface growth material, since both processes (nucleation and growth) compete for the participating elements (Ti, Si, O, Na, K, and Cl in this work).We consider the formation of 16 bulk materials (s=TiO NaCl [s]) that form from 11 elements (Mg, Si, Ti, O, Fe, Al, Ca, S, Na, K and Cl) by 132 surface reactions.In total, we are solving 31 ODEs to describe the formation of cloud condensation nuclei (J * (z) = i J i , i=TiO 2 , SiO, KCl, NaCl) that grow to macroscopic sized cloud particles made of a mix of materials which changes depending on the local atmospheric gas temperature and gas pressure.Fig. 3: Volume ratio of the cloud forming materials for WASP-96b in 2D terminator slices, for the same models as in Fig. 2, materials are grouped as in (Helling et al. 2021) Retrieval and atmospheric spectra modelling with ARCiS: We follow the method of Min et al. (2020) by running a constrained retrieval on the previously published VLT/HST/Spitzer transmission spectra of WASP-96b (collated by Nikolov et al. (2022)) using Bayesian retrieval code ARCiS (ARtful modelling Code for exoplanet Science).The constrained retrieval includes the cloud formation models of Ormel & Min (2019), with the variety of cloud species included in this work the same as in Min et al. (2020) with additional Na and K silicate clouds .It is a simplified approach in comparison to the previous hierarchical modelling stages of this work, but it enables us to compare some retrieved parameters to our kinetic cloud formation models, such as cloud particle size, and expected cloud optical thickness.
It also allows us to explore the effect on the spectra of some of these parameters, which we do using the forward modelling mode of ARCiS.We generate simulated transmission spectra with varied atmospheric parameters, which can be compared to observations.Here, we focus on the effects of clouds in the atmosphere and how they can be present but still let us observe strong molecular features.The distribution of hollow spheres (DHS) method is employed in ARCiS for computing the cloud opacities (Min et al. 2005).This allows for the inclusion of nonspherical cloud particles, and the ability to vary the porosity of the particles.A similar simple model for particle porosity incorporating also the impacts of increased surface area and reduced settling on the clouds as well as optical effects was recently explored in Samra et al. (2020).We also explore the impact of this microphysically consistent model on the cloud optical depth in Section 3.4.The chemistry for all retrievals is constrained, with elemental abundances computed after the cloud formation process in order to take into account elements which may have formed into clouds.The equilibrium chemistry code GGchem (Woitke et al. 2018) is used to compute the molecular abundances given these elemental abundances.A summary of the retrieved parameters and their priors is given in Table 1.For comparison, we also run a cloud-free retrieval, with all other parameters the same.We present the results of these retrievals in Section 4.1.

Cloud properties on the hot Saturn WASP-96b
3.1.The 3D GCM atmosphere structure The 3D GCM results for the hot Saturn WASP-96b suggest a maximum day/night temperature contrast of ≈ 400K occurring in the low-pressure regimes at p gas < 10 −2 bar. Figure 1 demonstrates this by the 120 1D profiles extracted from the 3D GCM solution.The morning terminator (grey dashed line) shows some moderate temperature inversion of ≈ 100 K resulting from the cold gas being transported from the nightside to the dayside.The temperature differences between the morning and the evening terminators are maximum ≈ 200 K.
Based on its global temperature T glob ≈ 1300 K, the hot Saturn WASP-96b falls onto the boundary of the cool planet climate regimes identified in Helling et al. (2022).The cool planet climate regimes are characterised by a globally homogeneous formation of cloud condensation nuclei that leads to a globally homogeneous cloud coverage.WASP-96b does have some day/night temperature difference as shown in Fig. 1.The day/night temperature difference does cause a temperature difference in the evening and morning terminator regions as shown in the 2D terminator slice plot for local gas temperature in Fig. 2 (top left).The evening terminator is affected by the warm gas being advected from the dayside to the nightside, the morning terminator by the cold gas coming from the nightside.This causes a slightly elliptically shaped high-temperature area in the equatorial regions (θ = 0 o ).

The mineral clouds on WASP-96b
In this section, we discuss the results for mineral cloud formation in the hot Saturn WASP-96b to enable a detailed insight into the cloud properties that are likely to shape the JWST spectrum of the planet.

Very mildly asymmetric cloud coverage and mildly non-homogeneous cloud properties
Firstly, we observe that cloud condensation nuclei (CCNs) form for all the 1D profiles used here to probe cloud formation in the 3D atmosphere of WASP-96b (Fig. 2, top right).This is in stark contrast to the ultra-hot Jupiters, like WASP-18b and HAT-P-7b, but in line with other hot Jupiters like HD 189733b and HD 209458 b.Secondly, the rate of seed formation differs between the day and the nightside because of the different thermodynamic conditions.The morning terminator regions have a somewhat larger atmospheric volume where efficient formation of CCNs takes place compared to the evening terminator (i.e. at the morning terminator the region of efficient nucleation extends to deeper pressures than at the evening terminator).The nucleation rate, J * , follows the asymmetric morning/evening terminator temperature differences.
Once the CCNs are formed, they grow to macroscopic cloud particles through condensation, these cloud particles fall into the atmosphere towards higher pressures where growth may accelerate, before the cloud particles evaporate where their materials become thermally unstable at the cloud base (at p gas ≈ 10 2.2 bar, consistently around the terminator).Therefore, the extension of the cloud is determined by the atmospheric volume that is populated by cloud particles (Fig. 2, bottom left), not by where the CCNs form (Fig. 2, top right).Mean cloud particle sizes range from 0.3 µm at p gas = 10 −2 bar, to 15µm at 1bar and 0.1 mm at the cloud base at p gas = 158 bar.
The cloud extension is affected by the vertical mixing efficiency but also by material properties like the cloud particle porosity.The vertical mixing counteracts the element depletion of the gas phase by cloud formation which would inhibit the continuous presence of clouds in an atmosphere (see Sect. 3.3. in Woitke & Helling 2003).Samra et al. (2020) demonstrated that highly porous (low density) cloud would expand the cloud layer upwards, due to the increased surface area provided by porous cloud particles.However, the cloud base is unaffected as it is determined by thermal stability.High-porosity cloud particle would also have a larger mean particles size and would result in a somewhat larger dust-to-gas ratio.
Figure 2 (lower left) shows that the cloud particle mass load (in terms of the dust-to-gas ration ρ d /ρ) of the WASP-96b model atmosphere appears rather symmetrically distributed as a result of the homogeneous thermodynamic structure that is suggested by the 3D GCM for WASP-96b.However, there is still a slight asymmetry between the two terminators in the peak cloud particle mass load, due to the increased particle surface area from the larger nucleation rate on the morning terminator.

General material composition of cloud particles
We do not consider carbon condensate species and have demonstrated the mineral CCN formation is more prominent in oxygenrich atmospheres than the formation of hydrocarbon hazes (Helling et al. 2020).
The 2D terminator slice plots in Fig. 3 visualise the cloud particle material compositions in units of material volume fraction, V s /V tot (Eqs. 25, 26 in Helling &Woitke 2006 andEq. 6 in Helling et al. 2008).The 16 materials that are considered to contribute to the bulk growth of the cloud particles for our WASP-96b model are collected into four groups: metal oxides, silicates, high temperature condensates and salts (for details see caption of Fig. 3).With a view to the possible observability of cloud particles through their spectral fingerprints, important results are: -Cloud particles are made of a mix of all materials that are thermally stable at the local (T gas , p gas ).
CaSiO 3 [s]) dominate almost everywhere the material composition of the cloud particles that form the clouds of WASP-96b.Their dominant contribution, however, varies between 40% in the low-pressure regions (p gas 10 −1.5 bar on the morning terminator, p gas 10 −2.5 bar on the evening terminator) and >90%.We therefore suggest that the cloud particles in the lowpressure regions accessible by transmission spectroscopy should be dominated of a mix of ≈ 40% metal oxides and 40% silicates.These particles have a mean particle size of 10 −2 µm across the terminator region accessible by transmission spectroscopy.

Mineral ratios and the mean molecular weight
Mineral ratios (element ratios) like Si/O, Mg/O, Fe/O and the carbon-to-oxygen ratio, C/O, may help to link observations to evolutionary stages of exoplanets (Khorshid et al. 2021;Öberg & Bergin 2016).The challenge with this inference is that observations have so far only provides snippets of the atmosphere spectrum, i.e. information about limited wavelength intervals.Consequently, inferences of element ratios often hinge on retrieved abundances of individual gas species.This situation, however, may change with the recent JWST observations, although the new JWST observations are still limited in range as they are only from NIRISS/SOSS.
The carbon abundance can only exceed the oxygen abundance if either carbon is produced (like in AGB stars) or inserted otherwise into the system (e.g.meteoritic influx), or if the oxygen is simply reduced below the original carbon-abundance level.The latter might occur if the planet already starts out with a primordial atmosphere that is oxygen-depleted due to the formation of H 2 O/CO ices inside the planet-forming disk.The accretion of the gaseous envelope from this oxygen-depleted gas and further processing of the atmospheric gas into clouds can then lead to a C/O>1 (Helling et al. 2014).
Compared to the the Si/O, Mg/O, and Fe/O ratios, C/O is not affected as strongly.The lower abundances of Si, Mg, and Fe are more significantly depleted by condensation compared to oxygen (Fig. 4).However, the C/O ratio is enhanced in regions of cloud formation, due to oxygen element depletion from the gas phase through condensation (Fig. 4, bottom).The maximum value suggested by our simulations is C/O=0.75, which also occurs in the upper modelled atmosphere and thus likely in the optically thin atmosphere observed in transmission.We note that also the set of Y and T-type brown dwarfs Guillot et al. (2022) appears to reach super-solar C/O.
A decrease below the undepleted (solar) C/O value of 0.55 would suggest an increased amount of oxygen in the atmosphere or a substantially decrease of the carbon in the atmosphere.Molaverdikhani et al. (2020) show for the ultra-hot Jupiter HAT-P-7b that advection driven non-equilibrium may occur for H 2 O but C/O would be determined by the oxygen depletion by cloud formation.
We note that in this section the local gas-phase C/O is discussed that includes the depletion of oxygen by cloud formation.The retrieval procedure in Section 4.1 addresses the bulk C/O, which is indicative of the formation history of the planet.It is, however, demonstrated in Sect.4.2 that the retrieval suggests two solutions as best (or most probable) fit to the observed data.
The mean molecular weight remains that of an H 2 -dominated gas until deep in the atmosphere because the upper atmosphere temperatures remain > 2000K (see Fig. 1).

Where clouds on WASP-96b get optically thick
Here we compute, based on our microphysical cloud model, the pressure where the optical depth of the clouds reaches unity (p gas (τ(λ) = 1)) -here after referred to as the 'cloud deck'.This definition aligns with the cloud top pressure recovered in However, given the homogeneity of cloud formation in these models of WASP-96b, the effects of slant observing geometry can be corrected for through a correction factor (Fortney 2005).Furthermore, the optical depths derived here depend on the computational domain of the GCM, as above this level we are unable to model cloud formation.Previously Helling et al. (2022) has shown that restricted GCM domains misses the 'mineral haze layer' at higher altitudes.
Taking four key points around the equator (θ = 0.0 • ) of WASP-96b (Fig. 5 top left), we see that for wavelengths < 4µm indeed the cloud deck is consistent.Nonetheless, in this range the cloud deck still shows a substantial slope with wavelength.For example between 0.3...5 µm for the morning terminator, the cloud deck pressure varies between 1.3×10 −4 ...6×10 −3 bar.This suggests flat cloud decks are a poor fit across this wavelength range, which becomes more important as we move towards an era of broader wavelength coverage with a single mission.In the mid-IR, it is a different story, spectral differences may be expected within the silicate feature λ-range, visible with NIRSPec, NIRCam, ARIEL and MIRI.Spitzer sits right on this silicate edge at ∼ 4 µm, which may complicate combining datasets.In particular focusing on the two terminator limbs there is substantial differences in the cloud deck from 4...8 µm of nearly an order of magnitude in pressure.
The high pressure cloud decks (or cloud-free atmospheres) previously retrieved by Nikolov et al. (2018) in order to fit the Na line (dashed blue vertical line) are inconsistent with microphysical cloud formation.At the Na line, the model produces a cloud deck at ∼ 2.5 × 10 −4 bar, whereas pressures as high as 10 −2 bar are needed to fit the pressure broadening observed.However, given the pressure-temperature structures shown in Fig. 1, simple thermal stability arguments will expect clouds to form.Re-trievals cannot simply infer condensable-favouring global temperatures and claim low cloud decks, or cloud-free atmospheres without proposing some mechanism for cloud suppression.
We therefore explore 3 parameters, vertical mixing efficiency, cloud particle porosity, and metallicity as to how they impact the cloud deck pressure level.Reduction of the vertical mixing efficiencies is achieved by increasing the mixing timescale.This is done by the introduction of a 'mixing scale factor' f mix , such that the new mixing timescale is where τ mix is the mixing timescale from the GCM derived the same way as in Helling et al. (2022).Alternatively, using K zz = H p τ mix implies that similarly K zz is reduced by factors of f mix .Reduction of the vertical mixing efficiencies ( f mix = 0.1, 0.01) reduce the cloud deck height by an order of magnitude in pressure, but otherwise maintains the wavelength dependence.The mechanism for such inefficient mixing requires further investigations.
Porosity of cloud particles is included as described in Samra et al. (2020), reducing the effective material density of cloud particles ρ eff s through the porosity factor f por , ρ eff s = ρ s (1 − f por ).Mass conservation ensures that correspondingly the radius of cloud particles are also increased.In addition, the porosity of the cloud particle is incorporated into opacity calculation through inclusion of vacuum refractive index into effective medium theory (Bruggeman 1935;Looyenga 1965), in proportion to the porosity factor of the particle.
Highly porous particles reduce the opacity of the clouds, even with enhanced cloud mass load and a higher onset of bulk growth, as was the case in Samra et al. (2020), although the effect here is substantially larger than in the hotter atmosphere (T eff = 1800 K) presented there.Unlike reductions in mixing efficiency, increased cloud particle porosity affects wavelengths > 4 µm non-uniformly, the cloud deck pressure level at the silicate feature wavelengths is particularly affected by compact vs porous cloud particles.
The porosity and shape of cloud particles (or more generally aerosols) in exoplanet atmospheres remains uncertain (Gao et al. 2021).However, there is some tentative evidence of porous cloud particles (generated through collisions) as an explanation for the flat spectrum of the mini-Neptune GJ-1214b (Ohno et al. 2020).Although, as shown in Samra et al. (2022), collisions in exoplanet atmospheres when including turbulence are expected to be destructive rather than producing larger aggregate particles.However, hazes of fractal aggregates have been inferred for a number of bodies, both Solar system and exoplanets (e.g.Adams et al. 2019).Porosity has also been examined in the other astrophysical contexts, such as protoplanetary disc dust (Woitke et al. 2016), albeit only with a factor of 25% vacuum assumed.Overall whether cloud particles could reach such high (90%) porosities is unknown, but it still remains an effect which is important to consider in interpreting exoplanet observations.
For this work we have so far assumed solar elemental abundances (Asplund et al. 2009) for the gas composition of the atmosphere.Nikolov et al. (2022) derive elemental abundances for the host star WASP-96 using FEROS spectra.These abundances are enhanced with respect to the solar abundances of (Asplund et al. 2009), at most around a factor of ∼ 3x solar (for example Mg and Mn -see their Table 4).However, it is unclear if stellar metallicities correlate with planetary atmosphere abundances (Teske et al. 2019).Nonetheless, we examine the effect of enhanced metallicity for cloud formation in the atmosphere using a 10x Solar metallicity for the planetary atmosphere.This is done by increasing abundance of everything except H and He and renormalising to H abundance.Metallicity for WASP-96b, based on formation arguments for lower mass gas-giant planets (e.g.Carone et al. 2021;Chachan et al. 2019;Schneider & Bitsch 2021), is expected to be super-stellar.The enhanced metallicity increases cloud formation and leads to a higher cloud deck (bottom right of Fig. 5), making it even harder to reconcile with super-solar retrievals.

Clouds affecting 0.3 − 5µm spectra
Solving the radiative transfer problem through a prescribed atmospheric structure with clouds allows us to explore the effect of cloud properties on the possible spectral appearance of WASP-96b within the wavelength range that has so far been explored with VLT/FORS2, HST/WFC3 and Spitzer.This procedure is known as forward modelling within the retrieval community.The ARCiS (Min et al. 2020;Ormel & Min 2019) frame work is applied here.We further investigate the possibility of finding a solution to the VLT/FORS2, HST/WFC3 and Spitzer data including the effect of clouds using a retrieval setup which includes a parameterised cloud formation routine.
Retrieval studies of the HST/WFC3 and Spitzer/IRAC transmission spectra of WASP-43b also previously struggled to find evidence for the presence of clouds in the atmospheric spectra (Kreidberg et al. 2014;Chubb et al. 2020), at least at the pressure levels probed by these observations.However works like Helling et al. (2020) studied the formation of mineral clouds and hydrocarbon hazes in the atmosphere of WASP-43b and demonstrate that it is rather unlikely that WASP-43b is cloud-free.Robbins-Blanch et al. ( 2022) performed a detailed study of the effects of clouds on emission phase curve spectra of WASP43b by computing various forward models for clear and cloudy atmospheres using PICASO (Batalha et al. 2020;Batalha & Rooney 2020).
They found that cloudy phase curves, produced from 3D models, provide much better agreement with the WFC3 and Spitzer nightside data than cloud-free models.

Name Description Priors
Parameters included for both cloud formation and clear retrievals b : We use a gaussian prior on the mass based on radial velocity measurements of WASP-96b (Hellier et al. 2014).We compute the mass based on R p and log(g p ) and then place the prior on this derived parameter.2015), C 2 H 2 (Chubb et al. 2020), CN (Brooke et al. 2014), and SiO (Barton et al. 2013;Yurchenko et al. 2021).We ran two different retrievals: A constrained cloud formation retrieval, following the method of Min et al. (2020), on the combined HST/VLT/Spitzer data set (including HST/VLT offset) of Nikolov et al. (2022), and a completely clear (no parameters describing clouds included) retrieval, for comparison.Both retrievals assume chemical equilibrium for the gas-phase chemistry, which is computed after the cloud formation computations.The parameters retrieved and their priors can be found in Table 1.The best fit retrieved spectra, with 1, 2, and 3 σ shading, can be found in Figure 6 Guillot (2010).These are used to construct the pressuretemperature profiles for the two retrievals given in Figure 6 (right).The pressure-temperature profiles from the GCM used in our atmospheric modelling and cloud formation models are plotted for comparison.These can also be compared to Figure 5 of Nikolov et al. (2022), with the clear retrieval pressuretemperature profile in our Figure 6 (right) appearing very similar to their retrieved pressure-temperature profile, which also assumes equilibrium chemistry.
The two panels of Figure 6 illustrate solutions from the cloud formation retrieval based on slightly different model outputs.The 'best fit' model is the one with the lowest χ 2 when compared to the observations.However, there are other possible solutions which can also fit the data similarly well.The 'Median Probability Model' (MPM) (Barbieri & Berger 2004;Barbieri et al. 2021) is based on the posterior distributions shown in Fig- ure B.1, and represents an alternative region of parameter space for the cloud formation retrieval presented in this work.In particular, the posterior distribution for the retrieved C/O has two potential solutions; one with high C/O (∼0.78; the 'MPM' solution), and one with lower C/O (∼0.15; the 'best fit' solution).It can be seen from the left panel of Figure 6 that the retrieved best fit and MPM spectra are very similar in the region covered by the VLT and HST observations, and fit the observed data similarly well.They deviate, however, in the region of the spectra for wavelengths longer than the HST data.The same is true for the clear retrieval; it deviates most in the region not covered by the presented observed data.It is therefore difficult to distinguish them based on the presented observed data.The bottom two panels of Figure C.1 give the molecular abundances (volume mixing ratios) as a function of pressure for the MPM (left) and best-fit (right) solutions.There are large variations in the CO and CH 4 abundances in particular between these two solutions.The primary spectral features of both these molecules are not in the region of the spectra covered by the HST and VLT observations; around 2.3 and 5 µm for CO, and around 3.3 µm for CH 4 .As illustrated by Gasman et al. (2022), other absorption features of CH 4 which are present in the region of the HST data are often masked by stronger H 2 O features.This highlights the benefits of the JWST observations of WASP-96b (Pontoppidan et al. 2022), which do cover this region.
It can be seen by the right-hand panel of Figure 6 that the cloud formation retrieval solution with the lower C/O has a higher atmospheric temperature than the solution with the higher C/O.This is due to constraints placed by the cloud formation scheme.The same does not apply to the clear atmosphere retrieval; there is no constraint on the temperature based on what types of clouds are expected to form and a lower temperature is therefore preferred.The cloud formation model with the lower C/O necessitates a higher atmospheric temperature at the pressure layers being probed by the observations.This is because a low C/O and therefore an abundance of oxygen would mean that K-and Na-silicate clouds are expected to form if the temperature is low enough.However, the strong gas-phase Na and K spectral features suggest that these species of clouds either have not formed efficiently in this atmosphere, at least not at a very large abundance to take all the Na and K away from the gas-phase, or they form but they are optically thin enough such that the gasphase Na and K is still observable.A recent observation of planetary mass companion VHS 1256-1257 b using JWST's NIR-Spec IFU and MIRI MRS modes revealed solid-states silicate features in the longer wavelength part of the atmospheric spectrum (Miles et al. 2022).They conclude this as strong evidence for small silicate particles in the brown dwarf's atmosphere; as demonstrated by works such as Min et al. (2004), small particles lead to less overall extinction across the spectra, but a more prominent silicate feature.The atmospheric mixing is described by a K zz = 10 8 . . . 10 9 cm 2 /s.These observations suggest that silicate clouds do form in sub-stellar hot atmospheres and it is possible to directly observe their spectral signatures with JWST.This is either due to a lower C/O and therefore less oxygen present for these clouds to form from, or a lower C/O but a high enough atmospheric temperature to stop them from forming in abundance.We also ran a retrieval with cloud formation, but without the cloud species which condense Na and K included.In this case, the lower C/O solution was preferred, as there were no constraints placed based on the temperature of formation for these types of clouds.It should be noted that the C/O for the cloud formation retrievals is the bulk C/O, i.e. before clouds are allowed to form.Clouds will generally deplete oxygen from the atmosphere, so the retrieved C/O without cloud formation constraints (but where clouds should be present based on atmospheric temperature, such as for the clear retrieval), will be expected to be retrieved higher than the bulk C/O.The models with the lowest C/O ratio show significant H 2 O features in the region not covered by the HST and VLT data, which is not the case for the high C/O case.Since the newly observed JWST ERO data appears to show large H 2 O spectral signatures (Pontoppidan et al. 2022), a situation with a lower C/O and higher atmospheric temperature could be preferred.A more detailed study of this JWST spectra of WASP-96b in the near-future should be very informative.
The reduced-χ 2 value was very similar for the clear and the cloud formation retrievals (1.08 and 1.09, respectively), indicating that both can fit the observed data equally well.The cloud formation retrieval is however not considered statistically significant in comparison to the clear retrieval, i.e. the inclusion of extra parameters to characterise the clouds is not justified by a significant increase in how well the model fits the data.This has previously led to the conclusion that the atmosphere must be clear.Based on expectations from detailed simulations of cloud formation in the atmosphere of WASP-96b, as demonstrated in this paper, it would be very surprising from a physical point of view if the atmosphere is completely clear.We therefore try to offer some explanations to confront this seeming contradiction between what models predict (a cloudy atmosphere) and what the observed data appears to show (a clear atmosphere).
Firstly, in Figure B.1, the cloud diffusivity K zz (which represents the vertical mixing strength of cloud particles) is constrained to be no higher than around 2×10 6 cm 2 s −1 .This corresponds to the reduced mixing scenario explored in Figure 5 (top right).A reduction in vertical mixing of cloud particles settles the cloud layer further down in the atmosphere, thus reducing its optical effects on blocking or muting gas-phase molecular or atomic signatures.We further illustrate this by computing a forward model using our best fit cloud formation retrieval, with varying values of K zz (Figure 7). Figure 8 shows the pressure level where the optical depth of the atmosphere reaches unity (i.e. it becomes optically thick) as a function of pressure for the models.It can be seen from Figure 7 that lower values of K zz , around 2×10 6 cm 2 s −1 , allow for a clear detection of the Na doublet at around 0.6 µm, with the atomic line broadening by H 2 and He clearly visible.Increasing K zz mutes this, and other gasphase, spectral features.Our retrieved value of vertical cloud diffusion parameter K zz is between 1.7 × 10 5 and 1.8 × 10 6 cm 2 s −1 , within 1 σ bounds, and between 2.1 × 10 4 and 2.5 × 10 7 cm 2 s −1 within 3 σ bounds.The priors on K zz were 10 5 ...10 12 cm 2 s −1 , so these values are at the lower edge of the priors.In comparison the K zz values used in our microphysical model range between 10 7 ...10 9 cm 2 s −1 , values comparable to the directly imaged brown dwarf result from Miles et al. (2022).Hence with a reduction by a factor of 10 −2 , we find some agreement between the retrieved K zz and a microphysical model.Although, the deeper cloud deck (as shown in Section 3.4) still is higher than shown in Figure 8.We again re-iterate that the mechanism for such reduced mixing is uncertain.
Figure B.1 shows the free parameter, the integrated cloud particle nucleation rate Σ [g cm −2 s −1 ] is not particularly wellconstrained by the retrieval, but a lower value is preferred.The nucleation rate, J i [cm −3 s −1 ] for i=TiO 2 , SiO, KCl, NaCl, in our kinetic cloud formation model is not a parameter, but rather we calculate the local nucleation rate, (see Section 2) according to the local thermodynamic properties (T gas , p gas ).The integrated cloud particle nucleation rate can be calculated a Σ = i m i J i dz. (2) Where the summation is over the four nucleation species considered (i =TiO 2 ,SiO,NaCl,KCl), with the mass of individual monomer species m i and the respective nucleation rates J i .The integration depends on the computational volume of the GCM.The cloud extension is hence, a hidden within the Σ parameter.Using this, the terminator profiles (both morning and evening) from the GCM give integrated nucleation rates, of between 1.49 × 10 −15 ...4.56 × 10 −14 g cm −2 s −1 .This is in broad agreement with the retrieved values for the integrated nucleation rate, although the retrieval does favour slightly lower values 5 × 10 −17 to 6 × 10 −14 g cm −2 s −1 , within 1 σ bounds.Examining Figure B.1, we see that the lower retrieved values of Σ are towards the bottom of the prior region (lower limit of 10 −17 g cm −2 s −1 ).This is because the retrieval is essentially trying to minimise the formation of clouds in an atmosphere that, as the kinetic cloud formation model shows, should form clouds.
Another factor which could influence the observed spectra is the porosity of the cloud particles.The cloud particles in our retrievals are all assumed to be non-porous (porosity = 0).Figure 7 shows, however, that increasing porosity of cloud particles acts in a similar way to reducing the vertical mixing of cloud particles.The clouds become more transparent, reducing their optical thickness as porosity is increased.This particularly affects wavelengths shorter than around 2 µm.This trend agrees with the deeper cloud deck pressure found in Section 3.4 when including the effects of cloud porosity.However the effect here is greatly reduced when compared to Figure 5 (lower left), due to the already deeper cloud in the best fit constrained model.The effect of cloud particle porosity has been explored further in works such as Samra et al. (2020).The plausibility of highly porous cloud particles has been discussed in Section 3.4.
Figure 9 shows the average particle size retrieved by ARCiS for the 'best fit' and 'MPM' solutions as well as the surface average particle size from the kinetic cloud formation model (Appendix A, Eq.A.6 Helling et al. 2020).The average cloud particle size from the kinetic cloud formation are larger than the retrieved average particle sizes at all pressures in the atmosphere.The kinetic cloud formation model experiences condensational growth for the cloud particles right from the top of the GCM domain 10 −5 bar, and hence larger particle sizes.This contrasts with the results from ARCiS, where for both solutions growth of average particle size does not start until deeper in the atmosphere.Both the kinetic cloud formation model and the ARCiS parameterised model use cloud condensation nuclei (CCN), the first step of cloud formation of size ∼ 10 −3 µm.
Figure 6 demonstrates that the retrieval result both (constrained and un-constrained) have large uncertainties in the retrieved pressure-temperature profiles.A temperature inversion could indicate clouds present deeper in the atmosphere, acting to heat the upper layers and thus making them too hot for further cloud formation.This offers another explanation for why the atmosphere appears relatively cloud free from observations, but could still contains an abundance of clouds as predicted by our kinetic models.

Discussion
The compatibility of combining ground-and space-based data has been questioned in the case of WASP-96b by Yip et al. (2020a), who find a significant offset between HST and VLT observations.The offset was also addressed by Nikolov et al. (2022), who performed a combined analysis of various observations (VLT, HST, and Spitzer), and also determine an offset between the HST and VLT observations.This dataset (including their offset) is the one used in the retrievals of the present work.Yip et al. (2020b) further demonstrate the potential issues with combining data from instruments such as HST/WFC3 and Spitzer/IRAC which have no overlapping wavelength coverage.As the Spizter data photometry points come with large uncertainties, use of them in retrievals to quantify abundances of species such as CO or CO 2 can lead to biases in retrievals, particularly when constraining elemental abundances such as C/O.This has been highlighted previously by works such as Irwin et al. (2020).We therefore caution our retrieved value of C/O and await reduced spectra observed using JWST for further analysis.Our results which include cloud formation tentatively show two potential regions of parameter space which fit the observed pre-JWST data well: one with high C/O and a lower atmospheric temperature, and one with low C/O and a higher atmospheric temperature (in the region probed by the observations).The presence of large H 2 O features in the region covered by the JWST observations could push towards the low C/O with higher atmospheric Best fit clear retrieval K zz = 1x10 6 cm 2 s 1 K zz = 2x10 6 cm 2 s 1 K zz = 5x10 6 cm 2 s 1 K zz = 1x10 7 cm 2 s 1 Observed spectra K zz = 5x10 6 cm 2 s 1 , porosity = 0 K zz = 5x10 6 cm 2 s 1 , porosity = 0.5 K zz = 5x10 6 cm 2 s 1 , porosity = 0.9 Observed spectra  Pressure (bar) Best fit clear retrieval Best fit constrained retrieval K zz = 1x10 6 cm 2 s 1 K zz = 2x10 6 cm 2 s 1 K zz = 5x10 6 cm 2 s 1  temperature solution.We note that the best-fit analysis and retrieval of Nikolov et al. (2018) find a sub-solar C/O in both cases.These retrieved C/O are bulk C/O (i.e. the C/O available before clouds are allowed to form), and therefore a retrieved C/O without taking cloud formation into account will be higher than this bulk value.The bulk C/O is useful to be able to make links to planet formation; a parameterised planet formation model such as that of Khorshid et al. (2021) can be coupled with retrieval codes like ARCiS.In Figure C.1 we show the cloud species formed by the two solutions found in our retrieval in the two upper panels.On the left column we show the 'best fit solution' with C/O = 0.15 and a higher atmospheric temperature structure, and on the right the 'MPM' solution with C/O = 0.78 and a lower atmospheric temperature structure.The corresponding retrieved volume mixing ratios (VMR) of the gas-phase molecular species are given in the panels below.In the 'MPM' case, the high C/O (planetary bulk C/O) indicates not much oxygen available to form clouds. Condensates such as Fe, SiO 2 , Al 2 O 3 (corundum), TiO 2 and VO are able to form, and the Na/K silicate and Mg/Ca silicate condensates do not have enough oxygen left to be able to form even though the temperatures would allow them to do so.The 'best fit' case, on the other hand, has an abundance of oxygen in the atmosphere (low C/O).This allows the lower-temperature silicate clouds to form.The higher temperature profile does mitigate how much can form to some extent.Even with the Na/K silicate clouds that do form, it can be seen from the bottom left panel of Figure C.1 that there is still an abundance of gas-phase Na and K available to produce the strong atomic features seen in the observed spectra.

Conclusion
Based on our hierarchical modelling approach that combines a 3D GCM atmosphere solution for WASP-96b with a kinetic, non-equilibrium formation model for mixed-material cloud particles, we suggest that WASP-96b is not cloud-free as previously instigated based on retrieval approaches.This suggestion is supported by applying the ARCiS retrieval framework to the same HST, Spitzer and VLT data which shows that cloudy solutions do reproduce the observed spectra.More than one cloudy solution provides good fits to the observational data such that retrieved solutions may not be unique, requiring more physical input which is an important step towards discussing retrieval results from present and future JWST data.Clouds in WASP-96b would cause the following effects within the JWST wavelength range: -The cloud top varies with wavelength within the JWST NIR-Spec and NIRISS spectral ranges for at least 1 orders of magnitude in pressure.-Clouds become optically thick at different pressures dependent on wavelengths.To achieve optically thin clouds down to p = 10 −2 bar, as implied by the sodium feature at ∼ 0.6 µm, a moderate mixing efficiency is required.-The long wavelength end of NIRSpec and short end of MIRI may probe atmospheric asymmetries between the limbs of the terminator on WASP-96b.-WASP-96b could only be cloud free in the unlikely case of a truly static atmosphere, and only cloud free in the region probed by observations in the case of extremely low mixing efficiency or if a temperature inversion confines the clouds to higher pressures.

D
Fig. 2: WASP-96b cloud maps in 2D terminator slices as a result of the microphysical cloud model: Top Left: Local atmospheric gas temperature and gas pressure (T gas , p gas ), Top Right: Total mineral seed formation (nucleation) rate, J * = i J i [cm −3 s −1 ] (i=TiO 2 , SiO, NaCl, KCl).Bottom left: Dust to gas ratio ρ d /ρ.Bottom right: Mean cloud particle radius a [µm].

Fig. 4 :
Fig. 4: Gas phase abundance ratios, for morning (φ = −90.0,dashed lines) and evening (φ = 90.0,solid lines) terminators.Some morning/evening terminator differences exists.Top: The mineral ratios Si/O, Mg/O, Fe/O are strongly affected by element depletion due to cloud formation.Bottom: C/O ratio appears lesser affected but oxygen is two orders of magnitude more abundant in a solar set of element abundances.C/O is therefore shown on a linear scale, where C/O= 0.55 is the Solar values.
lit. value a , units of R J , flat linear prior log 10 (g p ) Base-10 logarithm of the planet surface gravity b 1 to 5 (with g given in cgs units) 4 cm 2 g −1 , flat log prior T int Temperature at an optical depth τ = 2 3 as caused by internal heat from the planet 10 to 3000 K, flat log prior Parameters included in cloud formation retrieval only K zz Cloud diffusion coefficient 10 5 -10 12 cm 2 s −1 , flat log prior log 10 Σ Nucleation rate 10 −17 -10 −7 g cm −2 s −1 , flat log prior a : Hellier et al. (2014)

Fig. 6 :
Fig. 6: The best fit models (1, 2 and 3 σ bounds shading) retrieved with ARCiS from the Nikolov et al. (2022) pre-JWST observations.In purple is the ARCiS retrieval with no clouds included, and in orange with cloud formation included (as described in Ormel & Min 2019).Both retrievals assume equilibrium chemistry.In the cloud formation case the retrieved parameters are shown in Figure B.1, and for the no cloud case they are shown in Figure B.2. MPM refers to the 'Median Probability Model', which is slightly different to the 'best fit' model; see text.Left: Retrieved spectrum, Right: the retrieved pressure-temperature profile, for reference the GCM temperature-pressure structures are over-plotted, using the same colour scheme as Fig. 1.
(left).The observed spectra is plotted for comparison.Our retrieved posterior distributions can be found in Figure B.1 for the cloud formation case, and Figure B.2 for the clear atmosphere case.The first four parameters in Figures B.1 and B.2 describe the pressure-temperature parametrisation of Figure C.2 shows the volume mixing ratios of various molecules predicted to be present at the equator of the morning (left) and evening (right) terminators, based on the GCM and kinetic cloud models of this work.CO is expected to be abundant in both cases, along with H 2 O at most altitudes.These are more similar to the lower C/O retrieval solution in the bottom left panel of Figure C.1.It can also been seen from Figure C.1 that the difference in H 2 O abundance between the two retrieval solutions can have a noticeable effect on the spectra, again in particular in the region above 1.7 µm.

Fig. 7 :
Fig. 7: The best fit synthetic spectrum from the constrained retrieval of the Nikolov et al. (2022) data (VLT, HST, and Spitzer) with varying values of cloud diffusivity, K zz , (top), and for varying cloud particle porosity for atmospheric models with cloud diffusivity K zz = 5×10 6 cm 2 s −1 (bottom).

Fig. 8 :
Fig.8: The same cloud formation setups with varying cloud diffusivity K zz as for Figure7, but this time showing the pressure levels where optical depth τ reaches unity, as a function of wavelength.The sharp drop of the Na lines (green line) is artificial and caused by the cut-off in the opacity data available.

Fig. 9 :
Fig.9: Average particle sizes from the constrained retrieval with ARCiS for the MPM (yellow) and the best fit (blue) models.Over-plotted are surface averaged particle sizes a A results from the kinetic cloud formation model (grey), for the evening (φ = 90.0• , solid) and morning terminators (φ = −90.0• , dashed).