© The Authors 2024

1 Introduction

Transmission and emission spectra are powerful diagnostics to obtain the composition of a planet atmosphere. From the wavelength and strength of the spectral features, the presence of a variety of molecules has been identified and their abundance has been constrained. With the Hubble Space Telescope (HST) and Spitzer, one of the most well-observed molecules in the exoplanet spectrum is H₂O (e.g., Huitson et al. 2013; Deming et al. 2013). Before the James Webb Space Telescope (JWST), the detection of carbon-bearing molecules, such as CO and CH₄, were mostly achieved by emission spectroscopy and limited to giant planets (Konopacky et al. 2013; Barman et al. 2015; Sheppard et al. 2017). However, the constraint on the abundance of C- or O-bearing species aside from H₂O was limited by wavelength coverage and sensitivity (Madhusudhan 2019). With its unprecedented spectral resolution, JWST has opened a new era for atmospheric spectrum characterization. Covering the spectral range from 0.6 μm to 15 μm, JWST has robustly observed CO (e.g., Rustamkulov et al. 2023), CO₂ (e.g., Grant et al. 2023; Xue et al. 2024; Alderson et al. 2023), and CH₄ (e.g., Bell et al. 2023) on giant planets using transmission spectroscopy. Moreover, JWST has extended our spectroscopic understanding to super-Earth planets. It detected signals of CO₂ and CH₄ on hycean candidates TOI-270 d (Holmberg & Madhusudhan 2024) and K2-18 b (Madhusudhan et al. 2023). For terrestrial planets, the emission spectrum of 55 Cnc e revealed a feature from CO or CO₂ (Hu et al. 2024). Besides gas features, solid phase silicate features have been discovered on hot Jupiters such as WASP-107 b (Dyrek et al. 2024) and WASP-17 b (Grant et al. 2023) for the first time.

The composition of an exoplanet is the key indicator of its formation location and evolution pathway. This is because the composition of the building blocks that the planet has accreted varies with the local physical conditions (temperature and pressure) in the protoplanetary disk. For example, the location of the snowline regulates the existence of volatiles and the content of the corresponding vapor (Öberg et al. 2011). The ratio and composition of solid and gas material that a planet accretes determines the bulk properties of the planet, such as the metallicity, C-to-O ratio (C/O), and volatile ratio (Madhusudhan et al. 2014; Booth et al. 2017; Schneider & Bitsch 2021; Turrini et al. 2021; Khorshid et al. 2022; Danti et al. 2023). In particular, the metallicity and C/O are identified as important tracers of planet formation history, because C- or O-bearing molecules are most commonly constrained in atmospheres. In addition, N/O or S/N offer valuable indicators to constrain a planet’s formation history (Turrini et al. 2021; Khorshid et al. 2024). JWST has detected SO₂ on a giant exoplanet atmosphere (Alderson et al. 2023; Dyrek et al. 2024), providing the opportunity to trace back their formation (Schneider & Bitsch 2021; Crossfield 2023; Khorshid et al. 2024). This is because the main sulfur carrier, FeS, condenses at a temperature higher than carbon or oxygen (650 K) (Kama et al. 2019), which indicates that S has been accreted in solid form.

Inferring the chemical makeup and the abundance of an atmosphere from spectroscopic data is a process referred to as “retrieval”. In retrieval a parameter-dependent forward model is computed to obtain a synthetic spectrum that is compared to the observational data. When performing retrieval, Bayesian analysis such as Markov-chain Monte Carlo (MCMC, e.g., Madhusudhan et al. 2011; Benneke & Seager 2012) or the nested sampling method (e.g., Line & Parmentier 2016; Lavie et al. 2017) are applied on the observed spectrum, to test how different sets of parameters fit the data. Bayesian analysis is done by running the forward model thousands to millions of times (e.g., Chubb & Min 2022) to obtain the probability distribution of each parameter. As the parameter space is highly multidimensional and the forward computation is repeated many times, the CPU runtime of each forward run must be short, preferentially limited to the order of seconds. There are already a dozen of retrieval codes to date (MacDonald & Batalha 2023), for example NEMESIS (Irwin et al. 2008), TauREx (Al-Refaie et al. 2021), POSEIDON (MacDonald & Madhusudhan 2017), ATMO (Wakeford et al. 2017), petitRADTRANS (Mollière et al. 2019), ARCiS (Min et al. 2020) and PICASO (Mukherjee et al. 2021), with varying assumptions and physical complexity. This complexity occurs in many aspects of the forward model, such as the gas-phase chemistry, temperature profile, or molecular opacity model. For example, free-chemistry, where molecular abundances are free parameters, and equilibrium chemistry, in which the Gibbs free energy is at its minimum, are two common choices for gas-phase chemistry.

A particular aspect atmosphere forward models need to address is the treatment of clouds. Clouds are ubiquitously found on exoplanets and brown drawfs (e.g., Pont et al. 2008, 2013; Kreidberg et al. 2014). Clouds affect the spectrum in a variety of ways. First, because they obscure the deep atmosphere from view, they tend to smooth out molecular features. At near-IR-tovisible wavelengths, small cloud particles or hazes high up in the atmosphere can lead to scattering slopes (Pont et al. 2013; Sedaghati et al. 2017). Second, at mid-IR, clouds display distinct solid features (e.g., the 10 μm silicate feature; Miles et al. 2023; Dyrek et al. 2024). Third, the formation of clouds consumes condensable vapor, which thereby depletes the corresponding atoms in the atmosphere (Komacek et al. 2022; Lee 2023). The change in the molecular inventory is reflected in the spectra by weaker molecular features. Moreover, the presence of clouds can also change the heat transfer, influencing the temperature structure in the atmosphere (Roman & Rauscher 2019; Harada et al. 2021; Roman et al. 2021). Therefore, unless atmospheres are undeniably cloud-free, understanding clouds is essential to interpret exoplanet spectral data.

Most forward models in retrieval codes parameterize clouds (see Barstow 2020 for an overview of cloud parameterizations). In its simplest form, this entails a cloud deck pressure below which the atmosphere is optically thick (Benneke & Seager 2012). Some models design a parameterized function representing cloud opacity, while others impose monodisperse particles and calculate Mie scattering (Lavie et al. 2017). Importantly, in these models the properties of the cloud layers are not limited by physical principles. A more physically motivated cloud model is proposed by Ackerman & Marley (2001) (AM01 hereafter). This model assumes that all saturated material condense onto cloud particles, and solves for the condensate concentration by balancing the upward diffusive motion with the raindrop sedimentation. The AM01 cloud model is used extensively in retrieval codes, for example, petitRADTRANS and PICASO. While popular, the AM01 model relies on adhoc parameterizations. The particle sedimentation velocity is prescribed with an f_sed parameter and the particle sizes are not solved for.

Other cloud models are physically motivated. DRIFT (Helling & Woitke 2006; Helling et al. 2008) solves the moment equations of the particle size distribution. Nucleation, condensation, coagulation, and settling are modeled in detail. Further complication comes from computing the time-dependent particle size distribution. This approach is adopted by CARMA (Turco et al. 1979; Toon et al. 1979), and is applied to H₂SO₄ clouds on Venus and silicate clouds on hot Jupiters (Powell et al. 2018). Recently Kiefer et al. (2024a) studied time-dependent cloud formation under gas phase disequilibrium, and found nucleation and cloud formation can happen within 1 second. Though these models consider a variety of physical processes, and promote our understanding on exoplanet clouds, they consume significantly more computational resources than parameterized cloud models. Therefore, forward model with high level of physical complexity is usually ran on a limited number of parameter grids (e.g., Mukherjee et al. 2021), rather than in the MCMC or nested sampling simulations.

The ARCiS atmosphere model (Min et al. 2020; Ormel & Min 2019) attempts to balance physical consistency with computational efficiency. The philosophy of ARCiS is to model physical processes that affect the cloud formation most, while parameterizing processes for which we have a poor understanding or those that are too time-consuming for a robust retrieval. This approach has been successfully applied to retrievals using self-consistent cloud formation (Min et al. 2020), disequilibrium chemistry (Kawashima & Ikoma 2019), and 3D phase curve retrievals (Chubb & Min 2022). However, one aspect where ARCiS (and other) cloud models can be improved on is the multispecies nature of clouds.

It is indeed expected that the composition of the cloud particles changes as function of atmospheric depth along with the changing molecular inventory and physical conditions (temperature, pressure). Therefore, an efficient cloud model that simulates, to first order, the realistic condensation and transport process taking place in the atmosphere would be a great asset to the community. To this effect, we develop a multispecies cloud model, ExoLyn, that computes cloud condensation, sedimentation, and diffusion in a computationally efficient, yet physically consistent way. The input for ExoLyn is a temperature-pressure (T-P), nucleation profile, cloud forming reactions and vapor abundances at the lower atmosphere boundary. Then ExoLyn will compute the cloud concentration and components by solving the transport equation including condensation, sedimentation, and diffusion. ExoLyn has been designed to study cloud formation across a variety of planets: hot Jupiters, sub-Neptunes, and self-luminous planets. ExoLyn is publicly available, written mainly in Python, and compatible with popular retrieval frameworks such as petitRADTRANS.

This paper is structured as follows. In Section 2 we review the ARCiS cloud model and extend it to multiple condensate species. In Section 3 we present the ExoLyn results for hot Jupiter planets. Our results are benchmarked against previous models in Section 3.3. The resulting transmission spectra are shown in Section 3.5. The application of ExoLyn to self-luminous planets and sub-Neptunes are demonstrated in Section 4 to show the code’s flexibility. In Section 5, ExoLyn’s performance is assessed and its implication to future atmospheric retrieval is discussed. Finally we summarize the paper in Section 6.

2 Method

ExoLyn extends the single species cloud model by Ormel & Min (2019) (OM19 hereafter). In ExoLyn, a cloud particle is composed of one or more nuclei and multiple condensate species. The condensates form following actual chemical reactions, which we describe in detail in Section 2.2. In Section 2.1 we first introduce clouds made of nuclei and a single condensate species.

2.1 Single condensate species cloud

We briefly introduce the OM19 1D cloud model which is dedicated to a single condensate cloud species. OM19 considered a cloud particle composed of a nucleus, coated by single-species condensate formed from corresponding hypothetical vapor (e.g., MgSiO₃ vapor). The condensate is formed on the nuclei according to the saturation ratio, S = p_vap/p_sat, where p_vap is the vapor partial pressure and p_sat is its saturation pressure. Besides condensation, particles also grow by colliding and sticking with other particles. After formation, cloud particles are transported through the cloud by settling and turbulence, characterized by a turbulent diffusion coefficient K_zz. The combination of condensation, diffusion, settling, and coagulation is solved for a steady solution, where all processes are in equilibrium. Let x_c and x_v to be the mass mixing ratio of the condensate and vapor species, the transport equations for condensate and vapor are $$-\frac{\mathrm{\partial }}{\mathrm{\partial }z}\left(K_{zz}{\rho }_{\text{gas }}\frac{\mathrm{\partial }{x}_c}{\mathrm{\partial }z}\right)+\frac{\mathrm{\partial }}{\mathrm{\partial }z}\left(x_c{\rho }_{\text{gas }}{v}_{\text{sed,p}}\right)={\mathcal{S}}_c$$(1) $$-\frac{\mathrm{\partial }}{\mathrm{\partial }z}\left(K_{zz}{\rho }_{\mathrm{g}\mathrm{a}\mathrm{s}}\frac{\mathrm{\partial }{x}_v}{\mathrm{\partial }z}\right)=-{\mathcal{S}}_c.$$(2)

In the equations above, ρ_gas is the total gas density. The left hand side contains the transport terms and the right hand side of Eq. (2) is the source or sink term for condensates and vapor, respectively. The particle sedimentation (gravitational settling) velocity v_sed,p is assumed to be the terminal velocity according to the Epstein drag law and the Stokes drag law (Weidenschilling 1977), that is, v_sed,p = a_pρ_•g max{1, 4a_p/9λ}/v_thρ_gas, where a_p is the particle radius, ρ_• is the particle density, g is gravitational acceleration, λ is the gas mean free path, and v_th is gas thermal velocity.

At each location, we assume the particles have a single size. In addition to the solid mixing ratio, the particle number density is required to calculate the particle size, which determines the sedimentation velocity and the condensation rate 𝒮_c. Therefore another equation for the particle number density n_p is needed to supplement the set of Equations (1) and (2). To keep all of the unknown variables dimensionless, OM19 solved for the “nuclei mixing ratio” defined as particle number density times nucleus mass normalized by gas density, x_n = n_pm_p0/ρ_gas, where m_p0 is the mass of the nuclei particles. Besides diffusion and settling, collision among the nuclei results in a sink term for their number density: $$-\frac{\mathrm{\partial }}{\mathrm{\partial }z}\left(K_{zz}{\rho }_{\mathrm{g}\mathrm{a}\mathrm{s}}\frac{\mathrm{\partial }{x}_n}{\mathrm{\partial }z}\right)+\frac{\mathrm{\partial }}{\mathrm{\partial }z}\left(x_n{\rho }_{\mathrm{g}\mathrm{a}\mathrm{s}}{v}_{\mathrm{s}\mathrm{e}\mathrm{d},\mathrm{p}}\right)=({\mathcal{S}}_n-\frac{x_n{\rho }_{\mathrm{g}\mathrm{a}\mathrm{s}}}{t_{\text{coag }}}),$$(3)

where 𝒮_n is the nucleation rate and t_coag is the coagulation time scale. The nucleation rate profile depends on the nuclei species (Gail et al. 1984; Lee et al. 2015), curvature of nuclei (the socalled Kelvin effect), and temperature profile (Lee et al. 2018) and is highly uncertain (see review by Helling 2019). Thus, instead of deriving 𝒮_n from first principles, OM19 took a prescribed log-normal distribution along the z direction, where the location of nucleation maxima, the spread of log-normal distribution, and the total amount of nucleation are prescribed as parameters (see Table 1). With solid and nuclei mass mixing ratio in hand, the mass for a single particle is m_p = m_p0 × x_c/x_n. The particle radius is found from its mass and internal density ρ_•, assuming a spherical shape. Thereafter, the particle coagulation time scale t_coag can be calculated with Equation (12) of OM19: $$t_{\text{coag }}^{-1}=\frac{1}{2}\pi n_p{\left(2a_p\right)}^2\mathrm{\Delta }v+2\pi n_p{a}_{p}min(v_{\mathrm{B}\mathrm{M}}{a}_p,D_p),$$(4)

where $v_{\mathrm{B}\mathrm{M}}=\sqrt{16k_{\mathrm{B}}T/\pi m_p}$ is the brownian motion velocity, D_p = k_BT/6πν_molρ_gasa_p, and Δv = v_sed,p/2 comes from differential settling (Okuzumi et al. 2012). v_mol is the molecular viscosity.

For boundary condition, OM19 fixed the solid and nuclei concentration at the lower boundary to be zero, and the vapor concentration to be the envelope bulk concentration, as every cloud species should sublimate for hot planet with substantial envelopes at sufficient depths. For the upper boundary, OM19 fixed the solid, vapor, and nuclei fluxes across it to be 0, as the planet is assumed to be isolated and static, without material flowing into or escape from it.

The single species cloud model in OM19 was generalized to multiple species by Min et al. (2020). This was done by running the cloud formation model once at a time for each cloud species, and iterating over all cloud species. However the problem arises that the sequence of iteration may matter in the final cloud concentration, because the species condensing first has access to all gas phase molecules, leaving less reactant for the next condensate. In this work we overcome this problem by calculating the condensation rate for all condensate species in a thermal chemical approach, which is described in Section 2.2.

Fig. 1 Comparison of ExoLyn’s multispecies approach to condensation against the single-species approach of Ormel & Min (2019). ExoLyn considers the actual vapor species that are present in the atmosphere. As illustrated in the right panel, the condensation rate of MgSiO3 is determined by the amount of Mg – the species that limits the condensation by its vapor concentration and stoichiometric number in the reaction.

2.2 Cloud with multiple condensates

We extend the cloud formation model to multiple cloud components. Instead of assuming the presence of hypothetical vapor in the gas phase (e.g., MgSiO₃ vapor), we consider the realistic chemical reaction from which multiple vapor species condense. The difference between the condensation reaction in our work and OM19 is shown in Figure 1, where we take the formation reaction of MgSiO₃ as an example, $$\mathrm{M}\mathrm{g}+\mathrm{S}\mathrm{i}\mathrm{O}+2{\mathrm{H}}_2\mathrm{O}⟶{\mathrm{M}\mathrm{g}\mathrm{S}\mathrm{i}\mathrm{O}}_3(\mathrm{s})+2{\mathrm{H}}_2.$$(5)

Let subscript i, j, and k denote the i^th solid phase species, the j^th vapor species, and the k^th chemical reaction, the general form of which is $${\displaystyle \sum _{\text{all vapor }j\text{ involved in }k}}{\nu }_{k,j}{V}_j\to S_{i(k)},$$(6)

where v_k,j is the stoichiometric number of vapor species V_j in the k-th reaction. In Equation (6) we have normalized the stoichiometric numbers with respect to that of the solid. Thus v_k,j can be fractional. With this notation the stoichiometric number of any vapor product is negative, for example ${\nu }_{k,{\mathrm{H}}_2}=-2$ for Equation (5). The transport equations for condensate and vapor can be generalized as $$-\frac{\mathrm{\partial }}{\mathrm{\partial }z}\left(K_{zz}{\rho }_{\text{gas }}\frac{\mathrm{\partial }{x}_{c,i}}{\mathrm{\partial }z}\right)+\frac{\mathrm{\partial }}{\mathrm{\partial }z}\left(x_{c,i}{\rho }_{\text{gas}}{v}_{\text{sed,p}}\right)={\displaystyle \sum _k}{\alpha }_{k,i}{\mathcal{S}}_{c,i(k)}$$(7) $$-\frac{\mathrm{\partial }}{\mathrm{\partial }z}\left(K_{zz}{\rho }_{\text{gas }}\frac{\mathrm{\partial }{x}_{v,j}}{\mathrm{\partial }z}\right)=-{\displaystyle \sum _k}\frac{v_{k,j}{\mu }_j}{{\mu }_{i(k)}}{\mathcal{S}}_{c,i(k)},$$(8)

where μ is the molecular weight of vapor (μ_j) or condensate (μ_i), i(k) is the condensate formed in reaction k, and 𝒮_c,i(k) is the condensate i formation rate of reaction k. When reaction k does not contribute to condensate i, α_k,i is 0, otherwise α_k,i is 1. The nuclei mixing ratio x_n, proportional to the cloud particle number by definition, is still described by Equation (3) and the mass of a single particle is updated according to m_p = m_p0 × Σ_ix_c,i/x_n. We note that we do not follow the concentration of gas which is abundant in the atmosphere, for example, H₂, He, or which does not participate in cloud formation, such as CO. Their combined concentration is defined as x = 1, with mean molecular weight m_gas. With this assumption, the production or consumption of vapor does not change m_gas.

Particles grow because vapor condenses on their surfaces. For a given condensate species, there exists one vapor molecule that limits the growth of the condensate (see Figure 1), which we refer to as the “limiting vapor species”, denoted as j* for reaction k. Take, for example, the formation of MgSiO₃ (Figure 1). If the impinging rate of Mg is 1 per unit time, SiO is 2, and H₂O is 3, then the reactant limiting the reaction rate is Mg. Formally, the limiting species is determined by finding the vapor species j that minimizes {x_v,jv_j/μ_jν_k,j}, for all vapor species taking place in reaction k. We follow Woitke & Helling (2003) for the impinging velocity v_j* = min{v_th,j*, 4D/a_p}, where D refers to the molecular diffusivity. Generalizing Equation (9) in OM19, the condensation rate 𝒮_c,k can be expressed as $${\mathcal{S}}_{c,i(k)}=f_{\text{stick }}\frac{x_{v,j^{\ast }}{\rho }_{\mathrm{g}\mathrm{a}\mathrm{s}}{\mu }_{i(k)}}{{\mu }_{j^{\ast }}{\nu }_{k,j^{\ast }}}\pi a_p^2{n}_p{v}_{j^{\ast }}(1-\frac{b_{s,i(k)}}{S_k}),$$(9)

where f_stick is the sticking probability parameter. In the term (1 − b_s,i(k)/S_k), b_s,i(k) denotes the surface fraction covered by condensate i(k), and S_k denotes the saturation ratio of reaction k. The generalization to OM19 lies in the term (1 − b_s,i(k)/S_k), which arises because only the surface covered with condensate i(k) can evaporate, while vapor can condense onto the entire surface of the particle. This implies that different condensates to evenly mix within the particle, which is known as the “dirty ball” assumption (Helling & Woitke 2006; Helling et al. 2008).

The saturation ratio is calculated from the Gibbs energy of reactants and products and partial pressure of all reactants (Berline & Bricker 1969; Woitke et al. 2018), $$S_i={\displaystyle \prod _j}{\left(x_{v,j}\frac{p{m}_{\mathrm{g}\mathrm{a}\mathrm{s}}}{p^{\circ }{\mu }_j{m}_u}\right)}^{v_{k,j}}\exp [-\frac{{\displaystyle \sum _{\text{product }}}{G}^{\circ }(T)-{\displaystyle \sum _{\text{reactant }}}{G}^{\circ }(T)}{RT}],$$(10)

where p is the pressure, G° is the Gibbs formation energy measured at reference pressure p° = 1 bar, and R is the ideal gas constant. We take the NIST-JANAF value for the Gibbs formation energy data (Chase 1998). Again taking Equation (5) as an example, the saturation ratio of MgSiO₃ amounts to $$\begin{array}{ll}{S}_{{\mathrm{M}\mathrm{g}\mathrm{S}\mathrm{i}\mathrm{O}}_3}=& x_{\mathrm{M}\mathrm{g}}{x}_{\mathrm{S}\mathrm{i}\mathrm{O}}{x}_{{\mathrm{H}}_2\mathrm{O}}^2{\left(\frac{p}{p^{\circ }}\right)}^2\frac{{\left(m_{\mathrm{g}\mathrm{a}\mathrm{s}}/m_u\right)}^4}{3.4\times {10}^5}\\ & \times \exp \left[\frac{G_{\mathrm{M}\mathrm{g}}^{\circ }+G_{\mathrm{S}\mathrm{i}\mathrm{O}}^{\circ }+2G_{{\mathrm{H}}_2\mathrm{O}}^{\circ }-G_{{\mathrm{M}\mathrm{g}\mathrm{S}\mathrm{i}\mathrm{O}}_3}^{\circ }}{RT}\right].\\ \hfill =& 5.4\times {10}^{37}\times x_{\mathrm{M}\mathrm{g}}{x}_{\mathrm{S}\mathrm{i}\mathrm{O}}{x}_{{\mathrm{H}}_2\mathrm{O}}^2{\left(\frac{p}{p^{\circ }}\right)}^2(\text{ at }1000\mathrm{K}).\end{array}$$(11)

We note that in the evaluation of the above expression, ${\nu }_{{\mathrm{H}}_2}=-2,x_{{\mathrm{H}}_2}=1$ and $G_{{\mathrm{H}}_2}^{\circ }=0$ by definition. The square dependence of $S_{{\mathrm{M}\mathrm{g}\mathrm{S}\mathrm{i}\mathrm{O}}_3}$ on pressure comes from the summation of stoichiometric number: H₂ on the product side contributes negatively to it.

2.3 Algorithm

The transport Equations (7), (8), and (3) are solved on a grid uniform in log p, assuming hydrostatic balance $\frac{\mathrm{\partial }}{\mathrm{\partial }z}=-\frac{m_{\mathrm{g}\mathrm{a}\mathrm{s}}g}{kT}\frac{\mathrm{\partial }}{\mathrm{\partial }\log P}$. Our code ExoLyn solves the boundary value problem, given by these equations through the nonlinear relaxation method. The total number of equations is n_e = n_v + n_c + 1, where n_v and n_c are the number of vapor and condensate species considered. If the number of grid points is N_grid, the total number of unknowns is n_eN_grid. We follow OM19 for the zero-flux upper boundary condition, $K_{zz}\frac{\mathrm{\partial }{x}_{ci}}{\mathrm{\partial }z}-x_{c,i}{v}_{\text{sed,p }}=0$ for condensate, nuclei and $\frac{\mathrm{\partial }{x}_{v,j}}{\mathrm{\partial }z}=0$ for vapor. In contrast to OM19, the concentration at the lower boundary is taken to be the concentration corresponding to the equilibrium condensation condition (e.g., Woitke et al. 2018). This generalization alleviates the caveat of artificially setting solid concentration to 0, and allows setting the lower boundary condition at a more flexible location, e.g, within or below the cloud deck. The computation domain is automatically adjusted to cover the entire vertical extent of the clouds. The location of the upper boundary is determined at the point where the total cloud mixing ratio is 10⁻¹⁰ times the peak cloud mixing ratio. The lower boundary is set at a pressure level twice the pressure where all clouds evaporate.

Above the cloud base, the particles rain down and raise their concentration at the lower layer of the atmosphere. However, they sublimate below the cloud base, causing a sharp transition from high particle concentration to cloud-free, which potentially destabilizes the solution. Therefore, care must be taken when solving for an atmosphere with strong sedimentation.

The general strategy to obtain the steady state solution is to gradually increase the level of complexity. First, we solve for equilibrium conditions before adding diffusion and settling transport to the equations. Hence we put K_zz = 0 and v_sed = 0 in Equations (7), (8), and (3), which translates to 𝒮_c,k = 0 and 𝒮_n − x_nρ_gas/t_coag = 0. Applying Newton’s method to these equations under element conservation, the mass mixing ratios (x_c,i, x_v,i, and x_n) found are the equilibrium solutions. Next, we gradually add diffusion by multiplying a control parameter f_dif to K_zz. As f_dif gradually increases from 0 to 1, diffusion is added to the transport equation. Finally, the settling term is added to the equations by multiplying v_sed,p with f_set, which similarly increases from 0 to 1. The control parameters, f_dif and f_set, ensure the change between each step is not too large, that is, the solution gradually relaxes to the steady state, and successful convergence is achieved.

Accounting for particle and vapor transport (i.e., given f_dif and f_set), Equations (7), (8), and (3) are solved using a standard relaxation approach. The equations can be written as a general form f(x) = 0, where x = (x_c, x_v, x_n). We start from a guess solution x^old and evaluate the residual of each equation, f(x^old) on each grid point (thus n_e × N_grid in total). Expanding the solution at the guess solution, f(x) ≈ f(x^old) + J(x^old)Δx = 0, where J(x^old) is the n_eN_grid × n_e N_gird Jacobian matrix. Solving this linearized system, we obtain Δx = −J⁻¹(x^old)f(x^old) and x = x^old + Δx is used as the initial guess for next iteration, until the solution converges.

The free parameters in the code are shown in Table 1. We adopt the parameterized T-P profile from Guillot (2010), which takes 7 parameters as shown by the first 7 rows in Table 1. Guillot (2010) follows the visual and thermal flux and calculates the temperature profile under thermal equilibrium. Besides, we note that our code also supports any other T-P profile by interpolating the user-defined T-P grid.

2.4 Radiation transfer

Solid opacities originating from cloud particles are of vital importance in determining transmission or emission spectra. Effective medium theory averages the (complex) refractive index m_i of each component constituting the “dirty cloud particle”, given their volume mixing ratio f_V. Bruggeman theory (see, for example, Choy 2016) assuming randomly mixed spherical monomer made of homogeneous material, gives the expression from which the averaged refractive index m can be solved: $${\displaystyle \sum _i}{f}_{V,i}\frac{m_i^2-m^2}{m_i^2+2m^2}=0,$$(12)

where the refractive indices of each pure material, m_i, has been compiled by Kitzmann & Heng (2018) (Table 2). The Bruggeman mixing rule is widely adopted in the literature (Mishchenko et al. 2016; Lines et al. 2018; Min et al. 2020), but does rely on assumptions which may in certain situations be less appropriate. In particular, while we use the distribution of hollow spheres (DHS) model for the shape of the particles (next paragraph), the mixing rule in Eq. (12) assumes that each component within the particle is made up of tiny spherical constituents, homogeneously mixed throughout the particle. The spherical nature is sometimes problematic as the effective dielectric constant may be different from other shapes (e.g., Fe needles are far more opaque than Fe spherules). Recently, Kiefer et al. (2024b) have compared how different opacity mixing treatments affect the spectra, including different mixing rules and core-shell particle structure. It still remains an open question what the morphology of a heterogeneous cloud particle is in reality and which opacity treatment to take.

Besides the material refractive index, cloud particles’ shape and size also determine their absorption and scattering properties. We assume the particle shape follows DHS (Min et al. 2005), because, on average, the optical properties of DHS can mimic an ensemble of irregularly shaped particles. Min et al. (2005) and Mollière et al. (2017) showed that the scattering properties of DHS is different from spherical particles whose scattering properties are given by Mie theory at mid-IR wavelength. Based on DHS, we use software OpTool (Dominik et al. 2021) to calculate the particle scattering and absorption opacity.

Apart from clouds, gas molecules also contribute to the atmosphere opacity, with their spectral lines indicative of the atmospheric component. The gas opacity in our model includes all the vapor species involved in cloud-formation (Table 3), plus other gas molecules abundant in exoplanet atmosphere – CO, CO₂, and CH₄. We also consider Rayleigh scattering of H₂, He, and collision-induced absorption (H₂-H₂ and H₂-He). After cloud formation we force the gas phase concentration to local thermo-chemical equilibrium values. This is done by mixing the “leftover vapor” and background gas (H₂, He, CO, CO₂, and CH₄) locally and running FastChem (Stock et al. 2022) to compute gas composition assuming locally chemical equilibrium. How disequilibrium chemistry such as photoevaporation and diffusive quenching (Madhusudhan & Seager 2011; Tsai et al. 2017; Kiefer et al. 2024a) affects the spectrum, is beyond the scope of this paper. Finally petitRADTRANS (Mollière et al. 2019) is run to calculate the transmission and emission spectra, based on the equilibrium gas concentration and the cloud opacity described in the above paragraph. The reference pressure in petitRADTRANS is set to 10⁻⁴ bar for all cases.

3 Results

Similar to OM19, our default model investigates a generic hot Jupiter at a distance of 0.05 au from a solar-like star. The temperature is prescribed by Guillot (2010) accordingly (see Figure 2), with internal heat flux corresponding to T_int = 500 K. Taking a diffusivity of K_zz = 10⁸ cm²s⁻¹ and solar elemental abundance, we run ExoLyn to investigate the cloud forming on the above hot Jupiter. The envelope bulk elemental abundances at the lower boundary is listed in the first column of Table 4, and the reactions we consider in the default model are listed in Table 3. In Sections 3.1 and 3.2, we present the results of the default model, focusing on solid and gas, respectively. In Section 3.3, we compare ExoLyn with other cloud formation models. Then, in Section 3.4, we discuss how K_zz, the nucleation profile, the choice of cloud species, and the molecular abundance at the boundary condition influences the cloud formation. Following that, in Section 3.5, radiation transfer calculations are performed to investigate how the gas and cloud species influence the transmission spectrum.

Fig. 2 Default temperature-pressure profile of the hot Jupiter (Section 3), self-luminous (Section 4.1), and sub-Neptune planet (Section 4.2), following Guillot (2010). The irradiation temperatures of the hot Jupiter, self-luminous planet, and the sub-Neptune planet are taken to be 500 K, 1000 K, and 300 K, correspondingly. Condensation curves of Mg2SiO4, SiO2, Fe, Al2O3 (for hot Jupiter and self-luminous planet), NaCl and ZnS (for sub-Neptune) are shown as colored lines. The solid, dashed, and dotted saturation curves denote solar, 10 times solar, and 100 times solar metallicity, respectively.

3.1 Cloud structure

Figure 3 presents the solid, gas, nuclei mixing ratios, and particle size in the default cloud model. The total column density of the cloud amounts to 0.03 g cm⁻². The cloud shows a twolayer structure: at around 10⁻² bar, it is dominated by MgSiO₃, Mg₂SiO₄, and SiO₂ (refered to as magnesiumsilicate or Mg-Si hereafter), whereas beyond 2 × 10⁻² bar, Fe-dominant condensates take over. The layered structure is consistent with the condensation sequence of magnesiumsilicate (≈1460 K) and iron condensate (≈1600 K) species, respectively. As a result, the magnesiumsilicate condensate group evaporates higher in the atmosphere than Fe. There is no cloud below 1 bar, because even the most refractory species considered here, Al₂O₃, evaporates. The right panel of Figure 3 gives a graphical illustration of the components making up the cloud particles by mass. At the location of the magnesiumsilicate cloud (P < 0.02 bar), though the temperature is below the iron condensation point, the cloud particles are depleted of iron, because the iron element settles to deeper atmosphere. In the iron-dominated regions of the cloud, at around 0.2 bar, we still find 1% of MgO and SiO₂, even though equilibrium chemistry would predict their absence due to the temperature higher than the magnesiumsilicates evaporation value. This result arises from the rapid downward transport by settling and the “dirty ball” assumption: the magnesiumsilicates are transported to the level due to settling and they do not sublimate efficiently limited by the surface coverage. Due to settling, the condensate mass mixing ratio increases with depth and peaks near the bottom of the cloud (≈10⁻² bar for magnesiumsilicates and ≈0.2 bar for Fe condensates). During settling, coagulation, and condensation grow the particle larger, reaching 4 μm at 0.5 bar (Figure 3) and leading to faster settling at lower levels. As the temperature increases with pressure, the cloud quickly evaporates at the bottom, causing a steeper slope at the cloud bottom than the cloud top.

Nuclei mass mixing ratio x_n, which is proportional to particle number density, is subject to both production and coagulation. The nuclei mixing ratio peaks at a height corresponding to 6 × 10⁻⁵ bar (see Figure 3), where they are injected in the default simulation. From their formation location the particles diffuse up and settle down, leading to decreasing particle number density, toward both upper and lower atmosphere. Particles hardly grow beyond their initial size (1 nm) at the location where the nuclei are injected, because they diffuse away before coagulate. As the nuclei mixing ratio decreases toward higher pressure, the particle size generally increases, despite a peak at 10⁻² bar, which is due to the magnesiumsilicate cloud.

Fig. 3 Atmosphere structure of the default hot Jupiter model. The thick colored line shows the cloud mass mixing ratio summed over all species. The dominant species by mass is illustrated by the color of the thick line. Vapor mixing ratios are shown by dashed lines. The dotted black line denotes the nuclei mass mixing ratio and the purple thin line denotes the cloud particle size, which should be referred to the upper axis. In the figure we also label the evaporation temperature of magnesiumsilicates (1460 K) and iron (1600 K) with horizontal dashed lines. The gray band refers to where the cloud optical depth in the IR band (1 μm − 20 μm), integrated from the top of the atmosphere, reaches 1. The right panel shows the composition of cloud particles in mass.

3.2 Molecular gas profile

The vapor mass mixing ratios are shown by dashed lines in Figure 3. Near the cloud base, all vapor species approach their boundary abundances, due to the evaporation of all cloud species. Fe, TiO, Al, SiO, and Mg are significantly depleted in the upper layers, because they are the limiting species that form Fe(s), TiO₂, Al₂O₃, and magnesiumsilicates, respectively. The cloud formation (condensation) results in a minimum region in the vapor cencentration (e.g., 4 × 10⁻⁴ bar for H₂O), to which vapor from lower regions diffuses. Above p = 4 × 10⁻⁴ bar, vapor mixing ratio rises again with atmospheric height. The vapor here comes from the evaporation of the cloud particles, due to lower pressure. The vapor mass mixing ratio approaches a constant value at the upper boundary, because the diffusive flux is imposed to be 0. In addition, non-limiting vapor species also deplete, albeit less prominently. For example, H₂O vapor mass mixing ratio decreases by a factor of 2 in the upper atmosphere. Therefore, consumption of condensable vapor due to cloud formation may result in less distinct molecular lines than the cloud-free case. This quenching effect implies that assuming constant gas abundance during retrieval may underestimate the metallicity of an exoplanet. The C/O ratio of the atmosphere is also affected by the formation of cloud. As the cloud formation consumes O-bearing vapor such as H₂O, consequently the atmospheric C/O ratio increases from 0.52 to 0.54

3.3 Comparison with previous studies

In this section we compare ExoLyn with Ackerman & Marley (2001) (AM01) and OM19. The main advantage of ExoLyn lies in the self-consistent formation of multiple cloud species, compared to the single species cloud in AM01 and OM19. Min et al. (2020) extended the OM19 method, but the multispecies cloud formation sequence is prescribed. For a meaningful comparison, we only include MgSiO₃ condensate in three models. Assuming equilibrium between sedimentation and diffusion, AM01 solves the equation for x_t, the total molecular concentration in gas and solid phase, $$K_{zz}\frac{\mathrm{\partial }{x}_t}{\mathrm{\partial }z}+f_{\mathrm{s}\mathrm{e}\mathrm{d}}{w}_{\ast }{x}_c=0,$$(13)

where w_* is the convective velocity scale and f_sed is an order of unity settling coefficient. We take default f_sed = 1, w_* = K_zz/H, where H is the scale height of the atmosphere, and $$x_c=max\{0,x_t-x_s\},$$(14)

where x_s is the saturation mixing ratio of MgSiO₃. The MgSiO₃ saturation vapor pressure in both AM01 and OM19 is P_s = 1.04 × 10¹⁷ exp[58 663/T] dyn m⁻² (Ackerman & Marley 2001). After x_t is calculated, cloud condensate mixing ratio x_c can be deduced through Equation (14).

The results are shown in Figure 4. Generally a good match is seen between the three models. For a single species cloud, the transport equations we solve are the same as those in OM19, except the expression for saturation ratio and condensation rate (Figure 1). As the saturation pressure used in AM01 and OM19 can translate to a reasonable estimation to the saturation ratio calculated from a more formal way (Equation (11)), OM19 predicts a similar cloud profile to ours. However, we warn that unlike OM19’s assumption, the precise saturation ratio depends not only on temperature, but also on the molecular abundance (Equation (11)), such as metallicity or C/O ratio. Therefore, the saturation pressure of MgSiO₃ adopted in OM19 cannot be universally applied.

The settling velocity in AM01 model is characterized by the f_sed parameter, rather than calculated from physical principles. Though the default AM01 run assuming f_sed = 1 agrees with our model reasonably well, it generally underestimates the cloud mixing ratio by a factor of 2. Furthermore, if f_sed is raised to 2 (thin yellow line), the extent of the cloud in AM01 shrinks significantly and no longer matches our model. Therefore, the overall agreement between AM01 and our model seems to be coincidental in the sense that the settling velocity corresponding to f_sed = 1 happens to be similar to the settling velocity of a ~ 10⁻² μm particle. This can be tested through changing K_zz and other parameters. In particular, the AM01 cloud is independent of both K_zz, which cancels in Equation (13) through the choice of w* and the nucleation profile. In the right panel of Figure 4, f_sed = 2 provides a better estimation of the K_zz = 10⁷ cm²s⁻¹ cloud than f_sed = 1. However, as we subsequently see in Section 3.4, K_zz greatly affects the cloud thickness and concentration.

Fig. 4 Comparison of the full ExoLyn model (default), ExoLyn restricted with MgSiO3 condensates only, the Ormel & Min (2019) model, and Ackerman & Marley (2001). The thickness of the lines and the symbols on the lines indicates the particle size, except in the Ackerman & Marley (2001) runs, where the particle size is not solved for as it decouples with the cloud concentration in the model. Right: the same as left, but Kzz decreased to 107 cm2s−1. In this case, the fsed = 2 run (thin line) provides a better correspondence to other models.

3.4 Parameter study

3.4.1 Atmosphere diffusivity

The eddy diffusivity K_zz determines the extent of cloud, because the vertical turbulent motion of gas serves as an effective diffusivity term, blowing up the trace gas and solid. Formally K_zz is derived from mixing length theory (e.g., Wang et al. 2015) or from hydrodynamical simulations by measuring the mixing flux and calculate how it correlates with the concentration gradient. (e.g., Zhang & Showman 2018). K_zz depends on the atmospheric thermal structure, heating flux, precipitation, and phase change, as is suggested by Ge et al. (2024); Leconte et al. (2024). Due to the uncertainty in K_zz, in our study we use a constant K_zz as a representative value over the entire atmosphere.

We test how the choice of K_zz influences the clouds in Figure 5. The effect of increasing K_zz is twofold: (i) cloud particles are transported upward – counteracting settling; (ii) the cloud-forming vapor is replenished more efficiently. As a result, the total cloud mixing ratio increases with K_zz. A boxlike cloud, extending from 10⁻⁶ bar to 1 bar, appears for the K_zz = 10¹⁰ cm²s⁻¹ run. In contrast, the K_zz = 10⁴ cm²s⁻¹ or K_zz = 10⁶ cm²s⁻¹ runs exhibit a clear two-layered cloud structure, because settling is more effective. As a result of the thinner cloud layer, the cloud optical depth τ_cloud = 1 surface lies deeper in the atmosphere for lower K_zz. In the case of very low diffusivity, K_zz = 10⁴ cm²s⁻¹, the cloud opacity never reaches unity.

Comparing the K_zz = 10⁶ cm²s⁻¹ to K_zz = 10⁴ cm²s⁻¹ case, the particle number density roughly overlaps, because coagulation dominates the particle number density distribution when the transport is weak. As a result, the particle size increases mainly due to the increasing amount of condensates. When K_zz increases to 10⁸ cm²s⁻¹, more particles are transported downward up to 10⁻² bar due to stronger diffusion. The particles grow so large that they settle quickly near the cloud bottom (<10⁻² bar). Therefore, the particle nuclei near the cloud bottom deplete and in turn raise the particle sizes to 4 μm. At the other extreme, for K_zz > 10¹⁰ cm²s⁻¹ transport is dominated by diffusion, and the particle number density at the cloud lower boundary increases with stronger diffusive supply. Consequently, the particle size at the cloud bottom decreases when the diffusion is strong. For K_zz = 10¹⁰ cm²s⁻¹, the particle size reaches 0.1 μm throughout the cloud.

We also plot the Fe mixing ratio to examine how vapor content is affected by K_zz. With stronger diffusion, there is more Fe vapor diffused up, although its concentration above the cloud formation layer is negligible in all K_zz runs. Our result implies that in a strongly turbulent atmosphere, the molecular features become less apparent, or even indiscernible, given the thicker clouds.

Fig. 5 Cloud profiles under different eddy diffusivity parameters Kzz, which serve as effective vertical turbulence. Left: total cloud mixing ratio (solid lines) and Fe vapor mixing ratio (dashed lines) as an example. Right: nuclei mixing ratio (solid lines) and particle size (dashed lines). Similar to Figure 3, the location where the cloud optical depth reaches unity is shown in the left panel by the horizontal bands with corresponding colors. We note that there is no blue band corresponding to Kzz = 104 cm2s−1, because the cloud optical depth never reaches 1.

3.4.2 Nucleation

As discussed in Section 2.1, how much and where nuclei form on exoplanets is a controversial topic of atmosphere modeling. In ExoLyn, the profile of the nucleation production rate is entirely parameterized. We find that the nucleation profile plays a less important role in affecting the cloud than eddy diffusivity. In Figure 6 we vary the imposed nucleation profile, including nucleation strength ${\dot{\mathrm{\Sigma }}}_n$, nucleation position P_n, nucleation spread σ_n, and nuclei particle size a_p0. Increasing the nucleation rate by 4 orders of magnitude has little effect on the cloud structure. This is because the amount of cloud condensation primarily depends on the vapor concentration, rather than the amount of nuclei. When there are ample amount of nuclei, coagulation among particles mitigates the effect of further increasing the nuclei injection rate. On the other hand, when fewer nuclei are injected, the cloud particles’ sizes grow larger, resulting in a more compact cloud that lies deeper in the atmosphere.

Changing the nucleation location hardly influence the clouds below it (light and dark green line in Figure 6). This is because nuclei form higher in the atmosphere can settle down to lower levels. For the same reason, the cloud profile changes little when there is larger spread in the nucleation (the pink line). However, when nuclei are injected deeper in the atmosphere, the clouds will not extend much upward because the particles have grown too large to reach higher regions by diffusing up against gravity, leading to less nuclei concentration above the nucleation pressure. We do observe an increase of particle size at around 10⁻³ bar and a reduction of the cloud extent when nucleation proceeds at a pressure scale which is two orders of magnitude higher than in the default model. In hot Jupiters, (TiO₂)_N clustering can be a major source of nucleation. Deriving updated Gibbs energy of (TiO₂)_N clusters, Lee et al. (2015) found TiO₂ nucleation happens around a few 10⁻⁴ bar, broadly consistent with our fiducial choice of the nucleation profile. Photochemistry haze, another potential source of nucleation on hot Jupiter’s atmosphere, is found to form around 10⁻⁵ to 10⁻⁶ bar (Arfaux & Lavvas 2022; Lavvas & Koskinen 2017; Arfaux & Lavvas 2024). Therefore, though the actual location of nucleation on a hot Jupiter is uncertain, our fiducial clouds result can represent the real case.

Lastly, if the nuclei initial size is increased to 10 nm, the upper cloud extent is diminished. The effect is similar to decreasing ${\dot{\mathrm{\Sigma }}}_n$, which also acts to increase the cloud particles’sizes. A larger nuclei size is more appropriate if the Kelvin effect is strong. The Kelvin effect expresses the suppression of condensation due to the curvature of the nuclei. When the particle size (a) is small and the surface tension of the material (Γ) is high, it reduces the condensation exponentially by a factor $$\exp \left(\frac{2\mathrm{\Gamma }{\rho }_{•}\mu m_u}{k_{b}Ta}\right)$$(15)

where μm_u is the molecular weight and T is the temperature. To overcome the Kelvin effect, nuclei are likely to coagulate to a larger radius, an effect that can be modeled in ExoLyn by choosing a larger size for the nuclei. In summary, although ExoLyn cloud results are insensitive to nucleation, assumptions for large nucleation rates, its nucleation parameters provide the flexibility to mimic the real nucleation profile.

Fig. 6 Sensitivity of the nucleation profiles on the cloud profiles. In the figure we changed the nucleation strength ${\dot{\mathrm{\Sigma }}}_n$, nucleation location Pn, spread of nucleation σn, and nuclei size ap0. The parameters are changed with respect to the default value in Table 1. The vertical lines indicate the location of nucleation, with its length meaning 1σ dispersion of the nucleation profile.

Fig. 7 Cloud profiles with different condensation reactions. The full set of reactions Table 3 is shown as the black line. The green and yellow line only includes the magnesiumsilicate group versus the iron and refractory reaction group, that is, the first four rows and the last six rows in Table 3. The vertical lines represent where the near-IR optical depth accumulates to unity.

3.4.3 Condensation reactions

Figure 7 illustrates how the included condensation reactions affect the physical properties of the cloud. For “only MgSi”, we only include the first 4 reactions in Table 3, whereas for “only Fe” we only include the last 6 reactions, which is representative of the more refractory cloud species. With “only MgSi” reactions, one can recover the entire cloud mixing ratio above the magnesium silicate evaporation line, while with “only Fe” model one can recover the iron and refractory peaks at 10⁻² to 1 bar. Compared to the full model (Figure 3), in the “only Fe” model, the iron and refractory cloud also occupies the height dominated by magnesiumsilicate cloud in the default model. Though Fe is more opaque than MgSi, the τ = 1 of “only Fe” model is similar to that of the default model, because the Fe particles are smaller. As the τ = 1 surface (given by the vertical lines on the right) lies roughly in the magnesium silicate cloud, only including magnesium silicates condensates is still an acceptable estimation to the full reactions. However, this may not hold when the cloud optical depth would be lower, for instance, when K_zz ≲ 10⁶ cm²s⁻¹.

3.4.4 The effect of envelope bulk abundance on the cloud properties

For the default model, we assume solar values (Asplund et al. 2009) for the planet’s envelope elemental abundances. Based on that, we assume that CO consumes all C and H₂O takes the leftover O atoms. This gives $x_{{\mathrm{H}}_2\mathrm{O}}=1.4\times {10}^{-3}$ and x_CO = 3.2 × 10⁻³ at the lower boundary, which is consistent with equilibrium chemistry calculated with FastChem (Stock et al. 2022). In Table 4 and Figure 8 we consider two cases with different elemental abundances: increasing the metallicity to 10 times solar values (10×-[Z]); and increasing C/O to n_C/n_O = 0.9.

In the 10×-[Z] case, the cloud bottom extends deeper in the atmosphere up to 3 bar. This is because the condensation curves (dashed lines in Figure 2) shifts rightward when increasing vapor concentrations, intersecting the T-P profile at higher temperatures and pressures. For the same reason, the Mg₂SiO₄ cloud does not evaporate until 0.3 bar. The formation of the lowhanging Mg₂SiO₄ cloud consumes almost all Mg, causing low Mg concentration above 0.03 bar, (see Figure 8a). Therefore a cloud layer dominated by SiO₂ forms on top of this, with its maximum mixing ratio reaching 10⁻³. The cloud particles sizes in the 10×-[Z] run become over an order of magnitude larger than in the default model (Figure 3) with sizes in Fe-dominated region reaching ≈70 μm. The large size leads to fast settling, resulting in a thinner iron layer and a razer-thin refractory cloud deeper in the atmosphere (around 2 bar).

Though all molecular mixing ratios at the lower boundary increase by a factor of 10, we do not see a concomitant growth of the cloud in terms of column density. In fact only 3 times more cloud (in terms of column density) is formed compared to the default case. The reason is that settling provides a negative feedback to cloud growth: the more condensates formed force the cloud to settle faster to hotter, lower region, giving negative feedback to the cloud formation. Interestingly, above 2×10⁻² bar the Mg-Si cloud mixing ratio and particles sizes remain unchanged compared to the default solar metallicity run. Although the planet bulk vapor abundances are enhanced, the excess vapor is largely consumed by the “deep” dense cloud and is unable to reach the extended upper cloud (P < 2 × 10⁻² bar). The decoupling of metallicity and cloud concentration (especially the upper cloud concentration) suggests that even with the presence of clouds, a metal-rich planet may yet exhibit strong molecular lines, which we subsequently demonstrate further in Section 3.5.

We also run a model in which the C/O ratio is increased such that n_C/n_O = 0.9, but where the oxygen plus carbon elemental abundance by mass, x_C + x_O remains at solar values (Figure 8b). As a result, CO consumes most of the oxygen, decreasing the H₂O abundance at the lower boundary to $x_{{\mathrm{H}}_2\mathrm{O}}=6.8\times {10}^{-5}$, that is, less water is present in the atmosphere. The scarcity of H₂O suppresses the formation of Mg-Si, Al₂O₃, and TiO₂ cloud particles, which all require water vapor (Table 3). As a result, Fe becomes the dominant cloud species (Figure 8b). A similar trend has been reported by Helling et al. (2017) It is of course well-known that C/O=1 constitutes a clear boundary above which the chemistry of the atmosphere and composition of the cloud switches from O-rich to C-rich. But here we find that even before C/O=1 is reached the cloud composition shows a clear dependence with increasing envelope bulk C/O, becoming Fe-rich.

As the surface tension of Fe (1870 dyn cm⁻¹, Brillo & Egry 2005) is higher than silicates (307 dyn cm⁻¹ for SiO₂, Janz 1967), Fe vapor can only condense onto larger nuclei due to the Kelvin effect (Eq. (15)). Therefore, the nuclei particles need to coagulate before Fe can condense onto them. To mimic a more realistic nucleation profile when the particle is made primarily of Fe, we test the C/O=0.9 case where larger nuclei particles are injected into the deeper atmosphere. We take the nuclei particle size a_p0 = 5 nm, a factor of 5 higher than the default parameter and nucleation pressure P_n = 6 × 10⁻⁴ bar, an order of magnitude higher. The result is shown in Figure 8c. As explored in Section 3.4.2, the nucleation profile has minimal effect on the cloud formation, except that the clouds hang lower in the atmosphere. Larger particle size makes Fe particles harder to diffuse up to the Mg-Si condensing region. Therefore, the upper cloud is composed of comparatively more Mg-Si, though Fe is still the dominant species, unaffected by the choice of nucleation height and size.

Fig. 8 Effect of vapor abundances at the lower boundary of atmosphere on the solid, vapor, nuclei profile, particle size, and cloud composition. (a) 10×-Z hot Jupiter model; (b) C/O=0.9 hot Jupiter model; (c) C/O=0.9 hot Jupiter model, similar to (b), but with injection of larger nuclei deeper in the atmosphere, mimicking nucleation featuring high particle surface tension. The thick solid, dashed, and dotted lines indicate total cloud, vapor, and nuclei mixing ratios, respectively. The purple solid line shows the particle size. The right subpanel in each panel shows the composition of the cloud particles in mass. Compared to the default model (Figure 3) the upper cloud layers of the 10×-[Z] run become devoid in Mg (SiO2 dominates) while in the C/O=0.9 run Fe-bearing cloud particles make it into the upper regions.

3.5 Synthetic spectrum

In Figure 9 we plot the transmission spectra of the default hot Jupiter model (black curve), focusing on 0.5 μm to 20 μm, which covers the wavelength of HST, Spitzer, and JWST. To demonstrate the contribution of clouds, we also plot in gray the spectrum that would emerge when the solid particles are removed, while preserving the gas species. Compared with the clear model, the presence of clouds enhances the transit depth on the order of ~10² ppm. At the short wavelength end, the cloudy spectrum shows a slope consistent with Rayleigh scattering, reflecting the presence of sub-μm size particles. Furthermore, a distinct silicate feature shows up at 8 to 10 μm, unambiguously indicating the formation of silicate particles. Though Fe, Al₂O₃, and TiO₂ cloud layers also form (Figure 3), they lie under the τ = 1 surface and are thus obscured by the magnesiumsilicate cloud layer.

Without clouds the amplitude of the molecular lines can reach 100 to 200 ppm. The thick clouds in the atmosphere reduce the amplitude to ~50 ppm. From the spectrum we can identify the CO feature at 4.6 μm and the H₂O features at 2.5 and 5.6 μm, though cloud formation depletes the H₂O mixing ratio from its boundary value (1.4 × 10⁻³) to $x_{{\mathrm{H}}_2\mathrm{O}}\approx 8\times {10}^{-4}$ (Figure 3). Besides the 2.9 μm H₂O line, H₂S can still be marginally identified as a shoulder to the right, at 3.5 μm. On the other hand, the features of CH₄ are not found in the spectrum, because most of the carbon exists in the form of CO in the chemical equilibrium state for the temperatures considered here.

We further show in Figure 9 how K_zz affects the spectrum. As we have discussed in Section 3.4, the higher the K_zz, the thicker the cloud, suppressing the molecular features. Indeed we see that the scattering slope dominates the entire spectrum below 5 μm when K_zz = 10¹⁰ cm² s⁻¹. The 10 μm silicate bump also becomes more prominent. With decreasing K_zz, the water lines in the nearIR band become more prominent, together with the 3.5 μm H₂S feature because the clouds become thinner. Decreasing K_zz further to the extremely low 10⁴ cm² s⁻¹ even the 10 μm silicate feature vanishes. The silicate cloud lies so deep in the atmosphere that the gas optical depth from 2 to 20 μm above it already exceeds 1. However, clouds still make a difference in the transit depth at wavelengths smaller than 2 μm, because the cloud opacity still contributes to the total opacity, raising the height where the total optical depth reaches unity.

We show how the boundary elemental abundance affects the spectrum in Figure 10. When the metallicity is increased by a factor of 10, the H₂O, H₂S, and CO₂ concentrations are boosted, but not much more silicate cloud particles form in the upper atmosphere (Section 3.4.4). As a result, these molecular features become more prominent (see Fig. 10). The most prominent example is H₂O features at 1.4, 1.9, and 2.8 μm. As the cloud density and particle sizes in the optically thin region (above 10⁻² bar) closely resemble the default run (see Section 3.4.4), the 10 μm feature strength remain unchanged. Another remarkable result of the 10×-[Z] run is the formation of a SiO₂ cloud in the upper atmosphere (Figure 8a, also Woitke et al. 2020). As a result, a 12 μm solid feature arises due to SiO₂ (s).

For the C/O=0.9 case, there is little H₂O that can supply the formation of MgSi clouds (see Section 3.4.4). Therefore the MgSiO₃ and SiO₂ features around 10 μm disappear (Figure 10, red curve). The opacity increases due to higher Fe mixing ratio in the upper atmosphere. Compared to the default run, the reduced amount of Mg-Si solid leads to weaker settling, hence retaining more Fe in the upper atmosphere. In addition, because the Fe-to-MgSi ratio is higher than the default case, the effective refractive index of cloud particles is closer to that of Fe, showing strong absorption at near-IR wavelengths. As a result the τ = 1 surface is raised to 6 × 10⁻³ bar at wavelengths smaller than 5 μm where most molecular features lie. However, the opacity in the default model may be somewhat suppressed due to the assumptions inherent to the Bruggeman mixing rule with spherical inclusions (Section 2.4). Featureless near-IR transmission spectra of hot Jupiter planets could be attributed to the presence of Fe cloud particles or high turbulent diffusivity K_zz, but mid-IR spectra at 10 μm wavelength can help distinguish the two scenarios. When injecting nuclei deeper in the atmosphere, the near-IR part of the spectrum shows less transit depth, but the shape and amplitude of the spectral features is similar to the C/O=0.9 case (red curve). At 10 μm, the Mg-Si features becomes more apparent. The reason for these changes are that the cloud is located deeper in the atmosphere and that its particles are composed of more transparent Mg-Si material.

Fig. 9 Transmission spectrum of a cloudy hot Jupiter (black line). The gray line shows the transmission spectrum of a clear atmosphere, leaving the gas components unchanged. Kzz is varied to show how cloud formation changes the atmospheric transmission spectrum. The horizontal lines and arrows indicate molecular features. At long wavelengths the Kzz = 104 cm2 s−1 run curves are indistinguishable from the “clear” curve.

Fig. 10 Transmission spectrum of models varying metallicity and C/O ratio. The dark and light line colors refer to cloudy and clear atmosphere, respectively. The black and gray lines are the default cloudy and clear model (Section 3.1). 10×-[Z] and C/O=0.9 cases are shown by the blue lines and the red lines. For the C/O=0.9 run, we additionally conduct a run injecting larger (5 nm) nuclei particles deeper in the atmosphere to accounting for the Kelvin effect (brown curve).

Fig. 11 Cloud structure and emission spectrum of a self-luminous planet with physical properties similar to HR 8799 e. The left panel shows the cloud, gas, and nuclei mixing ratios, with every components the same meaning as Figure 8. The brown thick and thin lines present the temperature profiles we used in the clear and cloudy models, respectively. The right panel shows the synthetic emission spectrum of a cloudy atmosphere (the black line) and a clear atmosphere (the blue line), whose effective temperatures are both 1000 K. For reference, the spectrum of a 1000 K black body radiation is also plotted as the dashed line.

4 Application to self-luminous planets and sub-Neptunes

ExoLyn is designed to be broadly applicable. In this section, we demonstrate the adaptability of ExoLyn by applying it to a selfluminous planet and a sub-Neptune planet, with the aim to show how clouds affect the emission and transmission spectrum.

4.1 Self-luminous planet

Self-luminous planets are planets with a temperature structure entirely governed by their internal thermal flux. Detected through direct imaging, these planets are relatively massive and orbit far from their host stars (P ⪆ 10 yr). These planets are typically young (⪅100 Myr) and are still in the process of radiating away their formation energy. Self-luminous planets therefore have a hot interior and much cooler atmosphere. Consequently, we expect clouds to form. When the effective temperature is fixed, the presence of clouds would make the spectrum redder (Currie et al. 2023). For example, the strong emission of 51 Eri b in Lp (3.8 μm) band, probed by Gemini, lends support to a thick cloud with small particle size (Macintosh et al. 2015; Samland et al. 2017).

To explore how clouds influence the emission spectrum of self-luminous planets, we apply ExoLyn to a 10 M_J distant giant planet akin to HR 8799 e (GRAVITY Collaboration 2019). We take the Guillot (2010) temperature profile for a self-luminous planet, where the outgoing irradiation flux determines the temperature profile. The controlling parameters here is the internal irradiation temperature T_int and the constant, gray opacity at the irradiation IR band κ_IR. For the former we take the internal irradiation flux corresponding to a temperature of 1000 K (Bonnefoy et al. 2016; GRAVITY Collaboration 2019). If a self-luminous planet atmosphere, cloudy or not, is in thermal equilibrium, it will radiate away the same amount of energy it receives from the interior. To isolate the role clouds play in a self-luminous planet, we tune the IR opacity parameter in the Guillot model, such that the effective temperature T_eff at which the planet (cloudy or clear) radiates is the same as T_int = 1000 K. For a clear atmosphere, T_int = T_eff translates to κ_IR = 0.04 cm² g⁻¹ and for a cloudy atmosphere this gives κ_IR = 0.1 cm² g⁻¹, consistent with the addition of the cloud opacity. The temperature profile of the cloudy run is shown in Figure 2. As the lower atmosphere can be as hot as 2000 K, the condensate species and reactions are the same as listed in Table 3. Here, the goal is not to fit the observed spectrum of HR 8799 e in particular, but rather to demonstrate how ExoLyn can be applied toward self-luminous planets and how clouds change the planet’s emission spectrum in general.

The cloud structure and emission spectrum calculated from petitRADTRANS are shown in Figure 11. Similar to the hot Jupiter, the self-luminous planet clouds also show a layered structure. The upper cloud above 1.5 bar is dominated by magnesium silicate. Due to the lack of stellar irradiation heating, most of the region where clouds form remains cooler than 1000 K. Therefore magnesiumsilicate particles do not evaporate until 1.5 bar (1600 K). With increasing pressure, the temperature quickly increases with depth, reaching the iron evaporation temperature at 3.5 bar. Consequently, iron becomes the dominant cloud species in the 1.5 bar to 3.5 bar range. But we also note that due to diffusion, Fe is efficiently transported to the magnesiumsilicate cloud region, comprising 15% of the cloud particle mass there.

In the right panel we show the synthetic emission spectrum. The spectrum of a clear atmosphere devoid of clouds and a cloudy atmosphere are drawn with the blue and black lines, respectively. The emission spectrum of the clear self-luminous planet resembles a black-body with additional molecular absorption features. The absorption lines arise because the photons from the hotter interior get absorbed by molecules in the cooler upper region.

In the cloudy run, photons originating from the hot interior are absorbed by the cloud particles in the atmosphere and re-emitted at a cooler emission temperature. Therefore, the spectrum is redder compared to the clear one, even though the effective temperature of the two atmospheres are the same. Moreover, the molecular lines <2 μm are suppressed due to cloud absorption. We observe absorption lines of CO and H₂O only in the near-IR band. For example, the dip from 3 to 4 μm is due to absorption of CO. The color of the spectrum and the line strength at near-IR wavelengths allow us to distinguish a cloud-free self-luminous planet from a cloudy one. Though fitting a real spectrum is beyond the scope of this work, we note that SPHERE (Zurlo et al. 2016), GPI (Greenbaum et al. 2018), and VLTI/GRAVITY (GRAVITY Collaboration 2019; Mollière et al. 2020; Nasedkin et al. 2024) spectrum of HR8799 e lies in between our cloudy and clear cases, suggesting less thick cloud as our cloudy run, that is, smaller K_zz.

Fig. 12 Cloud structure and transmission spectrum of a sub-Neptune planet with physical properties similar to GJ 1214 b. The condensate species considered here is NaCl, KCl, ZnS, and Na2S. In the left panel, the meaning of each line is the same as Figure 8. The cloud is dominated by Na2S (upper part) and NaCl (cloud bottom). The right panel shows the synthetic spectrum and the blue line shows one set of parameters that leads to flat spectrum in 1 to 2 μm range, consistent with observation (Kreidberg et al. 2014).

4.2 Sub-Neptunes

In this section we focus on a warm sub-Neptune planet, whose temperature and mass is lower than that of a hot Jupiter. As the condensation temperatures of magnesiumsilicates are always far higher than the local atmospheric pressure (see Figure 2), the condensates in Table 3 are likely to form deeper in the interior or near the surface of the planet. The condensates forming in the atmosphere regions of sub-Neptune planets are usually salts, such as KCl and NaCl (Morley et al. 2013; Ohno & Okuzumi 2018).

Here we consider a sub-Neptune with physical properties similar to GJ 1214 b, a 8.17 M_⊕ and 2.74 R_⊕ planet 0.015 au away from a 0.216 R_⊙ star (Carter et al. 2011; Cloutier et al. 2021; Mahajan et al. 2024), with irradiation temperature at 300 K. JWST observations of GJ 1214 b reveal thick and reflective cloud layers (Kempton et al. 2023). Thus the condensate species chosen are NaCl, KCl, Na₂S, and ZnS, whose condensation temperatures lie in the atmospheric temperature range (Morley et al. 2013; Charnay et al. 2015; Gao & Benneke 2018). We take a 100 times solar metallicity (Table 4), with Na, K, Cl, Zn, and S existing in the form of Na, K, HCl, Zn, and H₂S vapor. The chemical reactions that form these species are listed in the lower part of Table 3. Other parameters, such as K_zz and the nucleation description, are kept the same as Table 1, except the mean molecular weight m_gas, for which 5m_u is taken, consistent with its metallicity.

The atmosphere structure and the resulting transmission spectrum for the sub-Neptune planet are shown in Figure 12. As a result of the solar element ratios, HCl is the limiting species for NaCl and KCl, and Na is limiting the formation of Na₂S. We find that Na-salts are the major cloud species, because Na is over one order of magnitude more abundant than K or Zn. In the right panel of Figure 12, we show the transmission spectrum from visual to mid-IR band. Although scattering masks the molecular features for wavelengths λ < 1 μm, various features from CO₂ and CH₄ unveil themselves at IR wavelengths. Unlike the magnesiumsilicates on hot Jupiters, Na-salts are mostly transparent at λ > 2 μm and thus the spectrum at longer wavelength is insensitive to the formation of the Na₂S and NaCl salts.

The Hubble Space Telescope found that the warm-Neptune GJ 1214 b shows a flat spectrum from 1 to 2 μm, with an amplitude of less than 100 ppm (Kreidberg et al. 2014). The HST observation is plotted in Figure 12. We further test how ExoLyn can reproduce a similarly flatness of the spectrum, though our intention is not to quantitatively retrieve the planet property in this paper. With higher K_zz and smaller a_p0, we find that thicker clouds form and extend to the upper atmosphere. In addition, when decreasing the nucleation rate by four orders of magnitude, the particle sizes becomes larger, reaching micron size throughout the cloud. We also changed the reference pressure of the transmission spectrum calculation to 10⁻⁵ bar The combined effect of decreasing ${\dot{\mathrm{\Sigma }}}_{\text{nuc }}$ and increasing K_zz is that clouds with larger particle sizes becomes dominant in the upper atmosphere, which is consistent with a flat spectrum with low amplitude (Kempton et al. 2023).

5 Discussion

5.1 Assessment on the performance and assumptions

In the preceding sections we have demonstrated the viability of ExoLyn to calculate the cloud distribution and to obtain realistic synthetic model spectra under different planet conditions. For ExoLyn the simplifications pertain primarily the prescription of nucleation at the expense of a number of free parameters (a_p0, ${\dot{\mathrm{\Sigma }}}_n$, and σ_n), in addition to the more standard assumptions of 1D T-P profile and fixed K_zz. Under these assumptions, ExoLyn computes the mixing ratio of cloud and vapor, according to realistic condensation prescription in tandem with particle transport. Indeed, ExoLyn is efficient. On average, each complete iteration for the default hot Jupiter (Table 3), with 100 grid points, takes 2.1 seconds on a laptop PC (Intel Core^TM i7-12700H @ 4.7 GHz). In the following two paragraphs, we assess the assumptions adopted by ExoLyn, aiming to answer the question: how reliable is ExoLyn’s result.

One caveat of this work is the assumption of a single representative size for the particles at each height. Though this assumption works for hot Jupiters where nuclei form in the higher atmosphere due to, for example, nucleation of refractory species such as TiO₂ (Lee et al. 2015; Helling et al. 2017; Helling 2019) or photochemistry (Zahnle et al. 2016; Kawashima & Ikoma 2019), it may not be compatible with the nucleation condition on a terrestrial planet characterized by a surface. On Earth, nuclei are supplied from the surface by the evaporation of ocean salt particles or the eruption of volcanoes (Woods 1993; Quinn et al. 2017). These nuclei diffuse up due to turbulent motion, until growing too large due to vapor condensation and precipitate down to the surface. Therefore, the ground level is characterized by two particle sizes: one for smaller nuclei and the other for the larger rain drops. At present, such a bimodal size distribution cannot be accounted for in the current version of ExoLyn. One possible solution is to model the size of the particles with two or more bins. For example, Ohno & Okuzumi (2017) followed the smaller “cloud particles” and larger “rain particles” separately, which by their definition, have upward and downward net velocity, respectively. The smaller particles are carried up by the background upwelling vapor, until they grow heavy enough by condensation to rain down, at which point they are “promoted” to rain particles. These large particles sweep up the smaller particles during their descent. In studies that model the size distribution of particles in protoplanetary disks, such a twopopulation size distribution has also been widely applied and has proven to be in good agreement with the full size-distribution simulations (Birnstiel et al. 2012, 2015). In future works, ExoLyn can be adjusted to a similar strategy, solving the transport and formation equation of both small nuclei and large particles.

Another key assumption of ExoLyn is the single spatial dimension. In reality, horizontal variation of temperature and transport are important factors that shape the appearance and properties of clouds. For example, hot Jupiters are usually tidally locked, with a hot day-side and a cooler night-side. It is found that due to lower irradiation intensity, clouds are prone to form in the higher latitude region or the night side (Lines et al. 2018; Roman & Rauscher 2019; Roman et al. 2021) As a result, caution should be exercised when modeling tidally locked planet spectra or phase curves with a horizontally isothermal temperature model. Alternative to modeling clouds in 2D or 3D, which takes significantly more computational resources, one can set up a collection of 1 D vertical grids spanning the equatorial plane or the entire surface (Webber et al. 2015; Feng et al. 2020; Helling et al. 2020; Samra et al. 2023). With independent temperature profiles for each vertical model, ExoLyn can generate a 2D or 3D cloud profile resulting in a more realistic spectrum. In a similar vain, the horizontal transport of cloud particles and vapor could be more efficient than vertical diffusion, as is suggested by selfconsistent 2D simulation (Powell & Zhang 2024). Therefore, the horizontal transport of vapor replenishes the cloud-depleted region, increasing the cloud formation efficiency. The horizontal advection can also transport cloud particles to hot regions where 1D models would predict clouds to be absent. Similar to Powell & Zhang (2024)’s approach, a collection of ExoLyn clouds on the equatorial plane can be post-processed with horizontal advection, under a velocity characterizing the longitudinal jet. Though not fully modeling the 3D nature, this approach serves as a reasonable first-order approximation that empowers the atmosphere retrieval with horizontal transport of the cloud.

5.2 The need for a self-consistent clouds in retrievals

Several works have modeled cloud formation with a higher degree of physical consistency. For example, DRIFT-PHOENIX (Witte et al. 2009) considers nucleation rates consistent with the the vapor abundance and calculates gas phase equilibrium chemistry on the fly. Following this approach, DIFFUDRIFT (Woitke et al. 2020) accounts for vertical transport of dust and gas similar to ExoLyn. Besides, CARMA (Turco et al. 1979; Toon et al. 1979) solves the particle size distribution with ~10 mass bins. Though designed for higher precision, the qualitative structure of clouds in hot Jupiters obtained by these codes is very similar to the computations by ExoLyn: Mg-Si at the top, Fe below it, and Al₂O₃ in the lower regions of the cloud. Despite the similar results, however, the above codes are usually too time-consuming to be directly used in atmosphere retrieval. One DRIFT-PHOENIX run can take 2 CPU minutes and one DIFFUDRIFT run can take 500 CPU hours (Woitke et al. 2020). Therefore, these models are typically used as pre-calculated forward model grids (e.g., Samra et al. 2020; Gaidos & Hirano 2023; Yang et al. 2024). In contrast, ExoLyn has been designed to trade CPU-consuming precision for efficiency, rendering it suitable for retrievals.

With ExoLyn, we are able to retrieve the transmission spectrum in a self-consistent way. As clouds can suppress molecular features and change the molecular abundances in the upper atmosphere, cloud-free retrieval results may not represent the bulk abundances of the planet. Specifically, the molecular abundances may be underestimated, or the mean molecular weights may be overestimated. In addition, retrievals where parameterized clouds are put in by hands introduce excessive parameters, which cause degeneracies among the parameters (e.g., Burningham et al. 2017), resulting in poorly constrained molecular abundances. A joint retrieval, where cloud formation is modeled simultaneously with the gas components in a physical model can lift the degeneracy. With an ever-increasing number of multi-band spectra on planet atmospheres, ExoLyn allows us to retrieve the cloud and gas simultaneously, providing a physically consistent posterior on metallicity and C/O ratio.

The multi species feature of ExoLyn also opens the opportunity to constrain the composition of cloud particles. This is because solid particles of different composition have different optical properties at mid-IR wavelength (Grant et al. 2023; Dyrek et al. 2024). For example, SiO₂ displays an additional absorption feature at 12 μm compared to MgSiO₃ (see Section 3.4.4 and Figure 10). Recently, JWST / MIRI observations on WASP-17 b directly observed silicate features at around 10 μm wavelength. The silicate peak in the spectrum clearly collocates with SiO₂ in quartz form (Grant et al. 2023), although a mixture with other silicates is not ruled out. A physically consistent cloud model is necessary to find out which cloud species form under the bulk condition of these hot planets and offers the opportunity to find a unique solution.

5.3 Future directions for model improvement

There are two directions where ExoLyn can be extended toward greater consistency. First, this paper calculates the T-P profile from the Guillot (2010) model which assumes a two-stream approach with fixed inward (visible) and outward (IR) opacities. However, cloud formation alters the optical properties and hence the temperature structure in the atmosphere (Roman & Rauscher 2019; Harada et al. 2021; Roman et al. 2021). Accounting for the cloud opacity during radiative transfer can influence the atmosphere temperature. For example, on the one hand, Roman & Rauscher (2019) found that scattering by cloud particles has an effect of cooling down the atmosphere. On the other hand, the absorption by cloud particles generally warms up the atmosphere (Phillips et al. 2020; Roman & Rauscher 2019).

Second, ExoLyn at present does not account for the energetics involved in phase-change processes. On colder planets, cloud formation releases latent heat and would change the temperature profile (Johansen et al. 2021; Ge et al. 2024). The fast computation of a cloud profile by ExoLyn paves the way for a retrieval coupling the T-P profile with clouds (Min et al. 2020). In a future work, iteration between cloud formation and thermal structure calculation through radiative transfer will be performed to comprehensively understand how clouds change the planet spectrum.

6 Conclusions

We have developed ExoLyn: an exoplanet cloud model, specifically tailored for atmospheric spectrum retrieval. ExoLyn solves the formation and transport equations of multispecies cloud condensates, vapor, and particle number (nuclei), under turbulent diffusion and sedimentation. For condensation, we consider actual chemical reactions that lead to the formation of cloud particles and the depletion of vapor. With an averaged computational time of ≲2 CPU s for 100 grid points, ExoLyn is ideal for MCMC simulations in atmospheric retrieval. We apply ExoLyn to hot Jupiters, self-luminous planets, and sub-Neptunes. Our conclusions are as follows:

1. On a hot Jupiter planet three compositionally distinct cloud layers form from top down. These are: the magnesiumsilicate cloud layer, the iron cloud layer, and the refractory cloud layer. The magnesiumsilicate layer is mainly composed of a mixture of MgSiO₃, SiO₂, and Mg₂SiO₄. The refractory layer is dominated by Al₂O₃. The layered structure is consistent with the condensation sequence of the condensate species;

2. The cloud thickness and extent of the cloud are chiefly controlled by the diffusivity (K_zz parameter) of the atmosphere. The stronger the diffusion, the thicker the clouds, the larger cloud particles grow, and the stronger cloud layers mix and overlap. The nucleation properties – which we prescribe in the model – has weak effects on the cloud profile. In particular, the effect of increasing the nucleation rate is counteracted by particle coagulation;

3. For hot-Jupiters at diffusivity levels of K_zz = 10⁸ cm²s⁻¹, iron and refractory clouds are not visible because they are located in the optically thick region in the atmosphere. When thicker clouds form at higher diffusivity, the molecular lines (H₂O and CO) are suppressed, but the magnesiumsilicate features at 8–10 μm become more prominent. On the other hand, when K_zz ⩽ 10⁴ cm² s⁻¹, clouds disappear from the IR transmission spectra (Figure 9), because thin clouds form in the molecular optically thick layer;

4. Increasing metallicity leads to significant changes in the composition of cloud particles and vapor components. The formation of Mg₂SiO₄ at a depth of 1 bar depletes the Mg vapor concentration in the upper atmosphere. As a result, the dominant cloud species becomes SiO₂ rather than MgSiO₃. This transition is reflected in an additional SiO₂ feature at 12 μm. Due to higher molecular abundances, metalrich planets exhibit stronger molecular lines, even with the presence of cloud;

5. For a hot Jupiter planet with C/O ratio close to unity, iron clouds dominate because the low H₂O abundance suppresses the formation of O-bearing species. The iron cloud increases the transit depth, suppresses the molecular lines, and the Mg-Si features disappears;

6. ExoLyn is also applicable to self-luminous planets and sub-Neptunes. Clouds that form on self-luminous planets redden the spectra and suppress the molecular lines. On sub-Neptunes Na₂S and NaCl clouds form, whose scattering obscures the molecular features below 2 μm. ExoLyn predicts a flat spectrum in the near-IR when a modest amount of nuclei form in a turbulent atmosphere (high K_zz), consistent with the observation of the sub-Neptune GJ 1214 b.

Data availability

The code ExoLyn is available at https://github.com/helonghuangastro/exolyn.

Acknowledgements

This work is supported by the National Science Foundation of China under grant Nos. 12250610189 and 12233004. The authors thank the anonymous referee for their comments, which have improved the quality of the manuscript. The authors gratefully thank Kazumasa Ohno, Peter Gao, Michael Zhang, Xi Zhang, and Paul Mollière for insightful discussion.

References

All Tables

All Figures

Fig. 1 Comparison of ExoLyn’s multispecies approach to condensation against the single-species approach of Ormel & Min (2019). ExoLyn considers the actual vapor species that are present in the atmosphere. As illustrated in the right panel, the condensation rate of MgSiO3 is determined by the amount of Mg – the species that limits the condensation by its vapor concentration and stoichiometric number in the reaction.

In the text

Fig. 2 Default temperature-pressure profile of the hot Jupiter (Section 3), self-luminous (Section 4.1), and sub-Neptune planet (Section 4.2), following Guillot (2010). The irradiation temperatures of the hot Jupiter, self-luminous planet, and the sub-Neptune planet are taken to be 500 K, 1000 K, and 300 K, correspondingly. Condensation curves of Mg2SiO4, SiO2, Fe, Al2O3 (for hot Jupiter and self-luminous planet), NaCl and ZnS (for sub-Neptune) are shown as colored lines. The solid, dashed, and dotted saturation curves denote solar, 10 times solar, and 100 times solar metallicity, respectively.

In the text

Fig. 3 Atmosphere structure of the default hot Jupiter model. The thick colored line shows the cloud mass mixing ratio summed over all species. The dominant species by mass is illustrated by the color of the thick line. Vapor mixing ratios are shown by dashed lines. The dotted black line denotes the nuclei mass mixing ratio and the purple thin line denotes the cloud particle size, which should be referred to the upper axis. In the figure we also label the evaporation temperature of magnesiumsilicates (1460 K) and iron (1600 K) with horizontal dashed lines. The gray band refers to where the cloud optical depth in the IR band (1 μm − 20 μm), integrated from the top of the atmosphere, reaches 1. The right panel shows the composition of cloud particles in mass.

In the text

Fig. 4 Comparison of the full ExoLyn model (default), ExoLyn restricted with MgSiO3 condensates only, the Ormel & Min (2019) model, and Ackerman & Marley (2001). The thickness of the lines and the symbols on the lines indicates the particle size, except in the Ackerman & Marley (2001) runs, where the particle size is not solved for as it decouples with the cloud concentration in the model. Right: the same as left, but Kzz decreased to 107 cm2s−1. In this case, the fsed = 2 run (thin line) provides a better correspondence to other models.

In the text

Fig. 5 Cloud profiles under different eddy diffusivity parameters Kzz, which serve as effective vertical turbulence. Left: total cloud mixing ratio (solid lines) and Fe vapor mixing ratio (dashed lines) as an example. Right: nuclei mixing ratio (solid lines) and particle size (dashed lines). Similar to Figure 3, the location where the cloud optical depth reaches unity is shown in the left panel by the horizontal bands with corresponding colors. We note that there is no blue band corresponding to Kzz = 104 cm2s−1, because the cloud optical depth never reaches 1.

In the text

Fig. 6 Sensitivity of the nucleation profiles on the cloud profiles. In the figure we changed the nucleation strength ${\dot{\mathrm{\Sigma }}}_n$, nucleation location Pn, spread of nucleation σn, and nuclei size ap0. The parameters are changed with respect to the default value in Table 1. The vertical lines indicate the location of nucleation, with its length meaning 1σ dispersion of the nucleation profile.

In the text

	Fig. 7 Cloud profiles with different condensation reactions. The full set of reactions Table 3 is shown as the black line. The green and yellow line only includes the magnesiumsilicate group versus the iron and refractory reaction group, that is, the first four rows and the last six rows in Table 3. The vertical lines represent where the near-IR optical depth accumulates to unity.
In the text

Fig. 8 Effect of vapor abundances at the lower boundary of atmosphere on the solid, vapor, nuclei profile, particle size, and cloud composition. (a) 10×-Z hot Jupiter model; (b) C/O=0.9 hot Jupiter model; (c) C/O=0.9 hot Jupiter model, similar to (b), but with injection of larger nuclei deeper in the atmosphere, mimicking nucleation featuring high particle surface tension. The thick solid, dashed, and dotted lines indicate total cloud, vapor, and nuclei mixing ratios, respectively. The purple solid line shows the particle size. The right subpanel in each panel shows the composition of the cloud particles in mass. Compared to the default model (Figure 3) the upper cloud layers of the 10×-[Z] run become devoid in Mg (SiO2 dominates) while in the C/O=0.9 run Fe-bearing cloud particles make it into the upper regions.

In the text

Fig. 9 Transmission spectrum of a cloudy hot Jupiter (black line). The gray line shows the transmission spectrum of a clear atmosphere, leaving the gas components unchanged. Kzz is varied to show how cloud formation changes the atmospheric transmission spectrum. The horizontal lines and arrows indicate molecular features. At long wavelengths the Kzz = 104 cm2 s−1 run curves are indistinguishable from the “clear” curve.

In the text

Fig. 10 Transmission spectrum of models varying metallicity and C/O ratio. The dark and light line colors refer to cloudy and clear atmosphere, respectively. The black and gray lines are the default cloudy and clear model (Section 3.1). 10×-[Z] and C/O=0.9 cases are shown by the blue lines and the red lines. For the C/O=0.9 run, we additionally conduct a run injecting larger (5 nm) nuclei particles deeper in the atmosphere to accounting for the Kelvin effect (brown curve).

In the text

Fig. 11 Cloud structure and emission spectrum of a self-luminous planet with physical properties similar to HR 8799 e. The left panel shows the cloud, gas, and nuclei mixing ratios, with every components the same meaning as Figure 8. The brown thick and thin lines present the temperature profiles we used in the clear and cloudy models, respectively. The right panel shows the synthetic emission spectrum of a cloudy atmosphere (the black line) and a clear atmosphere (the blue line), whose effective temperatures are both 1000 K. For reference, the spectrum of a 1000 K black body radiation is also plotted as the dashed line.

In the text

Fig. 12 Cloud structure and transmission spectrum of a sub-Neptune planet with physical properties similar to GJ 1214 b. The condensate species considered here is NaCl, KCl, ZnS, and Na2S. In the left panel, the meaning of each line is the same as Figure 8. The cloud is dominated by Na2S (upper part) and NaCl (cloud bottom). The right panel shows the synthetic spectrum and the blue line shows one set of parameters that leads to flat spectrum in 1 to 2 μm range, consistent with observation (Kreidberg et al. 2014).

In the text