© The Authors 2025

1 Introduction

The size distribution of cloud particles plays a critical role in shaping the observable properties of sub-stellar atmospheres and significantly influences the thermal balance through scattering and absorption opacities, which depend on the composition of the condensing material. Recent evidence for complex particle size distributions in sub-stellar atmospheres includes the inference of nanometre-sized particles in WASP-17b (Grant et al. 2023), consistent with a solid SiO₂ composition based on infrared absorption features. Retrieval modelling of brown dwarf atmospheres (Burningham et al. 2021) also suggests a cloud profile comprising various particle sizes, with ~0.1 μm silicate particles being dominant and lofted to high altitudes.

Typically, microphysical models of cloud formation in substellar atmospheres use equation sets that assume a monodisperse distribution, utilising a single representative average value (e.g. Woitke & Helling 2003; Helling et al. 2008; Ohno & Okuzumi 2018; Lee 2023). However, the monodisperse formulation neglects the effects of polydispersity on cloud formation processes such as condensation, evaporation, and collisional growth, where the variety of particle sizes can significantly influence the overall cloud structure. Simulations using bin-resolving models such as the CARMA model (e.g. Gao et al. 2018; Powell et al. 2019) applied to hot Jupiter and brown dwarf atmospheres suggest a range of particle sizes can be present in sub-stellar atmospheres, as well as potential multi-modality. These results highlight the importance of accounting for polydispersity when modelling cloud formation in sub-stellar atmospheres.

Moment (or bulk) cloud microphysical methods offer a useful middle ground between the more computationally expensive bin (or spectral) models (e.g. Gao et al. 2018; Powell & Zhang 2024) and simplified schemes such as tracer saturation adjustment (e.g. Tan & Showman 2019; Komacek et al. 2022). Consequently, the moment method is frequently coupled to large-scale 3D atmospheric simulations, such as general circulation models (GCMs), in the Earth science community (e.g. Morrison & Gettelman 2008), hot Jupiter exoplanets (e.g. Lee et al. 2016; Lines et al. 2018; Lee 2023) and brown dwarf atmospheres (e.g. Lee & Ohno 2025). However, to date, only monodisperse moment methods have been employed in GCM simulations. While these can predict a representative particle size, they cannot self-consistently predict the width of the size distribution without additional assumptions. Other schemes use additional assumptions on the properties of the size distribution and particle sizes, such as time-independent equilibrium balance between vertical transport and settling velocity (Ackerman & Marley 2001), or localised, instantaneous formation or rain out (e.g. Parmentier et al. 2016; Roman & Rauscher 2019). These approaches generally do not evolve a tracer within the GCM, but instead diagnose cloud properties directly from the local temperature and chemical structure of each vertical profile. Of particular relevance to the current study is the work of Christie et al. (2022), who derived a variant of the Ackerman & Marley (2001) model that assumes a gamma distribution rather than the commonly used lognormal.

To overcome this limitation, several studies in the Earth cloud science literature adopt the three-moment method that uses the shape of the size distribution as a prognostic variable (e.g. Milbrandt & Yau 2005; Naumann & Seifert 2016; Paukert et al. 2019). In this follow-up study to Lee & Ohno (2025) and Lee (2025) (Paper I), which presented two-moment methods with a monodisperse and exponential size distribution, respectively, we extended the model to a three-moment method that is capable of predicting the evolution of size distribution shape. Due to the uncertainty in the exact size distribution properties present in sub-stellar atmospheres, we attempted to derive size distribution properties self-consistently with minimal parameterisations. As often assumed in Earth cloud models, we adopted the gamma particle size distributions to describe the underlying size distribution of cloud particles. In Section 2, we introduced the observational and theoretical motivation to assume the gamma size distribution. In Section 2.3, we presented the properties of the gamma distribution, its moments and introduce the three-moment microphysical scheme. In Section 3, we derived the condensation and evaporation rates for the gamma distribution. In Section 4, we derived collisional growth rates, taking into account the polydispersity introduced by the gamma distributions. In Section 5, we applied our new methodology to the same Y-dwarf KCl cloud formation case as in Paper I. Sections 6 and 7 contain the discussion and conclusions, respectively, of our study.

2 Introduction of gamma distribution

The gamma particle mass distribution, f (m) [cm⁻³ g⁻¹], is given by $$f(m)=\frac{N_{\mathrm{c}}}{{\lambda }^{\nu }\mathrm{\Gamma }(\nu )}{m}^{\nu -1}{e}^{-m/\lambda },$$(1)

where N_c [cm⁻³] is the total number density of the distribution, λ [g] is the scale parameter, ν the shape parameter, $\mathrm{\Gamma }(z)\equiv {\int }_0^{\mathrm{\infty }}{t}^{z-1}{e}^{-t}\mathrm{d}t$ is the gamma function, and m [g] is the mass of an individual particle. When ν = 1, the gamma distribution reduces to an exponential distribution, which was explored in Paper I. λ and ν are related to the mean mass (expectation value E[m] [g]) and variance (Var[m] [g²]) of the distribution through $$\mathrm{E}[m]=\lambda \nu ,$$(2)

and $$\mathrm{V}\mathrm{a}\mathrm{r}[m]={\lambda }^2\nu ,$$(3)

respectively. In Figure 1, we illustrated how the gamma distribution changes with shape parameter ν. In general, the gamma distribution consists of a power-law tail at smaller particle sizes with an exponential cut off at large sizes. As ν → 0, the distribution widens, becoming more polydisperse, while as ν → ∞ the distribution narrows, showing more monodisperse characteristics.

Fig. 1 Gamma size distribution compared to the exponential, Rayleigh and monodisperse distributions for various values of ν. The representative particle size is rc = 1 μm and total number density is Nc = 1 cm−3.

2.1 Observational motivation

There are several observational motivations for adopting the gamma distribution. A seminal work of Marshall & Palmer (1948) showed that the size distribution of rain droplets in Earth’s atmosphere can be well fitted by an exponential function, commonly known as the Marshall-Palmer distribution¹, except at small particle sizes. Later studies showed that the gamma distribution is sufficiently versatile to fit the observed cloud properties (e.g. Ulbrich 1983; McFarquhar et al. 2015). Beyond Earth clouds, Hansen & Hovenier (1974) showed that sulfuric acid clouds modelled with the Hansen distribution², a variant of the gamma distribution, can reproduce the observed polarisation signatures of Venusian clouds. In brown dwarf atmospheres, Hiranaka et al. (2016) showed that the gamma distribution containing abundant sub-micron particles can well explain the near-infrared colours of some red L dwarfs. Retrieval modelling by Burningham et al. (2021) also found that a gamma distribution is strongly preferred over a lognormal distribution when fitting the emission spectrum of the red L dwarf 2M2224-0158.

2.2 Physical motivation

In addition to observational evidence, there are several physical motivations for adopting the gamma distribution. The fundamental equation describing the evolution of particle size distribution is given by (see, e.g. Appendix A of Lee & Ohno 2025) $$\mathrm{d}n/\mathrm{d}D=N_0\exp (-{\lambda }_{\mathrm{D}}D),$$(6)

where z [cm] is the altitude, ρ_a [g cm⁻³] is the atmospheric density, v_f [cm s⁻¹] is the terminal settling velocity, (dm/dt) [g s⁻¹] the rate of change in particle mass due to microphysical processes such as condensation and K_zz [cm² s⁻¹], the eddy diffusion coefficient. For illustrative purposes here, we have ignored collisional growth and adopt a 1D framework that approximates vertical transport by atmospheric circulation as an eddy diffusion term.

To gain physical insights into the gamma distribution, we conducted an order-of-magnitude analysis of Eq. (6). Focusing on a single vertical layer, the vertical transport sink term acts to remove cloud particles from the layer, while the second term acts to increase or decrease the mass of the particles inside the layer. Thus, one may crudely approximate ∂/∂z ~ 1/H_a, where H_a [cm] is the atmospheric altitude length scale, to rewrite the equation as $$\frac{\mathrm{\partial }}{\mathrm{\partial }m}[\left(\frac{\mathrm{d}m}{\mathrm{d}t}\right)f(m)]\sim -(\frac{K_{\mathrm{z}\mathrm{z}}}{H_{\mathrm{a}}^2}+\frac{v_{\mathrm{f}}(m)}{H_{\mathrm{a}}})f(m),$$(7)

where we have considered a steady state. As shown later, the mass growth rate and settling velocity can be often expressed as a power-law function of mass. Let (dm/dt) = Am^α and v_f(m) = Bm^β, we can further rewrite the equation as $$\frac{\mathrm{\partial }}{\mathrm{\partial }m}[m^{\alpha }f(m)]\sim -\frac{1}{A{m}^{\alpha }}(\frac{K_{\mathrm{z}\mathrm{z}}}{H_{\mathrm{a}}^2}+\frac{B{m}^{\beta }}{H_{\mathrm{a}}})m^{\alpha }f(m).$$(8)

One can analytically solve this equation to obtain the steady-state size distribution as $$f(m)\propto m^{-\alpha }\exp [-(\frac{K_{\mathrm{z}\mathrm{z}}}{(1-\alpha )H_{\mathrm{a}}^2}+\frac{B{m}^{\beta }}{(1+\beta -\alpha )H_{\mathrm{a}}})\frac{m^{1-\alpha }}{A}].$$(9)

The obtained analytic distribution consists of a power-law tail at small particle mass with an exponential cutoff at large mass, reproducing the salient characteristics of the gamma distribution. This analysis indicates that the power-law tail at low particle mass encapsulates the information of mass dependence of particle growth. By introducing a particle mass growth timescale τ_grow = m/(dm/dt) = m^1-α/A [s], an eddy diffusion timescale τ_edd = H²_a/K_zz [s], and a settling timescale τ_set = H_a/v_f = H_a/Bm^β [s], one can re-write the analytic distribution into a physically informative form $$f(m)\propto m^{-\alpha }\exp [-(\frac{{\tau }_{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}}}{(1-\alpha ){\tau }_{\mathrm{e}\mathrm{d}\mathrm{d}}}+\frac{{\tau }_{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}}}{(1+\beta -\alpha ){\tau }_{\mathrm{s}\mathrm{e}\mathrm{t}}})].$$(10)

Thus, the distribution follows a power-law behavior at masses at which the growth is faster than vertical transport, i.e. τ_grow ≪ τ_edd and τ_grow ≪ τ_set, whereas the distribution exponentially drops once the vertical transport is faster than the growth. Note that Garrett (2019) presented a similar analysis to provide a physical reason why rain droplets in the Earth closely follow the gamma distribution. Although the derivation presented above is potentially overly simplistic to accurately predict the real size distribution, it provides a natural physical motivation to adopt the gamma distribution when developing a polydisperse particle distribution model including cloud microphysics.

2.3 Moments with assumed gamma distribution

We present a derivation of the mass moments assuming a gamma distribution, along with the source terms required for evolving the cloud microphysics ordinary differential equation (ODE) system in time. The mass moments of the particle mass distribution, M^(k) [g^k cm⁻³], are given by $$M^{(k)}={\int }_0^{\mathrm{\infty }}{m}^{k}f(m)\mathrm{d}m,$$(11)

where k is the order of the moment, which can take an integer or non-integer value. Inserting the gamma distribution in the moment definition equation gives $$M^{(k)}=\frac{N_{\mathrm{c}}}{{\lambda }^{\nu }\mathrm{\Gamma }(\nu )}{\int }_0^{\mathrm{\infty }}{m}^{k+\nu -1}{e}^{-m/\lambda }\mathrm{d}m.$$(12)

After some algebra and applying the definition of the gamma integral, the moment generator for the gamma mass distribution is analytically expressed as $$M^{(k)}=N_{\mathrm{c}}{\lambda }^k\frac{\mathrm{\Gamma }(\nu +k)}{\mathrm{\Gamma }(\nu )}.$$(13)

Each integer moment represents bulk values of the particle size distribution, such as the total number density, N_c [cm⁻³], for k = 0 $$M^{(0)}=N_{\mathrm{c}},$$(14)

the k = 1 moment that represents the total mass density, ρ_c [g cm⁻³], of the cloud particles $$M^{(1)}={\rho }_{\mathrm{c}},$$(15)

and the k = 2, Z_c [g² cm⁻³], moment $$M^{(2)}=Z_{\mathrm{c}},$$(16)

that is related to the Rayleigh scattering cross-section of the particles.

The parameters λ and ν can be estimated from the computed integer moments. The number-weighted mean mass, m_c [g], is given by the ratio of the first to zeroth moment $$m_{\mathrm{c}}=\frac{{\rho }_{\mathrm{c}}}{N_{\mathrm{c}}}=\lambda \nu .$$(17)

The variance, σ²_c [g²], is given by $${\sigma }_{\mathrm{c}}^2=\frac{Z_{\mathrm{c}}}{N_{\mathrm{c}}}-{\left(\frac{{\rho }_{\mathrm{c}}}{N_{\mathrm{c}}}\right)}^2={\lambda }^2\nu .$$(18)

Rearranging gives expressions for the shape and scale parameters $$\nu =\frac{m_{\mathrm{c}}^2}{{\sigma }_{\mathrm{c}}^2},$$(19)

and $$\lambda =\frac{m_{\mathrm{c}}}{\nu },$$(20)

respectively. The representative particle radius, r_c [cm], corresponding to the number-weighted mean mass, m_c, is given by $$r_{\mathrm{c}}={\left(\frac{3m_{\mathrm{c}}}{4\pi {\rho }_{\mathrm{d}}}\right)}^{1/3},$$(21)

where ρ_d [g cm⁻³] is the particle bulk density. Hereafter, we refer to r_c as the ‘representative radius’. Similarly, the mass-weighted mean mass can be defined as $$m_{\mathrm{p}}=\frac{{\int }_0^{\mathrm{\infty }}m\cdot mf(m)\mathrm{d}m}{{\int }_0^{\mathrm{\infty }}mf(m)\mathrm{d}m}=\frac{M^{(2)}}{M^{(1)}}=\lambda (\nu +1).$$(22)

This expression highlights that the number-weighted mean mass, m_c, is representative in the ν ≫ 1 regime, whereas the mass-weighted mean mass, m_p, deviates from m_c for ν ≪ 1, where the distribution becomes broader and polydispersity plays a more significant role. The representative radius of the massdominating particles, r_p [cm], is given by $$r_{\mathrm{p}}={\left(\frac{3m_{\mathrm{p}}}{4\pi {\rho }_{\mathrm{d}}}\right)}^{1/3}=r_{\mathrm{c}}{\left(\frac{\nu +1}{\nu }\right)}^{1/3}.$$(23)

In this study, we considered only spherical, homogeneous cloud particles, allowing Eqs. (21) and (23) to be applied. Porous particles and fractal geometries, as considered in some previous studies (e.g. Ohno et al. 2020; Samra et al. 2020), are not included.

3 Condensation with gamma distribution

In this section, we derive the condensation and evaporation rates for a cloud particle population described by a gamma distribution. Unlike the exponential case (Paper I), the gamma formulation requires additional expressions involving the second moment, in addition to the first. In the method of moments (e.g. Gail & Sedlmayr 2013), the rate of change of the k-th moment, dM^(k)/dt, due to condensation or evaporation is given in Paper I. Following the definitions of the Knudsen number, Kn, and population weighted Knudsen number from Paper I., for three moment scheme, we required a weighting of m² in addition to the number and mass weighted equations presented in Paper I $$\overline{{\mathrm{K}\mathrm{n}}_{m^2}}=\frac{{\lambda }_{\mathrm{a}}}{r_{\mathrm{c}}}{\nu }^{1/3}\frac{\mathrm{\Gamma }(\nu +5/3)}{\mathrm{\Gamma }(\nu +2)}.$$(24)

Figure 2 shows how the population-averaged Knudsen numbers vary with ν. Notably, the number-weighted formulation exhibits singularities for ν ≤ 1/3 due to the divergence of the gamma function. To resolve this, we adopted an incomplete gamma formulation (Appendix A), which avoids the singularities and is implemented in the model code. The resulting solution is also shown in Figure 2. As ν becomes larger and the distribution narrows, all three population averaged Knudsen numbers converge towards the monodisperse value.

In the continuum (diffusive) regime (Kn ≪ 1), the condensation or evaporation rate, following the methods presented in Paper I with the gamma distribution, is given by $$\frac{\mathrm{d}{M}^{(k)}}{\mathrm{d}t}=k{C}_0{\lambda }^{k-2/3}{N}_{\mathrm{c}}\frac{\mathrm{\Gamma }(\nu +k-2/3)}{\mathrm{\Gamma }(\nu )}.$$(25)

We then used Eq. (13) again for the k - 1 moment to factor out N_c and get the expression in terms of the previous moment $$\frac{\mathrm{d}{M}^{(k)}}{\mathrm{d}t}=k{C}_0{\lambda }^{1/3}{M}^{(k-1)}\frac{\mathrm{\Gamma }(\nu +k-2/3)}{\mathrm{\Gamma }(\nu +k-1)},$$(26)

where C₀ is a kinetic pre-factor given in Paper I. We then used the relation λ^1/3 = m_c^1/3v^−1/3 to reintroduce the representative particle size, r_c , into the equation. The full expression for the first moment (k = 1) is then $${\left(\frac{\mathrm{d}{\rho }_{\mathrm{c}}}{\mathrm{d}t}\right)}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}=4\pi r_{\mathrm{c}}D{m}_0{n}_{\mathrm{v}}(1-\frac{1}{S}){\nu }^{-1/3}\frac{\mathrm{\Gamma }(\nu +1/3)}{\mathrm{\Gamma }(\nu )}{N}_{\mathrm{c}},$$(27)

and the second moment (k = 2) as $${\left(\frac{\mathrm{d}{Z}_{\mathrm{c}}}{\mathrm{d}t}\right)}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}=8\pi r_{\mathrm{c}}D{m}_0{n}_{\mathrm{v}}(1-\frac{1}{S}){\nu }^{-1/3}\frac{\mathrm{\Gamma }(\nu +4/3)}{\mathrm{\Gamma }(\nu +1)}{\rho }_{\mathrm{c}}.$$(28)

In the free molecular (or kinetic) regime (Kn ≫ 1), similar to the Kn ≪ 1 regime solution, using Eq. (13) the final expression is given as $$\frac{\mathrm{d}{M}^{(k)}}{\mathrm{d}t}=k{C}_0{\lambda }^{k-1/3}{N}_{\mathrm{c}}\frac{\mathrm{\Gamma }(k+\nu -1/3)}{\mathrm{\Gamma }(\nu )},$$(29)

and in terms of the k - 1 moment $$\frac{\mathrm{d}{M}^{(k)}}{\mathrm{d}t}=k{C}_0{\lambda }^{2/3}{M}^{(k-1)}\frac{\mathrm{\Gamma }(k+\nu -1/3)}{\mathrm{\Gamma }(k+\nu -1)},$$(30)

where C₀ is a kinetic pre-factor given in Paper I. As with the Kn ≪ 1 regime, using the relation λ^2/3 = m_c^2/3v^−2/3 returns the representative particle size, r_c, to the expression. This then gives for the full first moment (k = 1) equation $${\left(\frac{\mathrm{d}{\rho }_{\mathrm{c}}}{\mathrm{d}t}\right)}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}}^{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}=4\pi r_{\mathrm{c}}^2{v}_{\mathrm{t}\mathrm{h}}{m}_0{n}_{\mathrm{v}}(1-\frac{1}{S}){\nu }^{-2/3}\frac{\mathrm{\Gamma }(\nu +2/3)}{\mathrm{\Gamma }(\nu )}{N}_{\mathrm{c}},$$(31)

and second moment (k = 2) $${\left(\frac{\mathrm{d}{Z}_{\mathrm{c}}}{\mathrm{d}t}\right)}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}}^{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}=8\pi r_{\mathrm{c}}^2{v}_{\mathrm{t}\mathrm{h}}{m}_0{n}_{\mathrm{v}}(1-\frac{1}{S}){\nu }^{-2/3}\frac{\mathrm{\Gamma }(\nu +5/3)}{\mathrm{\Gamma }(\nu +1)}{\rho }_{\mathrm{c}}.$$(32)

We adopted the same tanh-based interpolation function, f (Kn′), from Paper I to smoothly transition between the Kn ≪ 1 and Kn ≫ 1 regimes.

Fig. 2 Ratio of the population averaged Knudsen numbers with the monodisperse Knudsen number. Due to the gamma function properties, a pole occurs at values of ν ≤ 1/3 for the number-weighted Knudsen number (blue solid line), leading to singularities at ν ≤ 1/3. The incomplete gamma function formulation using Eq. (A.11), which accounts for a cutoff at small particle sizes, avoids the singularities in the gamma function (dashed blue line) and is valid across the ν range.

4 Collisional growth with gamma distribution

In this section, we extend the derivation of collisional growth rates from Paper I to include expressions for the second moment in both Knudsen number regimes. Drake (1972) derived the mass moment generating form of the Smoluchowski equation for collisional growth which was presented in Paper I. We required the addition of the second (k = 2) moment expression $$\begin{gathered} \frac{\mathrm{d}{M}^{(2)}}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}K(m,m^{\prime })m{m}^{\prime }f(m)f(m^{\prime })\mathrm{d}m\mathrm{d}{m}^{\prime } \\ =\overline{K_m(m,m^{\prime })}{\rho }_{\mathrm{c}}^2, \end{gathered}$$(33)

with the mass-weighted population averaged kernel for the second moment $$\begin{array}{rl}\overline{K_m(m,m^{\prime })}& =\frac{{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}K(m,m^{\prime })m{m}^{\prime }f(m)f(m^{\prime })\mathrm{d}m\mathrm{d}{m}^{\prime }}{{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}m{m}^{\prime }f(m)f(m^{\prime })\mathrm{d}m\mathrm{d}{m}^{\prime }}\\ & =\frac{1}{{\rho }_{\mathrm{c}}^2}{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}K(m,m^{\prime })m{m}^{\prime }f(m)f(m^{\prime })\mathrm{d}m\mathrm{d}{m}^{\prime }.\end{array}$$(34)

4.1 Particle Kn ≪ 1 regime

We attempt a zeroth and second moment collisional rate solution assuming a gamma distribution following the methods outlined in Paper I. The number density weighted population averaged kernel in this regime is given Paper I, with the mass-weighted averaged kernel as $$\begin{gathered} \overline{K_m(m,m^{\prime })}=\frac{K_0}{{\rho }_{\mathrm{c}}^2}{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}(\beta m^{-1/3}+{\beta }^{\prime }{m}^{\prime -1/3})(m^{1/3}+m^{\prime 1/3}) \\ \cdot m{m}^{\prime }f(m)f(m^{\prime })\mathrm{d}m\mathrm{d}{m}^{\prime }, \end{gathered}$$(35)

with K₀ the kinetic pre-factor given in Paper I.

After some algebra, the change in the number density and second moment are then $$\begin{array}{r}{\left(\frac{\mathrm{d}{N}_{\mathrm{c}}}{\mathrm{d}t}\right)}_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{g}}\approx -\frac{2k_{\mathrm{b}}T}{3{\eta }_{\mathrm{a}}}{N}_{\mathrm{c}}^2[1+\frac{\mathrm{\Gamma }(\nu +1/3)\mathrm{\Gamma }(\nu -1/3)}{\mathrm{\Gamma }(\nu {)}^2}\\ +\frac{{\lambda }_{\mathrm{a}}}{r_{\mathrm{c}}}A{\nu }^{1/3}(\frac{\mathrm{\Gamma }(\nu -1/3)}{\mathrm{\Gamma }(\nu )}+\frac{\mathrm{\Gamma }(\nu +1/3)\mathrm{\Gamma }(\nu -2/3)}{\mathrm{\Gamma }(\nu {)}^2})],\end{array}$$(36)

and $$\begin{array}{r}{\left(\frac{\mathrm{d}{Z}_{\mathrm{c}}}{\mathrm{d}t}\right)}_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{g}}\approx \frac{4k_{\mathrm{b}}T}{3{\eta }_{\mathrm{a}}}{\rho }_{\mathrm{c}}^2[1+\frac{\mathrm{\Gamma }(\nu +4/3)\mathrm{\Gamma }(\nu +2/3)}{\mathrm{\Gamma }(\nu +1{)}^2}\\ +\frac{{\lambda }_{\mathrm{a}}}{r_{\mathrm{c}}}A{\nu }^{1/3}(\frac{\mathrm{\Gamma }(\nu +2/3)}{\mathrm{\Gamma }(\nu +1)}+\frac{\mathrm{\Gamma }(\nu +4/3)\mathrm{\Gamma }(\nu +1/3)}{\mathrm{\Gamma }(\nu +1{)}^2})],\end{array}$$(37)

respectively.

4.2 Particle Kn ≫ 1 regime

We can follow the same methods in Morán & Kholghy (2023) and Paper I to derive a population averaged kernel function for the gamma distribution in the Kn ≫ 1 regime. The additional second moment expression required for the three moment scheme is then $$\begin{gathered} \overline{K_m(m,m^{\prime })}=\frac{K_{0}H}{{\rho }_{\mathrm{c}}^2}{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}{(m^{1/3}+m^{\prime 1/3})}^2(m^{-1/2}+m^{\prime -1/2}) \\ \cdot m{m}^{\prime }f(m)f(m^{\prime })\mathrm{d}m\mathrm{d}{m}^{\prime }, \end{gathered}$$(38)

for the mass-weighted population averaged kernel, with K₀ a kinetic pre-factor given in Paper I. Here, we assumed H = 1/√2, revised from the H = 0.85 value from Paper I. Using the moment generator for the gamma mass distribution (Eq. (13)), results in $$\begin{array}{r}\overline{K(m,m^{\prime })}=2K_{0}H{\lambda }^{1/6}[\frac{\mathrm{\Gamma }(\nu +2/3)\mathrm{\Gamma }(\nu -1/2)}{\mathrm{\Gamma }(\nu {)}^2}\\ +\frac{2\mathrm{\Gamma }(\nu +1/3)\mathrm{\Gamma }(\nu -1/6)}{\mathrm{\Gamma }(\nu {)}^2}+\frac{\mathrm{\Gamma }(\nu +1/6)}{\mathrm{\Gamma }(\nu )}],\end{array}$$(39)

and $$\begin{array}{r}\overline{K_m(m,m^{\prime })}=2K_{0}H{\lambda }^{1/6}[\frac{\mathrm{\Gamma }(\nu +5/3)\mathrm{\Gamma }(\nu +1/2)}{\mathrm{\Gamma }(\nu +1{)}^2}\\ +\frac{2\mathrm{\Gamma }(\nu +4/3)\mathrm{\Gamma }(\nu +5/6)}{\mathrm{\Gamma }(\nu +1{)}^2}+\frac{\mathrm{\Gamma }(\nu +7/6)}{\mathrm{\Gamma }(\nu +1)}].\end{array}$$(40)

Using the relation ${\lambda }^{1/6}=m_{\mathrm{c}}^{1/6}{\nu }^{-1/6}=m_{\mathrm{c}}^{2/3}{m}_{\mathrm{c}}^{-1/2}{\nu }^{-1/6}$, the representative particle size, r_c, is returned, and the final expression for the change in number density is $$\begin{array}{r}{\left(\frac{\mathrm{d}{N}_{\mathrm{c}}}{\mathrm{d}t}\right)}_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{g}}\approx -H\sqrt{\frac{8\pi k_{\mathrm{b}}T}{m_{\mathrm{c}}}}{r}_{\mathrm{c}}^2{N}_{\mathrm{c}}^2{\nu }^{-1/6}[\frac{\mathrm{\Gamma }(\nu +2/3)\mathrm{\Gamma }(\nu -1/2)}{\mathrm{\Gamma }(\nu {)}^2}\\ +\frac{2\mathrm{\Gamma }(\nu +1/3)\mathrm{\Gamma }(\nu -1/6)}{\mathrm{\Gamma }(\nu {)}^2}+\frac{\mathrm{\Gamma }(\nu +1/6)}{\mathrm{\Gamma }(\nu )}],\end{array}$$(41)

and for the second moment $$\begin{array}{r}{\left(\frac{\mathrm{d}{Z}_{\mathrm{c}}}{\mathrm{d}t}\right)}_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{g}}\approx 2H\sqrt{\frac{8\pi k_{\mathrm{b}}T}{m_{\mathrm{c}}}}{r}_{\mathrm{c}}^2{\rho }_{\mathrm{c}}^2{\nu }^{-1/6}[\frac{\mathrm{\Gamma }(\nu +5/3)\mathrm{\Gamma }(\nu +1/2)}{\mathrm{\Gamma }(\nu +1{)}^2}\\ +\frac{2\mathrm{\Gamma }(\nu +4/3)\mathrm{\Gamma }(\nu +5/6)}{\mathrm{\Gamma }(\nu +1{)}^2}+\frac{\mathrm{\Gamma }(\nu +7/6)}{\mathrm{\Gamma }(\nu +1)}].\end{array}$$(42)

In Eqs. (36) and (41), negative arguments to the gamma function are possible when ν is accompanied with a negative fractional scaling. This leads to singularities in the solution, similar to the number-weighted Knudsen number (Fig. 2). In Appendix A, we provided a physical explanation for this diverging behaviour as well as a derivation utilising the incomplete gamma function which avoids the singularities when using the gamma function. This incomplete gamma function solution is ultimately used in the numerical model. We used the diffusive Knudsen number interpolation function similar to Paper I to calculate coagulation rates in the intermediate Knudsen number regime (Kn ~ 1) for the zeroth and second moments population averaged kernels.

4.3 Gravitational coalescence

We followed a similar analysis from Paper I for calculating the gravitational coalescence rates. For a gamma distribution, we made the assumption that ∆v_f and E are well defined at their respective representative particle size values for the zeroth and second moment, allowing them to be taken outside the integral. However, we note this assumption may induce significant inaccuracy in the solution since ∆v_f and E vary by a large amount dependent on the individual particle mass pairs, especially for broad distributions. Solving the integral from Paper I for the second moment results in $$\overline{K_m(m,m^{\prime })}=\frac{2K_0}{{\rho }_{\mathrm{c}}^2}(M^{(5/3)}{M}^{(1)}+M^{(4/3)}{M}^{(4/3)}),$$(43)

for the second moment population averaged kernel.

Using the moment generator (Eq. (13)), we get for the number density averaged kernel $$\overline{K(m,m^{\prime })}=2K_0{\lambda }^{2/3}[\frac{\mathrm{\Gamma }(\nu +2/3)}{\mathrm{\Gamma }(\nu )}+\frac{\mathrm{\Gamma }(\nu +1/3{)}^2}{\mathrm{\Gamma }(\nu {)}^2}],$$(44)

and for the mass-weighted kernel $$\overline{K_m(m,m^{\prime })}=2K_0{\lambda }^{2/3}[\frac{\mathrm{\Gamma }(\nu +5/3)}{\mathrm{\Gamma }(\nu +1)}+\frac{\mathrm{\Gamma }(\nu +4/3{)}^2}{\mathrm{\Gamma }(\nu +1{)}^2}].$$(45)

Using ${\lambda }^{2/3}=m_{\mathrm{c}}^{2/3}{\nu }^{-2/3}$, the representative particle size, r_c, can be returned to the final expression. The final equation for the zeroth moment is then $$\begin{gathered} {\left(\frac{\mathrm{d}{N}_{\mathrm{c}}}{\mathrm{d}t}\right)}_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{l}}\approx -\pi r_{\mathrm{c}}^2\mathrm{\Delta }{v}_{\mathrm{f}}\overline{E}{N}_{\mathrm{c}}^2 \\ {\nu }^{-2/3}[\frac{\mathrm{\Gamma }(\nu +2/3)}{\mathrm{\Gamma }(\nu )}+\frac{\mathrm{\Gamma }(\nu +1/3{)}^2}{\mathrm{\Gamma }(\nu {)}^2}], \end{gathered}$$(46)

and for the second moment is $$\begin{gathered} {\left(\frac{\mathrm{d}{Z}_{\mathrm{c}}}{\mathrm{d}t}\right)}_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{l}}\approx 2\pi r_{\mathrm{c}}^2\mathrm{\Delta }{v}_{\mathrm{f}}\overline{E}{\rho }_{\mathrm{c}}^2 \\ {\nu }^{-2/3}[\frac{\mathrm{\Gamma }(\nu +5/3)}{\mathrm{\Gamma }(\nu +1)}+\frac{\mathrm{\Gamma }(\nu +4/3{)}^2}{\mathrm{\Gamma }(\nu +1{)}^2}]. \end{gathered}$$(47)

Fig. 3 Relative factors between the gamma distribution and monodisperse distribution for the condensation and evaporation rates (left) and collisional growth rates (right). The vertical dotted line shows the ν =1 values, denoting an exponential distribution (Paper I). As ν → ∞, the rates tend towards the monodisperse limit. However, for coagulation in the Kn ≫ 1 regime, assuming H = 0.85 does not recover the monodisperse limit, even for large ν, while a value of H =1/√2 does recover the monodisperse limit.

4.4 Summary of condensation and collisional growth results

In Figure 3, we showed the fractional difference between our derived gamma distribution expressions and the monodisperse values from Lee & Ohno (2025). For the condensation and evaporation rates, we found that the rates are reduced in the gamma distribution formulation by around a factor of two when ν = 0.1 compared to the monodisperse rates. The rate converges towards the monodisperse solution when ν → ∞.

Figure 3 also shows the fractional change in the coagulation and coalescence rates using the gamma distribution compared to the monodisperse formulation in Lee & Ohno (2025). Again, we found the rates converge towards the monodisperse solution when ν → ∞. Generally, coalescence rates are reduced for broader distributions compared to the monodisperse case, and the coagulation rate is enhanced, sometimes by several magnitudes compared to the monodisperse case when the distribution is broad. As stated earlier, for the coagulation rate, this plot shows that assuming the value of H ≈ 0.85 does not recover the monodisperse limit, but the value of H ≈ 1/√2 does. We therefore preferred a value of H ≈ 1/ √2 to ensure the monodisperse limit is recovered. More rigorously, the optimal value of H should depend on the dominant range of masses and should be a function of ν and λ, ranging between 1/√2 and 1.

For completeness, in Appendix C, we also derived population averaged kernels for the turbulence induced shear collisional rate. We did not include the effects of turbulence induced collisions in the current study due to uncertainty in the atmospheric turbulent kinetic energy dissipation rate in sub-stellar atmospheres. However, we note that calculations performed in Samra et al. (2022) suggest turbulence driven collisions may play an important role in setting the overall cloud structure in brown dwarf atmospheres.

5 Y-dwarf KCl cloud simple 1D example

In this section, we apply the three-moment gamma mass distribution framework to a simple 1D KCl Y-dwarf cloud formation example. We used the same setup as Paper I, taking T_eff = 400 K, log g = 3.25 and log g = 4.25 Y-dwarf temperature-pressure (T-p) profiles from Gao et al. (2018) and assuming a constant K_zz = 10⁸ cm² s⁻¹. We included homogeneous nucleation of KCl cloud particles to form seed particles of size 1 nm following the same scheme as in Lee & Ohno (2025). In addition, we updated the monodisperse and exponential distribution models from Paper I, where appropriate, using the methods above in order to create a one-to-one comparison as much as possible. In particular, we have now included moment dependent settling velocities following Appendix B, as well as a more accurate second order implicit Crank-Nicolson diffusion scheme and second order explicit TVD-MUSCL advection scheme. We have plotted the absolute relative difference between the gamma and monodisperse results as $$|\delta N_{\mathrm{c}}|=|\frac{N_{\mathrm{c}}(\mathrm{g}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a})-N_{\mathrm{c}}(\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{o})}{N_{\mathrm{c}}(\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{o})}|,$$(48)

with the equivalent expression for the representative particle size.

It is important to note that numerical issues appear in the formulation for large values of ν due to the inherent scalings present in the equations. First, we converted all equations outlined in this paper to numerically use the ln Γ(z) function rather than the Γ(z) function in the code, which ensures numerical stability for larger ν values. To achieve quick numerical convergence of the lnΓ(z) function, we assumed a maximum of ν = 100 when the moment solutions give a highly monodisperse size distribution. The size distribution is already shows highly monodisperse properties at ν =10 (Fig. 3), suggesting this limit will have a negligible effect on the end results.

Noticeable differences are produced between the gamma and monodisperse results in some cases. Figure 4 shows our results for the 1D column model. In the log g = 3.25 case, the gamma distribution produces smaller particles by a factor of ~2-3 and a larger number density by a factor of ~ 10-20 compared to the monodisperse results. This shows the potentially strong effects of polydispersity on the microphysical rates, as the distribution is quite broad across the deeper pressure levels, but narrows at lower pressure (Figure 5). In the log g = 4.25 case, fewer differences are seen between the gamma and monodisperse assumptions, with both r_c and N_c within around only a relative difference of only ~0.3. The largest differences are seen at the cloud base, where the most polydisperse behavior is present, with the largest divergence between r_c and r_p occurring here.

We observed larger differences in the log g = 3.25 case, which is presumably because the gravitational separation of cloud particles is more effective. As seen in Figure 4, the mass mixing ratio of clouds decreases with pressure in the log g = 3.25 case, indicating that some cloud particles have too fast settling velocity to be transported by eddy diffusion. Meanwhile, we observed subtle differences (<20%) among monodisperse, exponential, and gamma distribution models for the log g = 4.25 case, in which the gravitational separation barely occurs as informed by the vertically constant cloud abundance. Thus, our results suggest that polydipersity of the size distribution does matter when the gravitational separation of cloud particles play a role, as it controls the fraction of small particles that can be transported by eddy diffusion. From Figure 4, the exponential distribution results are shifted towards the gamma distribution solution, compared to the monodisperse solution. This suggests the exponential distribution succeeds in recovering some of the behaviour of the gamma distribution, although the inflexibility of the exponential distribution limits the full realisation of the microphysical effects from polydispersity.

In Figure 5, we showed the reconstructed gamma size distribution, derived from the moment solution values, as a function of pressure in the atmosphere for both test cases. This shows how the size distribution is changing with altitude, with both cases showing a much broader size distribution near the cloud base which narrows to a more peaked distribution at lower pressure. This suggests a differentiation of cloud particle sizes is occurring in the atmosphere, brought on by the relative differences in particle size settling velocity and microphysical effects. Since the second moment settles faster than the first and zeroth moment, unless the distribution is very narrow, the variance of the distribution naturally reduces with pressure. This allows the three-moment model to reproduce the trend of mean particle sizes decreasing with height for the log g = 3.25 case, which is seen in bin based microphysical models (e.g. Gao et al. 2018), which is unable to be captured by the monodisperse model.

Fig. 4 KCl cloud structures using the log g = 3.25 (top row) and log g = 4.25 (bottom row) Y-dwarf temperature-pressure profiles from Gao et al. (2018) and assuming a constant Kzz = 108 cm2 s−1. The dashed lines show the monodisperse size distribution, the dotted lines the exponential distribution, and the solid lines the exponential size distribution results. The left panels show the temperature-pressure profiles with the mass mixing ratio of the condensate, q1, vapour, qv and saturation point, qs (pink dash-dot line). The middle panels show the representative particle sizes, rc and rp, and total number density, Nc. The right panels show the relative difference in rc and Nc between the monodisperse and gamma cloud structure.

Fig. 5 Reconstructed gamma particle size distributions from the simulations performed in Section 5 for the log g = 3.25 (left) and log g = 4.25 (right) cases. This shows the size distribution properties with pressure (colour bar) in the atmosphere. In general, the size distribution narrows with decreasing pressure, with the most broad distributions occurring near the cloud base.

6 Discussion

Our results suggest that strong polydispersity is present at the cloud base, but as the cloud structure reaches lower pressures the cloud size distribution narrows, eventually showing highly monodisperse characteristics in the upper atmosphere. We suggest this is a key benefit from the three-moment scheme, where the variance of the distribution is self-consistently calculated due to the addition of the second moment to the scheme. In addition, we suggest the variation in settling velocity between the moments, where the second moment settles at the highest rate, is an important part of this processes, as reducing the value of Z_c relative to ρ_c and N_c decreases the variance of the distribution with height. This is seen most strongly in Figure 4 for the log g = 3.25 case, where a gradient towards smaller particles is produced in the representative particle size in the upper atmosphere for the gamma distribution, while the monodisperse results show a constant value.

Woitke et al. (2020) outline an approach to limit the coupling of cloud particles to the diffusive mixing of the atmosphere dependent on the local Stokes number of the particle relative to the largest turbulent eddy size. In this study, we assumed all particles were equally affected by the vertical diffusive motions through the K_zz value. We suggest that this Stokes number would also be moment dependent, resulting in further differentiation between the moment solutions, due to the first and second moments being more decoupled from the diffusive vertical mixing rate than the zeroth. This would likely result in a further decrease in the variance and representative particle mass in the upper atmosphere compared to the results presented in this study.

Compared to Paper I, our gamma distribution results show better agreement with the microphysical bin-resolving model CARMA performed in Gao et al. (2018). For the log g = 3.25 case, the mixing ratio profile shows excellent agreement with the CARMA, as well as the representative particle size for the gamma distribution case now being smaller than the exponential and monodisperse case, more in line with the CARMA results. However, for the log g = 4.25 case, CARMA recovers particle sizes in the ~3 micron range, slightly larger than that produced here. This may suggest strong multi-model behavior in this case, with the moment method capturing the smaller particle size population while CARMA contains a secondary larger particle sized mode. We note that due to differences in the microphysical equations used, direct one-to-one comparison between our results and CARMA should be treated with caution. However, due to the similarly in overall cloud structures between the bin and bulk models, this suggests the bulk model is capturing well the complex behaviour of the bin model. Despite this, we note that many bin models applied to exoplanet atmosphere clouds show strong evidence of bi- or multi-modal size distributions (e.g. Powell et al. 2019), typically with a nucleation secondary component composed of small particles and a primary component of particles that have undergone condensational and collisional growth. Possibly, to capture this effect using the moment method, multiple moments can be used to describe different parts of the distribution, such as the ‘cloud’ and ‘rain’ components used in Ohno & Okuzumi (2017). We note that such multi-moment schemes are standard practice for Earth water cloud microphysical models (e.g. Seifert & Beheng 2001).

Our results from Figure 5 suggest that the particle size distribution is likely to be well differentiated between the cloud base and cloud top, with a broad distribution at the cloud base and more peaked particle sizes at the cloud top. We suggest that this behaviour fits well with current retrieval results for clouds in brown dwarf atmospheres (e.g. Burningham et al. 2021), which best fit a broad particle size at the cloud base with a small particle component higher in the upper atmosphere. With a broader particle size and increased number density, the cloud opacity is generally larger at the cloud base which decreases at the distribution becomes thinner and more monodisperse with a smaller particle size. A broader distribution also reduces the strength of the infrared absorption features present in hot Jupiter and brown dwarf data (e.g. Burningham et al. 2021; Grant et al. 2023). This is generally due to the contribution of larger particles in the distribution as shown in Wakeford & Sing (2015). Should a more peaked size distribution be present in the observable upper atmosphere, emergent from a deeper atmosphere, broader distribution at the cloud base, this would enhance size specific cloud signatures at infrared wavelengths, as well as help produce a size-dependent scattering slope at optical wavelengths.

7 Conclusions

In this follow-up study to Lee & Ohno (2025) and Lee (2025), we presented a self-consistent three-moment scheme for sub-stellar cloud microphysics including condensation, nucleation and collisional growth assuming a gamma size distribution. Importantly, the present three-moment scheme is capable of simulating the evolution of the size distribution shape, which was fixed to the more inflexible monodisperse or exponential distributions throughout atmospheres in the two-moment scheme. Overall, our derived equation set results in a modification to the monodisperse rates with the variance of the distribution also taken into account, going beyond the constant variance of the exponential distribution. As a check, we showed that our equation set trends towards the monodisperse solution when very narrow size distribution parameters are produced. These modifications come in the form of simple additional gamma function terms and parameters that depend on the size distribution properties. Here, we aimed to directly compare the monodisperse and polydisperse assumptions and isolate their effects on the cloud condensation, evaporation, and collisional growth rates.

In our 1D KCl Y-dwarf cloud test, we found differences between the monodisperse and gamma distribution results, stemming from the effects of polydispersity altering condensation, Brownian coagulation and gravitational coalescence rates compared to the monodisperse rates. Generally, Brownian coagulation rates are enhanced with polydispersity, while the gravitational coalescence rates are reduced. We showed that the variance of the distribution changes with height, with the size distribution narrowing with altitude. The log g = 3.25 case shows noticeable differences between the three-moment scheme and previous two-moment scheme with monodisperse and exponential distribution, suggesting that polydispersity plays a more important role for low gravity objects where the gravitational separation of cloud particles is more pronounced.

With consideration of additional first moments for different condensation species, mixed material grains can be taken into account similar to those in Helling et al. (2008), which is left to future studies. Other three moment schemes utilising other size distribution assumptions, such as the lognormal distribution, can also be easily derived using our framework.

Acknowledgements

E.K.H.L is supported by the CSH through the Bernoulli Fellowship. K.O. is supported by the JSPS KAKENHI Grant Number JP23K19072.

References

Appendix A Cause of gamma function divergence issue

In Section 4, we show that the coagulation rate of Eqs. (36) and (41) derived from the gamma distribution includes a gamma function that potentially involves a zero or negative argument. This seems problematic, as the coagulation rate diverges to infinity at ν → 1/3 and ν → 2/3 in Eq. (36), which is clearly unrealistic. This odd behavior arises because we have ignored the fact that the particle size distribution actually has a cut off at a certain small particle size. For example, cloud particles cannot be smaller than the size of constituting molecules. In practice, enhanced evaporation due to particle curvature (Kelvin effect) would prohibit the presence such small particles.

To make the problem clearer, let us consider a thought experiment: a huge single particle is embedded in swarm of cloud particles. Then, the rate of decrease in total number density of cloud particles due to collision between a single particle with a mass of M and cloud particles. For the Brownian coagulation in the continuum regime (Kn ≫ 1), one can describe the rate as $$\begin{array}{r}\frac{d{N}_{\mathrm{c}}}{dt}=-{\int }_{{\mathrm{m}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{{\mathrm{m}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}K(m,M)f(m)dm\\ =-\frac{K_0{N}_{\mathrm{c}}{M}^{1/3}}{{\lambda }^{\nu }\mathrm{\Gamma }(\nu )}{\int }_{{\mathrm{m}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{{\mathrm{m}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{m}^{\nu -4/3}{e}^{-m/\lambda }dm,\end{array}$$(A.1)

where the particle size distribution is postulated to have maximum and minimum masses of m_max [g] and m_min [g], and we have used Eqs. (1) with M≫m_max. We can see that the contribution from the collisions of particles with masses of m is then $$\mathrm{\Delta }{N}_{\mathrm{c}}(m)=m^{\nu -4/3}{e}^{-m/\lambda }\mathrm{\Delta }m\sim m^{\nu -1/3}{e}^{-m/\lambda }.$$(A.2)

In this equation, the contribution from smallest particles is completely different between the cases of ν > 1/3 and ν < 1/3. ∆N_c decreases with decreasing particle mass at ν > 1/3, implying that smallest particles barely contribute to the total collision rate, which justifies m_min = 0 for ν > 1/3. However, for ν < 1/3, ∆N_c turns out to increase with decreasing particle mass, implying that the collisions of smallest particles with m = m_min have the largest contribution to the total collision rate. In this case, one cannot set m_min = 0, as the collision rate does not converge to zero and diverges to infinity at smallest particle mass. Physically, this behavior happens because collision velocity is faster at smaller particles, and the size distribution needs to drop faster than the rise of collision velocity to make the collision rate being zero at m = 0. Note that the exponential cutoff $e^{-m/\lambda }$ guarantees that the collision rate always converges to zero at large size limit, justifying to set m_max = ∞ for any ν.

A.1 Prescription to correct gamma function negative arguments

To fix the divergence issue, one can introduce a prescription that accounts for the finite value of the minimum particle mass m_min. Since divergence issues are present for the population-averaged kernel K(m, m′) due to the reasons explained above, we need only to apply the prescription to terms in K(m, m′) that diverge. Supposing that the collisions of particles with m < m_min do not occur in reality, we can rewrite the population averaged collision kernel as $$\overline{K(m,m^{\prime })}=\frac{1}{N_{\mathrm{c}}^2}{\int }_{m_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{\mathrm{\infty }}{\int }_{m_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{\mathrm{\infty }}K(m,m^{\prime })f(m)f(m^{\prime })dmd{m}^{\prime },$$(A.3)

For the later convenience, we define the function $$I^{(k)}={\int }_{{\mathrm{m}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{\mathrm{\infty }}{m}^{k}f(m)dm=N_{\mathrm{c}}{\lambda }^k\frac{\mathrm{\Gamma }(\nu +k,x_{\mathrm{m}\mathrm{i}\mathrm{n}})}{\mathrm{\Gamma }(\nu )},$$(A.4)

where x_min = m_min/d is the dimensionless minimum mass, and $\mathrm{\Gamma }(s,x)\equiv {\int }_{\mathrm{x}}^{\mathrm{\infty }}{t}^{s-1}{e}^{-t}dt$ is the upper incomplete gamma function. This function returns the k-th moment from Eq. (13) as x_min → 0. In our current models, x_min is taken to be the mass of particles at the seed particle size of 1 nm.

Here, we repeat the derivation of the averaged kernel with the above prescription. For the Brownian coagulation in the particle diffusive regime (Kn ≪ 1), the averaged kernel can be rewritten as $$\begin{gathered} \overline{K(m,m^{\prime })}=\frac{K_0}{N_{\mathrm{c}}^2}{\int }_{m_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{\mathrm{\infty }}{\int }_{m_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{\mathrm{\infty }}[\beta m^{-1/3}+{\beta }^{\prime }{m}^{\prime -1/3}](m^{1/3}+m^{\prime 1/3}) \\ \cdot f(m)f(m^{\prime })dmd{m}^{\prime }. \end{gathered}$$(A.5)

Adopting the linear slip factor of β = 1 + AKn, the population-averaged kernel is given by $$\begin{array}{r}\overline{K(m,m^{\prime })}=\frac{2K_0}{N_{\mathrm{c}}^2}[I^{(0)}{I}^{(0)}+I^{(1/3)}{I}^{(-1/3)}\\ +{\lambda }_{\mathrm{a}}A{\left(\frac{3}{4\pi {\rho }_{\mathrm{d}}}\right)}^{-1/3}(I^{(0)}{I}^{(-1/3)}+I^{(1/3)}{I}^{(-2/3)})].\end{array}$$(A.6)

Applying the generator function and inserting the population-averaged kernel, the rate of change in number density is given by $$\begin{array}{r}{\left(\frac{d{N}_{\mathrm{c}}}{dt}\right)}_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{g}}\approx -\frac{2k_{\mathrm{b}}T}{3{\eta }_{\mathrm{a}}}{N}_{\mathrm{c}}^2[\frac{\mathrm{\Gamma }(\nu ,x_{\mathrm{m}\mathrm{i}\mathrm{n}}{)}^2}{\mathrm{\Gamma }(\nu {)}^2}+\frac{\mathrm{\Gamma }(\nu +1/3,x_{\mathrm{m}\mathrm{i}\mathrm{n}})\mathrm{\Gamma }(\nu -1/3,x_{\mathrm{m}\mathrm{i}\mathrm{n}})}{\mathrm{\Gamma }(\nu {)}^2}\\ +\frac{{\lambda }_{\mathrm{a}}}{r_{\mathrm{c}}}A{\nu }^{1/3}(\frac{\mathrm{\Gamma }(\nu ,x_{\mathrm{m}\mathrm{i}\mathrm{n}})\mathrm{\Gamma }(\nu -1/3,x_{\mathrm{m}\mathrm{i}\mathrm{n}})}{\mathrm{\Gamma }(\nu {)}^2}\\ +\frac{\mathrm{\Gamma }(\nu +1/3,x_{\mathrm{m}\mathrm{i}\mathrm{n}})\mathrm{\Gamma }(\nu -2/3,x_{\mathrm{m}\mathrm{i}\mathrm{n}})}{\mathrm{\Gamma }(\nu {)}^2})].\end{array}$$(A.7)

For the kinetic regime of Brownian motion, the population-averaged kernel can be rewritten as $$\begin{array}{r}\overline{K(m,m^{\prime })}=\frac{K_{0}H}{N_{\mathrm{c}}^2}{\int }_{m_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{\mathrm{\infty }}{\int }_{m_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{\mathrm{\infty }}{(m^{1/3}+m^{\prime 1/3})}^2(m^{-1/2}+m^{\prime -1/2})\\ \cdot f(m)f(m^{\prime })dmd{m}^{\prime }.\end{array}$$(A.8)

This equation yields $$\begin{array}{r}\overline{K(m,m^{\prime })}=\frac{2K_{0}H}{N_{\mathrm{c}}^2}(I^{(2/3)}{I}^{(-1/2)}+2I^{(1/3)}{I}^{(-1/6)}+I^{(1/6)}{I}^{(0)}).\end{array}$$(A.9)

Again, using the generator function of Eq. (A.4) and inserting the population-averaged kernel, the rate of change in number density is given by $$\begin{array}{r}{\left(\frac{d{N}_{\mathrm{c}}}{dt}\right)}_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{g}}\approx -H\sqrt{\frac{8\pi k_{\mathrm{b}}T}{m_{\mathrm{c}}}}{r}_{\mathrm{c}}^2{N}_{\mathrm{c}}^2{\nu }^{-1/6}[\frac{\mathrm{\Gamma }(\nu +2/3,x_{\mathrm{m}\mathrm{i}\mathrm{n}})\mathrm{\Gamma }(\nu -1/2,x_{\mathrm{m}\mathrm{i}\mathrm{n}})}{\mathrm{\Gamma }(\nu {)}^2}\\ +\frac{2\mathrm{\Gamma }(\nu +1/3,x_{\mathrm{m}\mathrm{i}\mathrm{n}})\mathrm{\Gamma }(\nu -1/6,x_{\mathrm{m}\mathrm{i}\mathrm{n}})}{\mathrm{\Gamma }(\nu {)}^2}+\frac{\mathrm{\Gamma }(\nu ,x_{\mathrm{m}\mathrm{i}\mathrm{n}})\mathrm{\Gamma }(\nu +1/6,x_{\mathrm{m}\mathrm{i}\mathrm{n}})}{\mathrm{\Gamma }(\nu {)}^2}].\end{array}$$(A.10)

Due to the reasons outlined above, in the model we only need to apply the incomplete gamma function for terms that may produce a negative or zero ν argument. For positive arguments to the gamma function, the incomplete gamma distribution is not required to be used and the equations outlined in the main text can be applied. This is shown in Figure 2, where, when a positive argument is given for both the incomplete gamma function and gamma function, the results are near identical. Therefore there is no significant loss of accuracy from using the gamma distribution for a positive ν argument than the incomplete gamma function.

For the number-weighted Knudsen number, Kn_N, the divergence issues are also solved by applying the incomplete gamma function rather than the gamma function, leading to $$\overline{{\mathrm{K}\mathrm{n}}_N}=\frac{{\lambda }_{\mathrm{a}}}{r_{\mathrm{c}}}{\nu }^{1/3}\frac{\mathrm{\Gamma }(\nu -1/3,x_{\mathrm{m}\mathrm{i}\mathrm{n}})}{\mathrm{\Gamma }(\nu )}.$$(A.11)

As shown in Figure 2, this avoids the singularities present in the gamma distribution formulation.

At first glance, it may seem that using the incomplete gamma function does not solve the negative argument to the gamma function problem as the possibility still exists in Eqs. (A.7) and (A.10). However, we can utilise the recurrence relation for the incomplete gamma function $$\mathrm{\Gamma }(\nu ,x_{\mathrm{m}\mathrm{i}\mathrm{n}})=\frac{\mathrm{\Gamma }(\nu +1,x_{\mathrm{m}\mathrm{i}\mathrm{n}})-x_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\nu }{e}^{-x_{\mathrm{m}\mathrm{i}\mathrm{n}}}}{\nu },$$(A.12)

which ensures that a positive argument is always present when computing the incomplete gamma function. We note that using the gamma function recurrence relation directly $$\mathrm{\Gamma }(\nu )=\frac{\mathrm{\Gamma }(\nu +1)}{\nu },$$(A.13)

would not solve this issue, as the exponential term in the incomplete gamma recurrence relation is required to properly bound the solution. We use the numerical implementation of the upper incomplete gamma function from Press et al. (1986), which we find is computationally stable and efficient enough for the microphysics code.

Appendix B Distribution weighted settling velocity

The population averaged settling velocity for a moment k can be given by (e.g. Straka 2009) $$\overline{v_{\mathrm{f}}^{(k)}}=\frac{{\int }_0^{\mathrm{\infty }}{v}_{\mathrm{f}}(m)m^{k}f(m)dm}{{\int }_0^{\mathrm{\infty }}{m}^{k}f(m)dm}=\frac{1}{M^{(k)}}{\int }_0^{\mathrm{\infty }}{v}_{\mathrm{f}}(m)m^{k}f(m)dm,$$(B.1)

where v_f(m) [cm s⁻¹] is the settling velocity of a particle with mass m. For the Stokes regime (Kn ≪ 1), we take the v_f(r) expression from Ohno & Okuzumi (2017, 2018) $$v_{\mathrm{f}}(r)=\frac{2\beta g{r}^2({\rho }_{\mathrm{d}}-{\rho }_{\mathrm{a}})}{9{\eta }_{\mathrm{a}}}{[1+{\left(\frac{0.45g{r}^3{\rho }_{\mathrm{a}}{\rho }_{\mathrm{d}}}{54{\eta }_{\mathrm{a}}^2}\right)}^{2/5}]}^{-5/4}.$$(B.2)

Ignoring the Cunningham slip factor and Reynolds number correction term, the population averaged velocity expression for the Stokes regime is $$\begin{array}{r}\overline{v_{\mathrm{f},\mathrm{s}\mathrm{t}\mathrm{k}}^{(k)}}=\frac{qs}{M^{(k)}}{\int }_0^{\mathrm{\infty }}{m}^{k+2/3}f(m)dm=\\ qs\frac{M^{(k+2/3)}}{M^{(k)}}=qs{\lambda }^{k+2/3}\frac{\mathrm{\Gamma }(\nu +k+2/3)}{\mathrm{\Gamma }(\nu +k)},\end{array}$$(B.3)

where $q=2g({\rho }_{\mathrm{d}}-{\rho }_{\mathrm{a}})/9{\eta }_{\mathrm{a}}$ and $s=(3/4\pi {\rho }_{\mathrm{d}}{)}^{2/3}$. This results in the k = 0, 1 and 2 population averaged settling velocities being $$\overline{v_{\mathrm{f},\mathrm{s}\mathrm{t}\mathrm{k}}^{(0)}}=\frac{2g({\rho }_{\mathrm{d}}-{\rho }_{\mathrm{a}})}{9{\eta }_{\mathrm{a}}}{r}_{\mathrm{c}}^2{\nu }^{-2/3}\frac{\mathrm{\Gamma }(\nu +2/3)}{\mathrm{\Gamma }\left(\nu \right)},$$(B.4) $$\overline{v_{\mathrm{f},\mathrm{s}\mathrm{t}\mathrm{k}}^{(1)}}=\frac{2g({\rho }_{\mathrm{d}}-{\rho }_{\mathrm{a}})}{9{\eta }_{\mathrm{a}}}{r}_{\mathrm{c}}^2{\nu }^{-2/3}\frac{\mathrm{\Gamma }(\nu +5/3)}{\mathrm{\Gamma }(\nu +1)},$$(B.5)

and $$\overline{v_{\mathrm{f},\mathrm{s}\mathrm{t}\mathrm{k}}^{(2)}}=\frac{2g({\rho }_{\mathrm{d}}-{\rho }_{\mathrm{a}})}{9{\eta }_{\mathrm{a}}}{r}_{\mathrm{c}}^2{\nu }^{-2/3}\frac{\mathrm{\Gamma }(\nu +8/3)}{\mathrm{\Gamma }(\nu +2)},$$(B.6)

respectively. Assuming a linear form of β, the above equations are augmented with a moment dependent slip correction term, leading to $$\overline{v_{\mathrm{f}}^{(0)}}=\overline{v_{\mathrm{f},\mathrm{s}\mathrm{t}\mathrm{k}}^{(0)}}[1+A\frac{{\lambda }_{\mathrm{a}}}{r_{\mathrm{c}}}{\nu }^{1/3}\frac{\mathrm{\Gamma }(\nu +1/3)}{\mathrm{\Gamma }(\nu +2/3)}],$$(B.7) $$\overline{v_{\mathrm{f}}^{(1)}}=\overline{v_{\mathrm{f},\mathrm{s}\mathrm{t}\mathrm{k}}^{(1)}}[1+A\frac{{\lambda }_{\mathrm{a}}}{r_{\mathrm{c}}}{\nu }^{1/3}\frac{\mathrm{\Gamma }(\nu +4/3)}{\mathrm{\Gamma }(\nu +5/3)}],$$(B.8)

and $$\overline{v_{\mathrm{f}}^{(2)}}=\overline{v_{\mathrm{f},\mathrm{s}\mathrm{t}\mathrm{k}}^{(2)}}[1+A\frac{{\lambda }_{\mathrm{a}}}{r_{\mathrm{c}}}{\nu }^{1/3}\frac{\mathrm{\Gamma }(\nu +7/3)}{\mathrm{\Gamma }(\nu +8/3)}],$$(B.9)

respectively.

Fig. B.1 Relative difference between the settling velocity of the gamma distribution and the monodisperse distribution for the zeroth moment (solid lines), first moment (dashed lines) and second moment (dash-dot lines) for the Kn ≪ 1 regime (blue colour) and Kn ≫ 1 regime (orange colour).

Lastly, including the Reynolds number correction factor in a fully analytical solution is more difficult due to the non-linearity of the equation. We therefore estimate the Reynolds number factor in Eq. (B.2) as a constant, given at the respective population weighted particle size for each moments settling velocity. We can use the gamma distribution parameters to estimate the respective population weighted mean radii. The number-weighted mean radius, (r) [cm], is estimated from $$⟨r⟩=r_{\mathrm{c}}{\nu }^{-1/3}\frac{\mathrm{\Gamma }(\nu +1/3)}{\mathrm{\Gamma }\left(\nu \right)},$$(B.10)

the number-weighted mass representative radius, r_c [cm], given by Eq. (21), and the mass-weighted mass representative radius, r_p [cm], given by Eq. (23).

The β slip factor used in the Stokes regime is only valid between values of Kn ≈ 0.5-83 (Kim et al. 2005), with no guarantee that a universal asymptotic value is reached in the Kn ≫ 1 limit. In addition, we require to use a linear form of β in our formulations, which, due to absence of the exponential factor, does not coverage to a constant solution for large Kn. We therefore include a solution for particles in the Epstein (Kn ≫ 1) regime. We use the settling velocity equation from Woitke & Helling (2003) $$v_{\mathrm{f}}(r)=\frac{\sqrt{\pi }g{\rho }_{\mathrm{d}}r}{2{\rho }_{\mathrm{a}}{c}_{\mathrm{T}}},$$(B.11)

where c_T [cm s⁻¹] is the atmospheric thermal velocity (Woitke & Helling 2003). Defining $q=\sqrt{\pi }g{\rho }_{\mathrm{d}}/2{\rho }_{\mathrm{a}}{c}_{\mathrm{T}}$ and $s=(3/4\pi {\rho }_{\mathrm{d}}{)}^{1/3}$, the population averaged velocity expression is $$\begin{array}{r}\overline{v_{\mathrm{f}}^{(k)}}=\frac{qs}{M^{(k)}}{\int }_0^{\mathrm{\infty }}{m}^{k+1/3}f(m)dm=\frac{qs{M}^{(k+1/3)}}{M^{(k)}}\\ =qs{\lambda }^{k+1/3}\frac{\mathrm{\Gamma }(\nu +k+1/3)}{\mathrm{\Gamma }(\nu +k)}.\end{array}$$(B.12)

This results in the k = 0, 1 and 2 population averaged settling velocities being $$\overline{v_{\mathrm{f}}^{(0)}}=\frac{\sqrt{\pi }g{\rho }_{\mathrm{d}}}{2{\rho }_{\mathrm{a}}{c}_{\mathrm{T}}}{r}_{\mathrm{c}}{\nu }^{-1/3}\frac{\mathrm{\Gamma }(\nu +1/3)}{\mathrm{\Gamma }\left(\nu \right)},$$(B.13) $$\overline{v_{\mathrm{f}}^{(1)}}=\frac{\sqrt{\pi }g{\rho }_{\mathrm{d}}}{2{\rho }_{\mathrm{a}}{c}_{\mathrm{T}}}{r}_{\mathrm{c}}{\nu }^{-1/3}\frac{\mathrm{\Gamma }(\nu +4/3)}{\mathrm{\Gamma }(\nu +1)},$$(B.14)

and $$\overline{v_{\mathrm{f}}^{(2)}}=\frac{\sqrt{\pi }g{\rho }_{\mathrm{d}}}{2{\rho }_{\mathrm{a}}{c}_{\mathrm{T}}}{r}_{\mathrm{c}}{\nu }^{-1/3}\frac{\mathrm{\Gamma }(\nu +7/3)}{\mathrm{\Gamma }(\nu +2)},$$(B.15)

respectively.

In Figure B.1, we show the fractional difference between the gamma distribution population averaged settling velocity and the equivalent monodisperse settling velocity for both the Stokes and Epstein regimes. From this, we see that the moments acquire different settling velocities to each other, with the zeroth moment settling slower than the first and second, and the second settling the fastest. This is expected, as each moment represents a different integrated value and characteristic of the distribution which would naturally have different settling velocities (e.g. Woitke et al. 2020). For large ν, the distribution becomes more monodisperse, resulting in all the moment dependent settling velocities converging on the monodisperse value. For intermediate Knudsen number regimes (Kn ~ 1) the same tanh function can be used from Paper I for each moment dependent settling velocity to interpolate between the Stokes and Epstein expressions.

Appendix C Turbulent shear kernel

The turbulence driven wind shear kernel is given by (Saffman & Turner 1956) $$K(r,r^{\prime })={\left(\frac{8\pi {ϵ}_{\mathrm{d}}}{15{\nu }_{\mathrm{a}}}\right)}^{1/2}(r+r^{\prime }{)}^3,$$(C.1)

where e_d [cm² s⁻³] is the dissipation of turbulent kinetic energy rate of the atmosphere and ν_a [cm² s⁻¹] the atmospheric kinematic viscosity. Using the Morán (2022) population averaged kernel method, the resulting population averaged kernels for the zeroth and second moment are $$\overline{K(m,m^{\prime })}=\frac{2K_0}{N_{\mathrm{c}}^2}(M^{(0)}{M}^{(1)}+3M^{(2/3)}{M}^{(1/3)}),$$(C.1)

and $$\overline{K_m(m,m^{\prime })}=\frac{2K_0}{{\rho }_{\mathrm{c}}^2}(M^{(1)}{M}^{(2)}+3M^{(5/3)}{M}^{(4/3)}),$$(C.3)

respectively, where K₀ is the kinetic pre-factor $$K_0={\left(\frac{8\pi {ϵ}_{\mathrm{d}}}{15{\nu }_{\mathrm{a}}}\right)}^{1/2}\left(\frac{3}{4\pi {\rho }_{\mathrm{d}}}\right).$$(C.4)

For the gamma distribution, the number density averaged kernel is then $$\overline{K(m,m^{\prime })}=2K_0\lambda [\frac{\mathrm{\Gamma }(\nu +1)}{\mathrm{\Gamma }(\nu )}+\frac{3\mathrm{\Gamma }(\nu +2/3)\mathrm{\Gamma }(\nu +1/3)}{\mathrm{\Gamma }(\nu {)}^2}],$$(C.5)

and the mass-weighted kernel $$\overline{K(m,m^{\prime })}=2K_0\lambda [\frac{\mathrm{\Gamma }(\nu +2)}{\mathrm{\Gamma }(\nu +1)}+\frac{3\mathrm{\Gamma }(\nu +5/3)\mathrm{\Gamma }(\nu +4/3)}{\mathrm{\Gamma }(\nu +1{)}^2}],$$(C.6)

The expression for the zeroth moment is then $$\begin{array}{r}{\left(\frac{d{N}_{\mathrm{c}}}{dt}\right)}_{\mathrm{T}\mathrm{S}}\approx -{\left(\frac{8\pi {ϵ}_{\mathrm{d}}}{15{\nu }_{\mathrm{a}}}\right)}^{1/2}{r}_{\mathrm{c}}^3{N}_{\mathrm{c}}^2{\nu }^{-1}\\ \cdot [\frac{\mathrm{\Gamma }(\nu +1)}{\mathrm{\Gamma }(\nu )}+\frac{3\mathrm{\Gamma }(\nu +2/3)\mathrm{\Gamma }(\nu +1/3)}{\mathrm{\Gamma }(\nu {)}^2}],\end{array}$$(C.7)

and for the second moment is $$\begin{array}{r}{\left(\frac{d{Z}_{\mathrm{c}}}{dt}\right)}_{\mathrm{T}\mathrm{S}}\approx 2{\left(\frac{8\pi {ϵ}_{\mathrm{d}}}{15{\nu }_{\mathrm{a}}}\right)}^{1/2}{r}_{\mathrm{c}}^3{\rho }_{\mathrm{c}}^2{\nu }^{-1}\\ \cdot [\frac{\mathrm{\Gamma }(\nu +2)}{\mathrm{\Gamma }(\nu +1)}+\frac{3\mathrm{\Gamma }(\nu +5/3)\mathrm{\Gamma }(\nu +4/3)}{\mathrm{\Gamma }(\nu +1{)}^2}].\end{array}$$(C.8)