© The Authors 2025

1 Introduction

The cloud particle size distribution plays a pivotal role in setting the observational properties of atmospheres that contain a cloud component, as the total opacity of the cloud is weighted by the individual contribution of each particle size to the size distribution. Further knock-on effects affect the thermal structure of the atmosphere through feedback with the radiative-transfer properties of the atmospheres, which, for example, in the case of brown dwarf atmospheres, can trigger local convection (e.g. Tan & Showman 2019), leading to changes in local chemical compositions (e.g. Lee et al. 2024). This makes the consideration and accurate modelling of cloud properties a key component in contemporary exoplanet and brown dwarf atmospheric modelling.

Exponential size distributions have been used for modelling Earth-based clouds, particularly raindrops, for decades, since Marshall & Palmer (1948) proposed an exponential distribution as a best fit to observed raindrop data. Exponential distributions have consistently been applied to model raindrops and snow in the Earth atmosphere literature (e.g. Pruppacher & Klett 1978; Lin et al. 1983; Hong & Lim 2006; Seifert 2008; Morrison et al. 2020), with these schemes being used in general circulation models (GCMs) and other atmosphere simulation efforts such as mesoscale and large eddy simulations. However, for exoplanet and brown dwarf atmospheres the exact particle size distribution is unknown and is currently not directly measurable with sufficient accuracy. Therefore, instead of parametrisation and fitting from in situ data, as is typically performed for Earth water clouds, it is advisable to develop theories and models in a manner that is generally applicable and self-consistent within their specific assumptions.

The moment (or bulk) method evolves integrated quantities of the particle size distribution over time, considering how cloud microphysical processes such as nucleation, condensation, evaporation, and collisions affect these quantities. The moment method is in contrast to the bin or spectral method, which directly integrates a discretised grid of particle sizes or masses. A well-used example of a bin model for sub-stellar objects is the 1D CARMA model (e.g. Gao et al. 2018; Powell et al. 2019). However, the bin method is typically too computationally expensive to couple to large-scale hydrodynamic models such as GCMs with contemporary resources, making the development of moment schemes a useful middle ground between computational efficiency and understanding the spatial time evolution of cloud microphysics.

Examples of moment methods used in the sub-stellar atmosphere literature context range from 1D static approaches (e.g. Woitke & Helling 2004; Helling et al. 2008) and 1D time-dependent approaches (e.g. Ohno & Okuzumi 2017; Woitke et al. 2020) to 3D time-dependent methods (e.g. Lee et al. 2016; Lines et al. 2018; Lee 2023). Several studies on sub-stellar atmospheres have employed multiple moment methods derived for monodisperse size distributions and subsequently reconstructed the particle size distribution using the moment solutions (e.g. Helling et al. 2008; Stark et al. 2015; Lee 2023).

Previous efforts typically assumed a delta function for a monodisperse size distribution to derive a mass moment equation set for the cloud microphysical system. In this short study, we go beyond a monodisperse size distribution and explore how the assumption of an exponential particle size distribution affects the condensation, evaporation and collisional growth rates of cloud particles for the mass moment method. In Section 2, we introduce the exponential distribution and its useful properties when applied to the moment method framework. In Section 3, we derive modifications to the monodisperse size distribution mass condensation and evaporation rates for an exponential distribution. We also propose a tanh function that smoothly interpolates between the continuum (diffusive) and free molecular (kinetic) regime limits near the critical Knudsen number. In Section 4, we derive expressions for collisional growth for Brownian coagulation and gravitational coalescence. In Section 5, we propose an interpolation function for the settling velocity between the continuum and free molecular regime limits. In Section 6, we apply our new theory to simple 1D models for a Y-dwarf KCl cloud formation case. Lastly, Sections 7 and 8, respectively, contain the discussion and conclusion of our study. The modelling and code from this paper can be found as part of the ‘mini-cloud’ framework in the GitHub repository¹.

2 Moments with an assumed exponential distribution

In this section, we present a derivation of the mass moments assuming an exponential distribution, and derive the source terms required for evolving the cloud formation microphysics ordinary differential equation (ODE) system in time. The exponential distribution, f(m) [cm⁻³ g⁻¹], is a special single parameter form of the gamma distribution given by $$f(m)=\frac{N_{\mathrm{c}}}{\lambda }{e}^{-m/\lambda },$$(1)

where N_c [cm⁻³] is the total number density of the distribution, λ [g] is the mean particle mass, which is equal to the expectation value, E[m], of the distribution, and m [g] is an individual mass in the distribution. The exponential distribution is a favourable distribution when attempting to derive a set of moment equations for a cloud microphysics system due its single parameter and useful analytical properties.

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)dm,$$(2)

where k is the moment power. Inserting the exponential distribution in the moment equation gives $$M^{(k)}=\frac{N_{\mathrm{c}}}{\lambda }{\int }_0^{\mathrm{\infty }}{m}^k{e}^{-m/\lambda }dm.$$(3)

After some algebra and applying the definition of the gamma integral², the moment generator for the exponential size distribution is analytically expressed as $$M^{(k)}=N_{\mathrm{c}}{\lambda }^k\mathrm{\Gamma }(k+1),$$(4)

where Γ is the gamma function. The moment generator is useful to close the moment system when required through estimating the higher moment values from lower ones; for example, when calculating the moment dependent settling velocity (e.g. Woitke et al. 2020). 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}},$$(5)

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

The value of λ can be readily estimated from the moment solutions as $$\lambda =\frac{M^{(1)}}{M^{(0)}}=m_{\mathrm{c}},$$(7)

which closes the system of equations. Here, following the notation of Lee & Ohno (2025), m_c [g] is the average mass of the particles. The mean mass-weighted average particle size, r_c [cm], is then $$r_{\mathrm{c}}={\left(\frac{3m_{\mathrm{c}}}{4\pi {\rho }_{\mathrm{d}}}\right)}^{1/3},$$(8)

where ρ_d [g cm⁻³] is the particle bulk density. Throughout this study, we assume spherical particles of a homogeneous composition using Eq. (8). We do not consider non-spherical or porous particles.

In Fig. 1, we show visually the difference between the monodisperse and exponential distributions. From this, it is clear that the monodisperse distribution through its assumption of a single representative grain size places all emphasis on the mean particle size, whereas the exponential distribution is more ‘realistic’ in that there is a variance around the mean particle size. In addition, from Fig. 1 the particle mass distribution is generally spread across a large range of values compared to the particle size distribution. In the exponential distribution, this variance is constant and simply given by the square of the average value, $$\mathrm{Var}[m]={\lambda }^2.$$(9)

This single parameter property makes the exponential distribution useful for two-moment methods where only the mean particle mass or size can be estimated from the moment solutions. The exponential distribution is characterised by a gradual increase up to the mean particle size, after which the particle number density drops quickly. This places greater emphasis on the small particle population than the monodisperse distribution, a factor that must be considered when modelling physical processes affecting the size distribution. However, this single parameter approach has significant limitations on the flexibility of the size distribution (Section 7).

Fig. 1 Visualisation of the exponential size distribution compared against a monodisperse size distribution in size units (left) and mass units (right). This shows exponential distributions with (solid lines) λ = 10−3 (blue), 10−2 (orange), 0.1 (green), and 1 (red) μm, assuming Nc = 1 cm−3 and ρd = 1 g cm−3, with the equivalent monodisperse distribution (dashed lines). Overall, the exponential distribution enhances the representation of smaller particle sizes relative to the monodisperse distribution. The particle mass distribution is spread over more orders of magnitude than the particle size distribution.

3 Condensation with exponential distribution

In this section, we derive the rate of condensation and evaporation of mass from the cloud bulk when assuming an exponential distribution and compare it to the monodisperse derivation. Following the differentiation with time for the method of moments (e.g. Gail & Sedlmayr 2013), the change in time of the moments for the condensation and evaporation process is given by $$\frac{d{M}^{(k)}}{dt}=k{\int }_0^{\mathrm{\infty }}{m}^{k-1}\frac{dm}{dt}f(m)dm,$$(10)

where dm/dt [g s⁻¹] is the rate of change in the cloud particle mass. In the monodisperse size distribution, the particle size distribution is assumed to be a delta function at the mean particle size expressed formally as $$f(m)=N_{\mathrm{c}}\delta (m-\lambda ).$$(11)

dm/dt is then independent of mass, removing it from the integral in Eq. (10). This leads to the monodisperse moment method typically used in the literature (e.g. Woitke & Helling 2003; Helling et al. 2008; Gail & Sedlmayr 2013; Ohno & Okuzumi 2017; Lee 2023), whereby the rate of change in a moment, k, is $$\frac{d{M}^{(k)}}{dt}=k\frac{dm}{dt}{\int }_0^{\mathrm{\infty }}{m}^{k-1}f(m)dm=k\frac{dm}{dt}{M}^{(k-1)}.$$(12)

For example, the monodisperse size distribution source term for condensational growth or evaporation for the first moment (k = 1) is $$\frac{d{\rho }_{\mathrm{c}}}{dt}={\frac{dm}{dt}|}_{r_{\mathrm{c}}}{N}_{\mathrm{c}},$$(13)

where dm/dt is given by Eq. (21) or (30), dependent on the Knudsen number regime and evaluated at the mean particle size, r_c.

We define the Knudsen number, Kn, as $$\mathrm{K}\mathrm{n}=\frac{{\lambda }_{\mathrm{a}}}{r},$$(14)

where r [cm] is the particle radius, and λ_a [cm] is the atmospheric mean path length, which can be given by $${\lambda }_{\mathrm{a}}=\frac{2{\eta }_{\mathrm{a}}}{{\rho }_{\mathrm{a}}}\sqrt{\frac{\pi {\overline{\mu }}_{\mathrm{a}}}{8RT}}.$$(15)

η_a [g cm⁻¹ s⁻¹] is the atmospheric dynamical viscosity, which can be calculated following Lee & Ohno (2025) using the Rosner (2012) fitting function and Davidson (1993) mixing rule, ρ_a [g cm⁻³] the atmospheric mass density, and ${\overline{\mu }}_{\mathrm{a}}\left[\mathrm{g}{\mathrm{m}\mathrm{o}\mathrm{l}}^{-1}\right]$ the mean molecular weight of the atmosphere.

The monodisperse Knudsen number is given by $$\overline{\mathrm{K}\mathrm{n}}=\frac{{\lambda }_{\mathrm{a}}}{r_{\mathrm{c}}},$$(16)

and we can define a population number averaged Knudsen number, $\overline{{\mathrm{K}\mathrm{n}}_N}$, as $$\overline{{\mathrm{K}\mathrm{n}}_N}=\frac{{\int }_0^{\mathrm{\infty }}\mathrm{K}\mathrm{n}(m)f(m)dm}{{\int }_0^{\mathrm{\infty }}f(m)dm}=\frac{{\lambda }_{\mathrm{a}}s{\int }_0^{\mathrm{\infty }}{m}^{-1/3}f(m)dm}{N_{\mathrm{c}}},$$(17)

where s = (3/4πρ_d)^−1/3. Using the definition of the moments and the Eq. (4) moment generator, this results in the expression $$\overline{{\mathrm{K}\mathrm{n}}_N}=\frac{{\lambda }_{\mathrm{a}}s{M}^{(-1/3)}}{N_{\mathrm{c}}}=\frac{{\lambda }_{\mathrm{a}}}{r_{\mathrm{c}}}\mathrm{\Gamma }\left(\frac{2}{3}\right)=\overline{\mathrm{K}\mathrm{n}}\mathrm{\Gamma }\left(\frac{2}{3}\right),$$(18)

where M^(−1/3) is the k = −1/3 non-integer moment power. Similarly, the mass-weighted Knudsen number, ${\overline{\mathrm{K}\mathrm{n}}}_m$, is given by $$\overline{{\mathrm{K}\mathrm{n}}_m}=\frac{{\int }_0^{\mathrm{\infty }}\mathrm{K}\mathrm{n}(m)mf(m)dm}{{\int }_0^{\mathrm{\infty }}mf(m)dm}=\frac{{\lambda }_{\mathrm{a}}s{\int }_0^{\mathrm{\infty }}{m}^{2/3}f(m)dm}{{\rho }_{\mathrm{c}}},$$(19)

with $$\overline{{\mathrm{K}\mathrm{n}}_m}=\frac{{\lambda }_{\mathrm{a}}s{M}^{(2/3)}}{{\rho }_{\mathrm{c}}}=\frac{{\lambda }_{\mathrm{a}}}{r_{\mathrm{c}}}\mathrm{\Gamma }\left(\frac{5}{3}\right)=\overline{\mathrm{K}\mathrm{n}}\mathrm{\Gamma }\left(\frac{5}{3}\right),$$(20)

where M^(2/3) is the k = 2/3 non-integer moment power.

3.1 Kn ≪ 1 regime

For the exponential distribution (and others), unlike for the monodisperse distribution, the mass dependence in the growth rate cannot be ignored in the integral of Eq. (10). The change in mass for condensation and evaporation in the continuum regime (Kn ≪ 1) is given by (e.g. Woitke & Helling 2003) $${\left(\frac{dm}{dt}\right)}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}=4\pi rD{m}_{\mathrm{v}}{n}_{\mathrm{v}}(1-\frac{1}{S}),$$(21)

where r [cm] is the radius of the particle, m_v [g] the mass of one unit of vapour, n_v [cm⁻³] the number density of the vapour, and S the supersaturation ratio of the vapour. D [cm² s⁻¹] is the diffusive rate of the vapour onto the particle surface and can be given by (e.g. Jacobson 2005) $$D=\frac{5}{16N_{\mathrm{A}}{d}_i^2{\rho }_{\mathrm{a}}}\sqrt{\frac{RT{\overline{\mu }}_{\mathrm{a}}}{2\pi }\left(\frac{{\mu }_{\mathrm{v}}+{\overline{\mu }}_{\mathrm{a}}}{{\mu }_{\mathrm{v}}}\right)},$$(22)

where N_A is Avogadro’s number, d_i [cm] the collision diameter, and μ_v [g mol⁻¹] the molecular weight of the condensable vapour.

The particle radius is converted to particle mass through the relation $$r={\left(\frac{3m}{4\pi {\rho }_{\mathrm{d}}}\right)}^{1/3},$$(23)

which reveals the mass dependence, ∝m^1/3, of condensational growth and evaporation in the continuum regime. Defining the mass independent kinetic pre-factor as C₀, $$C_0=4\pi {\left(\frac{3}{4\pi {\rho }_{\mathrm{d}}}\right)}^{1/3}D{m}_{\mathrm{v}}{n}_{\mathrm{v}}(1-\frac{1}{S}),$$(24)

Eq. (10) becomes $$\frac{d{M}^{(k)}}{dt}=k{C}_0{\int }_0^{\mathrm{\infty }}{m}^{1/3}{m}^{k-1}f(m)dm=k{C}_0{M}^{(k-2/3)}.$$(25)

3.1.1 Monodisperse solution

For the monodisperse assumption, Eq. (12) was used. This leads to the change in the first moment (k = 1) being $${\left(\frac{d{\rho }_{\mathrm{c}}}{dt}\right)}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}=4\pi r_{\mathrm{c}}D{m}_{\mathrm{v}}{n}_{\mathrm{v}}(1-\frac{1}{S})N_{\mathrm{c}}.$$(26)

3.1.2 Exponential distribution solution

Using the exponential distribution moment generator (Eq. (4)), Eq. (10) is now $$\frac{d{M}^{(k)}}{dt}=k{C}_0{\lambda }^{k-2/3}{N}_{\mathrm{c}}\mathrm{\Gamma }(k+1/3).$$(27)

We can then use Eq. (4) for the k − 1 moment to factor out N_c and get the expression in terms of the previous moment. $$\frac{d{M}^{(k)}}{dt}=k{C}_0{\lambda }^{1/3}{M}^{(k-1)}\frac{\mathrm{\Gamma }(k+1/3)}{\mathrm{\Gamma }(k)}.$$(28)

The full expression for the first moment (k = 1) is then $${\left(\frac{d{\rho }_{\mathrm{c}}}{dt}\right)}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}=4\pi r_{\mathrm{c}}D{m}_{\mathrm{v}}{n}_{\mathrm{v}}(1-\frac{1}{S})\mathrm{\Gamma }\left(\frac{4}{3}\right)N_{\mathrm{c}},$$(29)

which reveals that the rate of change in the first moment in the Kn ≪ 1 regime for the exponential distribution is a factor of Γ(4/3) different to the monodisperse rate (Eq. (27)).

3.2 Kn ≫ 1 regime

In the large Knudsen number regime, Kn ≫ 1, the free molecular or kinetic regime, the mass dependence is ∝m^2/3 as the condensation rate expression is given by (e.g. Woitke & Helling 2003), $${\left(\frac{dm}{dt}\right)}^{\text{free\hspace{0.17em}}}=4\pi r^2{v}_{\text{th\hspace{0.17em}}}{m}_{\mathrm{v}}{n}_{\mathrm{v}}\alpha (1-\frac{1}{S}),$$(30)

where $v_{\text{th\hspace{0.17em}}}=\sqrt{k_{\mathrm{b}}T/2\pi m_v}\left[\mathrm{c}\mathrm{m}{\mathrm{s}}^{-1}\right]$ is the thermal velocity of the condensable vapour, with m_v [g] the mass of the condensable vapour. α = [0,1] is the ‘sticking efficiency’ of the vapour to the particle surface. Again, defining the constant kinetic pre-factor, C₀, as $$C_0=4\pi {\left(\frac{3}{4\pi {\rho }_{\mathrm{d}}}\right)}^{2/3}{v}_{\mathrm{t}\mathrm{h}}{m}_{\mathrm{v}}{n}_{\mathrm{v}}\alpha (1-\frac{1}{S}),$$(31)

Eq. (10) becomes $$\frac{d{M}^{(k)}}{dt}=k{C}_0{\int }_0^{\mathrm{\infty }}{m}^{2/3}{m}^{k-1}f(m)dm=k{C}_0{M}^{(k-1/3)}.$$(32)

3.2.1 Monodisperse solution

As in the Kn ≪ 1 regime, (12) was used for the monodisperse distribution. This leads to the change in the first moment (k = 1) being $${\left(\frac{d{\rho }_{\mathrm{c}}}{dt}\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}_{\mathrm{v}}{n}_{\mathrm{v}}(1-\frac{1}{S})N_{\mathrm{c}}.$$(33)

3.2.2 Exponential distribution solution

Similar to the Kn ≪ 1 regime solution, after using the moment generator equation (Eq. (4)) the moment expression is given as $$\frac{d{M}^{(k)}}{dt}=k{C}_0{\lambda }^{k-1/3}{N}_{\mathrm{c}}\mathrm{\Gamma }(k+2/3),$$(34)

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

This gives for the full first moment (k = 1) expression $${\left(\frac{d{\rho }_{\mathrm{c}}}{dt}\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}_{\mathrm{v}}{n}_{\mathrm{v}}(1-\frac{1}{S})\mathrm{\Gamma }\left(\frac{5}{3}\right)N_{\mathrm{c}},$$(36)

which reveals that the first moment in the Kn ≫ 1 regime is a factor of Γ(5/3) different when assuming an exponential size distribution compared to the monodisperse distribution.

Overall, the difference between the monodisperse rate of change and exponential assumed distribution for the first moment is the factors of gamma, Γ(4/3) and Γ(5/3), which are both <1. The rate of condensation in the exponential distribution is therefore a factor of ~0.9 slower than the monodisperse distribution. This makes intuitive sense, as the monodisperse distribution will weight the number of particles at the mean value, as opposed to the exponential distribution, which uses weights in such a way as to get a better representation of smaller particles in its formulation. Combined with the exponential distributions’ strong drop-off in particle size after the mean value, the condensation rate for the exponential distribution is expected to be lower than the monodisperse distribution.

Our derived analytical solution makes it trivial to include an exponential particle size distribution dependence in current exoplanet atmosphere cloud condensation codes that use the moment method. Analogous derivations in particle radius or volume space instead of mass are a possibility for theories and models that apply the moments in radius or volume space (e.g. Helling et al. 2008; Gail & Sedlmayr 2013).

3.3 Condensation regime interpolation function

Equations (29) and (36) are both valid in their respective Knudsen number regimes of Kn ≪ 1 and Kn ≫ 1, respectively. However, in the intermediate regime, Kn ~ 1, a direct expression for the mass condensation and evaporation rate equation is difficult to derive.

To address this, following Woitke & Helling (2003), we applied an interpolation function that transitions smoothly between the diffusive and free molecular regimes. When the expressions Eqs. (29) and (36) are equal, this results in (Woitke & Helling 2003) $$\alpha v_{\mathrm{t}\mathrm{h}}r=D.$$(37)

Using the Knudsen number definition, the critical Knudsen number, Kn_cr, for the monodisperse distribution (r = r_c) is given by $${\mathrm{K}\mathrm{n}}_{\mathrm{c}\mathrm{r}}=\frac{{\lambda }_{\mathrm{a}}\alpha v_{\mathrm{t}\mathrm{h}}}{D},$$(38)

and for the population-averaged exponential distribution is $$\overline{{\mathrm{K}\mathrm{n}}_{\mathrm{c}\mathrm{r}}}=\frac{{\lambda }_{\mathrm{a}}\alpha v_{\mathrm{t}\mathrm{h}}\mathrm{\Gamma }(5/3)}{D\mathrm{\Gamma }(4/3)}.$$(39)

The critical Knudsen number can then be used to interpolate between each regime near the critical value (e.g. Woitke & Helling 2003). We propose a symmetric smooth tanh function of the form $$f\left({\mathrm{K}\mathrm{n}}^{\mathrm{\prime }}\right)=\frac{1}{2}[1-\tanh (a{\log }_{10}{\mathrm{K}\mathrm{n}}^{\mathrm{\prime }})],$$(40)

where ${\mathrm{K}\mathrm{n}}^{\mathrm{\prime }}=\overline{\mathrm{K}\mathrm{n}}/{\mathrm{K}\mathrm{n}}_{\mathrm{c}\mathrm{r}}$ for the monodisperse formulation, or $\overline{{\mathrm{K}\mathrm{n}}^{\mathrm{\prime }}}=\overline{{\mathrm{K}\mathrm{n}}_m}/\overline{{\mathrm{K}\mathrm{n}}_{\mathrm{c}\mathrm{r}}}$ for the population-averaged value, and has a range between 0 and 1, with $\overline{{\mathrm{K}\mathrm{n}}_m}$ the mass-weighted population-averaged Knudsen number (Eq. (20)). For the monodisperse distribution, the net change in the first moment is then a linear combination using this interpolation function, $${\left(\frac{d{\rho }_{\mathrm{c}}}{dt}\right)}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}}={f\left({\mathrm{K}\mathrm{n}}^{\mathrm{\prime }}\right)\frac{dm}{dt}|}_{r_{\mathrm{c}}}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}{N}_{\mathrm{c}}+{[1-f\left({\mathrm{K}\mathrm{n}}^{\mathrm{\prime }}\right)]\frac{dm}{dt}|}_{r_{\mathrm{c}}}^{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}{N}_{\mathrm{c}}.$$(41)

We set a scaling value of a = 2, which smoothly transitions between the asymptotic expressions at around two orders of magnitude in the Knudsen number away from the critical number. At one order of magnitude away from the critical number, the interpolation function gives a value of f(Kn′) ≈ 0.1 or 0.9. For consideration of the exponential distribution, the respective gamma factors for each regime can be easily re-included in Eq. (41), $$\begin{array}{rl}{\left(\frac{d{\rho }_{\mathrm{c}}}{dt}\right)}_{\text{cond\hspace{0.17em}}}=& {f\left(\overline{{\mathrm{K}\mathrm{n}}^{\mathrm{\prime }}}\right)\frac{dm}{dt}|}_{r_{\mathrm{c}}}^{\text{diff\hspace{0.17em}}}\mathrm{\Gamma }\left(\frac{4}{3}\right)N_{\mathrm{c}}\\ & +{[1-f\left(\overline{{\mathrm{K}\mathrm{n}}^{\mathrm{\prime }}}\right)]\frac{dm}{dt}|}_{r_{\mathrm{c}}}^{\text{free\hspace{0.17em}}}\mathrm{\Gamma }\left(\frac{5}{3}\right)N_{\mathrm{c}}.\end{array}$$(42)

4 Collisional growth with exponential distribution

Another important process in shaping the particle size distribution is the collision and subsequent growth of particles, either through Brownian motion of particles (coagulation) or relative settling velocities of different sized particles (coalescence). In Drake (1972), the mass moment generator expression of the Smoluchowski equation for collisional growth was derived as $$\begin{array}{rl}\frac{d{M}^{(k)}}{dt}=& \frac{1}{2}{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}K(m,m^{\mathrm{\prime }})f(m)f\left(m^{\mathrm{\prime }}\right)\\ & \times [{(m+m^{\mathrm{\prime }})}^k-m^k-m^{\mathrm{\prime }k}]dmd{m}^{\mathrm{\prime }},\end{array}$$(43)

where m and m′ represent different masses in the particle mass distribution and K(m, m′) [cm³ s⁻¹] the kernel function. We can immediately find the expression for the zeroth moment as $$\frac{d{M}^{(0)}}{dt}=-\frac{1}{2}{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}K(m,m^{\mathrm{\prime }})f(m)f\left(m^{\mathrm{\prime }}\right)dmd{m}^{\mathrm{\prime }},$$(44)

which in general serves as the basis for deriving the collisional growth rate expressions using moment methods.

In Lee & Ohno (2025), following arguments in Rossow (1978), it was assumed that the size distribution is monodisperse, allowing the kernel function to be moved outside the integral. In this study, we applied the method found in Morán & Kholghy (2023, Appendix C) to find size distribution averaged collision kernels, which we repeat below for completeness. Following Eq. (44), in Morán & Kholghy (2023) the number density population-averaged collision kernel is given as $$\begin{array}{rl}\overline{K(m,m^{\mathrm{\prime }})}& =\frac{{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}K(m,m^{\mathrm{\prime }})f(m)f\left(m^{\mathrm{\prime }}\right)dmd{m}^{\mathrm{\prime }}}{{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}f(m)f\left(m^{\mathrm{\prime }}\right)dmd{m}^{\mathrm{\prime }}}\\ & =\frac{1}{N_{\mathrm{c}}^2}{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}K(m,m^{\mathrm{\prime }})f(m)f\left(m^{\mathrm{\prime }}\right)dmd{m}^{\mathrm{\prime }}.\end{array}$$(45)

In the following sections, we use the Morán & Kholghy (2023) technique to derive population average collision kernels for both the Kn ≪ 1 and Kn ≫ 1 regimes. We then find collisional growth rate expressions for the monodisperse and exponential distribution cases.

4.1 Kn ≪ 1 regime

For collisional growth driven by Brownian motion (coagulation) in the Kn ≪ 1 regime, the collisional kernel is given by (Chandrasekhar 1943) $$K(r,r^{\mathrm{\prime }})=4\pi [D(r)+D\left(r^{\mathrm{\prime }}\right)](r+r^{\mathrm{\prime }}),$$(46)

where r and r′ represent different particle radii in the distribution and D(r) is the particle diffusion factor given by (Chandrasekhar 1943) $$D(r)=\frac{k_{\mathrm{b}}T\beta }{6\pi {\eta }_{\mathrm{a}}r}.$$(47)

In this, β, the Cunningham slip factor, generally depends on the radius (and therefore mass) through scaling with the Knudsen number in the functional form $$\beta =1+\mathrm{Kn}[A+B{e}^{-C/\mathrm{Kn}}].$$(48)

In this study, we used the parameters experimentally measured between Kn ≈ 0.5–83 from Kim et al. (2005) $$\beta =1+\mathrm{K}\mathrm{n}[1.165+0.483e^{-0.997/\mathrm{K}\mathrm{n}}],$$(49)

which is approximated to within ≈10% in the linear form $$\beta =1+A\mathrm{K}\mathrm{n},$$(50)

where A = 1.639, calculated from an optimised least-squares fitting to the Kim et al. (2005) coefficients.

Converting Eq. (46) to mass units gives the mass dependence of the kernel, $$K(m,m^{\mathrm{\prime }})=K_0(\beta m^{-1/3}+{\beta }^{\mathrm{\prime }}{m}^{\mathrm{\prime }-1/3})(m^{1/3}+m^{\mathrm{\prime }1/3}),$$(51)

where K₀ is the constant kinetic pre-factor given by $$K_0=\frac{2k_{\mathrm{b}}T}{3{\eta }_{\mathrm{a}}}.$$(52)

4.1.1 Monodisperse distribution solution

If we make the monodisperse assumption for the particle mass distribution, where f(m) ≈ N_cδ(m − λ) and f(m′) ≈ N_cδ(m′ − λ), the kernel function becomes $$K(m,m^{\mathrm{\prime }})=4K_0\overline{\beta },$$(53)

where $\overline{\beta }$ is Eq. (49) evaluated at the mean particle size, leading to a zeroth moment expression given by $$\frac{d{M}^{(0)}}{dt}=-2K_0\overline{\beta }{N}_{\mathrm{c}}^2.$$(54)

This results in the same equation from Lee & Ohno (2025) for the number density collisional growth rate, $${\left(\frac{d{N}_{\mathrm{c}}}{dt}\right)}_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{g}}\approx -\frac{4k_{\mathrm{b}}T\overline{\beta }}{3{\eta }_{\mathrm{a}}}{N}_{\mathrm{c}}^2,$$(55)

recovering the monodisperse limit.

4.1.2 Exponential distribution solution

We now attempt a zeroth moment solution, assuming an exponential size distribution. Following Morán & Kholghy (2023), the population-averaged kernel in this regime is given by $$\begin{array}{rl}\overline{K(m,m^{\mathrm{\prime }})}=& \frac{1}{N_{\mathrm{c}}^2}{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}K(m,m^{\mathrm{\prime }})f(m)f\left(m^{\mathrm{\prime }}\right)dmd{m}^{\mathrm{\prime }}\\ =& \frac{K_0}{N_{\mathrm{c}}^2}{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}(\beta m^{-1/3}+{\beta }^{\mathrm{\prime }}{m}^{\mathrm{\prime }-1/3})(m^{1/3}+m^{\mathrm{\prime }1/3})\\ & \cdot f(m)f\left(m^{\mathrm{\prime }}\right)dmd{m}^{\mathrm{\prime }}.\end{array}$$(56)

Assuming that β and β′ are well represented at the mean particle size, this can be converted into a function of non-integer moment powers as (Morán & Kholghy 2023) $$\overline{K(m,m^{\mathrm{\prime }})}=\frac{2K_0\overline{\beta }}{N_{\mathrm{c}}^2}(M^{(0)}{M}^{(0)}+M^{(1/3)}{M}^{(-1/3)}),$$(57)

where M^(1/3) and M^(−1/3) are the k = 1/3 and k = −1/3 non-integer power moment, respectively. Using the moment generator for the exponential mass distribution (Eq. (4)), we can immediately find the population-averaged kernel for the exponential distribution in the Kn ≪ 1 regime, $$\overline{K(m,m^{\mathrm{\prime }})}=2K_0\overline{\beta }[1+\mathrm{\Gamma }\left(\frac{4}{3}\right)\mathrm{\Gamma }\left(\frac{2}{3}\right)].$$(58)

Putting this averaged kernel into Eq. (44) gives the change in the number density as $${\left(\frac{d{N}_{\mathrm{c}}}{dt}\right)}_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{g}}\approx -\frac{2k_{\mathrm{b}}T\overline{\beta }}{3{\eta }_{\mathrm{a}}}[1+\mathrm{\Gamma }\left(\frac{4}{3}\right)\mathrm{\Gamma }\left(\frac{2}{3}\right)]N_{\mathrm{c}}^2.$$(59)

Evaluating the bracket containing the gamma functions gives ≈2.21, meaning that the rate of change in number density from collisional growth assuming an exponential distribution is therefore around ≈1.1 times larger than the monodisperse assumption in the Kn ≪ 1 regime.

4.1.3 Accounting for the Cunningham slip factor

Assuming the Cunningham slip factor in its linear form (Eq. (50)) enables a moment analytical solution including the effects of slip correction. Integrating including the β and β′ factors gives the population-averaged kernel as $$\begin{array}{rl}\overline{K(m,m^{\mathrm{\prime }})}=& \frac{2K_0}{N_{\mathrm{c}}^2}[M^{(0)}{M}^{(0)}+M^{(1/3)}{M}^{(-1/3)}\\ & +{\lambda }_{\mathrm{a}}A{\left(\frac{3}{4\pi {\rho }_{\mathrm{d}}}\right)}^{-1/3}(M^{(0)}{M}^{(-1/3)}+M^{(1/3)}{M}^{(-2/3)})],\end{array}$$(60)

which is Eq. (58) but with an additional moment-dependent correction term. The change in number density with time is then $$\begin{array}{rl}{\left(\frac{d{N}_{\mathrm{c}}}{dt}\right)}_{\text{coag\hspace{0.17em}}}\approx & -\frac{2k_{\mathrm{b}}T}{3{\eta }_{\mathrm{a}}}\{1+\mathrm{\Gamma }\left(\frac{4}{3}\right)\mathrm{\Gamma }\left(\frac{2}{3}\right)\\ & +\frac{{\lambda }_{\mathrm{a}}}{r_{\mathrm{c}}}A[\mathrm{\Gamma }\left(\frac{2}{3}\right)+\mathrm{\Gamma }\left(\frac{4}{3}\right)\mathrm{\Gamma }\left(\frac{1}{3}\right)]\}{N}_{\mathrm{c}}^2.\end{array}$$(61)

It is more difficult to ascertain the exact difference between the monodisperse rate and this formulation. However, assuming Kn ≈ 1, the rate is increased by a factor of ≈1.78 from the monodisperse formulation.

4.2 Kn ≫ 1 regime

For Brownian coagulation in the Kn ≫ 1 regime, the collisional kernel is given by (Jacobson 2005) $$K(r,r^{\mathrm{\prime }})=\pi {(r+r^{\mathrm{\prime }})}^2\sqrt{(v(m{)}^2+v{\left(m^{\mathrm{\prime }}\right)}^2)},$$(62)

where v is the thermal velocity of the particle given by $$v=\sqrt{\frac{8k_{\mathrm{b}}T}{\pi m}}.$$(63)

Converting the kernel to mass units gives $$K(m,m^{\mathrm{\prime }})=K_0{(m^{1/3}+m^{\mathrm{\prime }1/3})}^2{(m^{-1}+m^{\mathrm{\prime }-1})}^{1/2},$$(64)

where K₀ is the kinetic pre-factor, $$K_0={\left(\frac{3}{4\pi {\rho }_{\mathrm{d}}}\right)}^{2/3}\sqrt{8\pi k_{\mathrm{b}}T}.$$(65)

To make the integration analytically tractable later, we approximated the square root term as separable $${(m^{-1}+m^{\mathrm{\prime }-1})}^{1/2}\approx H(m^{-1/2}+m^{\mathrm{\prime }-1/2}),$$(66)

where H is a fitting factor. H should be in the range of [$1/\sqrt{2}$, 1], where $1/\sqrt{2}$ corresponds to the monodisperse case and 1 a highly broadened size distribution. We used a Brent–Dekker root-finding method to find the value of H that minimises the relative error across a typical Kn ≫ 1 mass range of 10⁻²⁰−10⁻¹⁴ g (r ~ 1 nm–0.1 μm), yielding H ≈ 0.85, which is approximately the mid value between $1/\sqrt{2}$ and 1.

4.2.1 Monodisperse solution

For the monodisperse solution, the mass distribution takes the form f(m) ≈ N_cδ(m − λ) and f(m′) ≈ N_cδ(m′ − λ). The kernel now becomes $$K(m,m^{\mathrm{\prime }})=4\sqrt{\frac{2}{\lambda }}{K}_0{\lambda }^{2/3},$$(67)

leading to a zeroth moment expression as $$\frac{d{M}^{(0)}}{dt}=-2\sqrt{\frac{2}{\lambda }}{K}_0{\lambda }^{2/3}.$$(68)

The full equation for the coagulation rate is then $${\left(\frac{d{N}_{\mathrm{c}}}{dt}\right)}_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{g}}\approx -2\sqrt{\frac{2}{m_{\mathrm{c}}}}{\left(\frac{3}{4\pi {\rho }_{\mathrm{d}}}\right)}^{2/3}\sqrt{8\pi k_{\mathrm{b}}T}{m}_{\mathrm{c}}^{2/3}{N}_{\mathrm{c}}^2,$$(69)

and simplifies to $${\left(\frac{d{N}_{\mathrm{c}}}{dt}\right)}_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{g}}\approx -8\sqrt{\frac{\pi k_{\mathrm{b}}T}{m_{\mathrm{c}}}}{r}_{\mathrm{c}}^2{N}_{\mathrm{c}}^2,$$(70)

which is the same expression found in previous monodisperse studies (e.g. Lee & Ohno 2025).

4.2.2 Exponential distribution solution

Finding a fully consistent analytical solution assuming an exponential distribution is difficult in this regime, but we can follow the same method in Morán & Kholghy (2023) to derive a population-averaged kernel function. Placing Eq. (64) into Eq. (45) gives $$\begin{array}{rl}\overline{K(m,m^{\mathrm{\prime }})}=& \frac{1}{N_{\mathrm{c}}^2}{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}K(m,m^{\mathrm{\prime }})f(m)f\left(m^{\mathrm{\prime }}\right)dmd{m}^{\mathrm{\prime }}\\ =& \frac{K_{0}H}{N_{\mathrm{c}}^2}{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}{(m^{1/3}+m^{\mathrm{\prime }1/3})}^2(m^{-1/2}+m^{\mathrm{\prime }-1/2})\\ & \cdot f(m)f\left(m^{\mathrm{\prime }}\right)dmd{m}^{\mathrm{\prime }}.\end{array}$$(71)

This results in a set of six integrals, the solutions of which are non-integer powers of the moments. After integrating and some algebra, the end result is $$\begin{array}{rl}\overline{K(m,m^{\mathrm{\prime }})}& =\frac{2K_{0}H}{N_{\mathrm{c}}^2}\\ & \cdot (M^{(2/3)}{M}^{(-1/2)}+2M^{(1/3)}{M}^{(-1/6)}+M^{(1/6)}{M}^{(0)}),\end{array}$$(72)

where M^(2/3), M^(−1/2), M^(1/3), M^(−1/6), and M^(1/6) are the k = 2/3, k = −1/2, k = 1/3, k = −1/6, and k = 1/6 non-integer powers of the moments, respectively. Using the moment generator for the exponential mass distribution (Eq. (4)) results in $$\overline{K(m,m^{\mathrm{\prime }})}=2K_{0}H{\lambda }^{1/6}[\mathrm{\Gamma }\left(\frac{5}{3}\right)\mathrm{\Gamma }\left(\frac{1}{2}\right)+2\mathrm{\Gamma }\left(\frac{4}{3}\right)\mathrm{\Gamma }\left(\frac{5}{6}\right)+\mathrm{\Gamma }\left(\frac{7}{6}\right)].$$(73)

Using the relation λ^1/6 = λ^2/3 λ^−1/2, the final expression for the change in number density can be given in a similar form to the monodisperse relation: $$\begin{array}{rl}{\left(\frac{d{N}_{\mathrm{c}}}{dt}\right)}_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{g}}\approx & -0.85\sqrt{\frac{8\pi k_{\mathrm{b}}T}{m_{\mathrm{c}}}}{r}_{\mathrm{c}}^2\\ & \cdot [\mathrm{\Gamma }\left(\frac{5}{3}\right)\mathrm{\Gamma }\left(\frac{1}{2}\right)+2\mathrm{\Gamma }\left(\frac{4}{3}\right)\mathrm{\Gamma }\left(\frac{5}{6}\right)+\mathrm{\Gamma }\left(\frac{7}{6}\right)]N_{\mathrm{c}}^2.\end{array}$$(74)

Overall, the exponential distribution collision rate is a factor of ≈1.37 faster than the monodisperse formulation in the Kn ≫ 1 regime.

4.3 Collisional regime interpolation function

Equations (59) and (74) are valid in the Kn ≪ 1 and Kn ≫ 1 limits, respectively; however, in many cases particles occupy an intermediate regime Kn ~ 1 where it is difficult to derive an analytical Brownian coagulation kernel. To solve this problem, many studies apply an interpolation function between the limiting regimes, with the most common one following the Fuchs (1964) prescription (e.g. Lavvas et al. 2011; Gao et al. 2018).

In this study, we followed Morán (2022), who propose a simple and accurate interpolation function between the continuum or diffusive (Kn ≪ 1) and free molecular or kinetic (Kn ≫ 1) regimes. This is designed as a modification to the diffusive regime kernel as $$K(r,r^{\mathrm{\prime }})=4\pi [D(r)+D\left(r^{\mathrm{\prime }}\right)](r+r^{\mathrm{\prime }})g\left({\mathrm{K}\mathrm{n}}_{\mathrm{D}}\right),$$(75)

where Kn_D is the diffusive particle Knudsen number given by $${\mathrm{K}\mathrm{n}}_{\mathrm{D}}=\frac{8\sqrt{2}[D(r)+D\left(r^{\mathrm{\prime }}\right)]}{\pi \sqrt{v(r{)}^2+v{\left(r^{\mathrm{\prime }}\right)}^2}(r+r^{\mathrm{\prime }})},$$(76)

for a particle size pair. The interpolation function, g(Kn_D), derived in Morán (2022) is stated as $$g\left({\mathrm{K}\mathrm{n}}_{\mathrm{D}}\right)={(1+\frac{{\pi }^2}{8}{\mathrm{K}\mathrm{n}}_{\mathrm{D}}{}^2)}^{-1/2},$$(77)

which varies between a value of 1 for the diffusive/continuum regime as Kn_D → 0 and recovers the kinetic or free molecular kernel as Kn_D → ∞. We can calculate a population-averaged diffusive particle Knudsen number, $\overline{{\mathrm{K}\mathrm{n}}_{\mathrm{D}}}$, using the population-averaged kernels, $$\overline{{\mathrm{K}\mathrm{n}}_{\mathrm{D}}}=\frac{2\sqrt{2}}{\pi }\frac{\overline{K_1(m,m^{\mathrm{\prime }})}}{\overline{K_{\mathrm{h}}(m,m^{\mathrm{\prime }})}},$$(78)

where $\overline{K_1(m,m^{\mathrm{\prime }})}$ is the population-averaged kernel for the Kn ≪ 1 regime from Eq. (60) and $\overline{K_{\mathrm{h}}(m,m^{\mathrm{\prime }})}$ is the population-averaged kernel for the Kn ≫ 1 regime from Eq. (73).

For a monodisperse size distribution, this interpolation function can be readily used to modify Eq. (55): $${\left(\frac{d{N}_{\mathrm{c}}}{dt}\right)}_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{g}}\approx -\frac{4k_{\mathrm{b}}T\overline{\beta }}{3{\eta }_{\mathrm{a}}}{N}_{\mathrm{c}}^{2}g\left({\mathrm{K}\mathrm{n}}_{\mathrm{D}}\right),$$(79)

where Kn_D is evaluated at the mean grain size. For the exponential distribution, the overall coagulation rate is given by $$\begin{array}{rl}{\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}}}\{1+\mathrm{\Gamma }\left(\frac{2}{3}\right)\mathrm{\Gamma }\left(\frac{4}{3}\right)\\ & +\frac{{\lambda }_{\mathrm{a}}}{r_{\mathrm{c}}}A[\mathrm{\Gamma }\left(\frac{2}{3}\right)+\mathrm{\Gamma }\left(\frac{4}{3}\right)\mathrm{\Gamma }\left(\frac{1}{3}\right)]\}{N}_{\mathrm{c}}^{2}g\left(\overline{{\mathrm{K}\mathrm{n}}_{\mathrm{D}}}\right),\end{array}$$(80)

which recovers the Kn ≫ 1 and Kn ≪ 1 collisional rate limits for the exponential distribution.

4.4 Gravitational coalescence

For gravitational coalescence, the collisional growth kernel is given by (e.g. Jacobson 2005) $$K(r,r^{\mathrm{\prime }})=\pi {(r+r^{\mathrm{\prime }})}^2|v_{\mathrm{f}}(r)-v_{\mathrm{f}}\left(r^{\mathrm{\prime }}\right)|E,$$(81)

where v_f [cm s⁻¹] is the settling velocity of the particle. Defining Δv_f = |v_f(r) − v_f(r′)|, we followed Ohno & Okuzumi (2017) and evaluated Δv_f with a parameter, ϵ, which estimates the relative velocity of the particles from the mean particle size settling velocity, v_f(r_c), giving Δv_f ≈ ϵv_f(r_c). This was taken as ϵ ≈ 0.5, following the results of Sato et al. (2016), who found this value to best reproduce the results of a collisional bin resolving model for grain growth in protoplanetary discs. E is a collisional efficiency factor, dependent on the particle Stokes number.

4.4.1 Monodisperse solution

For the monodisperse distribution, putting the kernel from Eq. (81) into Eq. (44) and evaluating at the mean particle size leads to the expression found in previous studies (e.g. Lee & Ohno 2025), $${\left(\frac{d{N}_{\mathrm{c}}}{dt}\right)}_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{l}}\approx -2\pi r_{\mathrm{c}}^2\mathrm{\Delta }{v}_{\mathrm{f}}E{N}_{\mathrm{c}}^2.$$(82)

In this case, E is given by the monodisperse properties with the Stokes number, Stk, given by $$\text{\hspace{0.17em}Stk\hspace{0.17em}}=\frac{\mathrm{\Delta }{v}_{\mathrm{f}}{v}_{\mathrm{f}}\left(r_{\mathrm{c}}\right)}{g{r}_{\mathrm{c}}}=\frac{ϵv_{\mathrm{f}}^2\left(r_{\mathrm{c}}\right)}{g{r}_{\mathrm{c}}}.$$(83)

E is then (Guillot et al. 2014) $$E=\{\begin{array}{ll}max[0,1-0.42{\mathrm{S}\mathrm{t}\mathrm{k}}^{-0.75}]& \overline{\mathrm{K}\mathrm{n}}<1\\ 1.& \overline{\mathrm{K}\mathrm{n}}\ge 1\end{array}.$$(84)

4.4.2 Exponential distribution

For an exponential distribution, we made the assumption that Δv_f and E are well defined at the mean particle size value, allowing them to be taken outside the integral, making the integral analytically more tractable. Defining the kinetic pre-factor as $$K_0={\left(\frac{3}{4\pi {\rho }_{\mathrm{d}}}\right)}^{2/3}\pi \mathrm{\Delta }{v}_{\mathrm{f}}\overline{E},$$(85)

the population-averaged kernel is expressed as $$\begin{array}{rl}\overline{K(m,m^{\mathrm{\prime }})}& =\frac{1}{N_{\mathrm{c}}^2}{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}K(m,m^{\mathrm{\prime }})f(m)f\left(m^{\mathrm{\prime }}\right)dmd{m}^{\mathrm{\prime }}\\ & =\frac{K_0}{N_{\mathrm{c}}^2}{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}{(m^{1/3}+m^{\mathrm{\prime }1/3})}^{2}f(m)f\left(m^{\mathrm{\prime }}\right)dmd{m}^{\mathrm{\prime }}.\end{array}$$(86)

Solving the integral results in $$\overline{K(m,m^{\mathrm{\prime }})}=\frac{2K_0}{N_{\mathrm{c}}^2}(M^{(2/3)}{M}^{(0)}+M^{(1/3)}{M}^{(1/3)}).$$(87)

In the same procedure as the coagulation kernels, assuming an exponential distribution results in the final expression being $${\left(\frac{d{N}_{\mathrm{c}}}{dt}\right)}_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{l}}\approx -\pi r_{\mathrm{c}}^2\mathrm{\Delta }{v}_{\mathrm{f}}\overline{E}[\mathrm{\Gamma }\left(\frac{5}{3}\right)+\mathrm{\Gamma }{\left(\frac{4}{3}\right)}^2]N_{\mathrm{c}}^2.$$(88)

For the polydisperse case, the population-averaged Stokes number, $\overline{\text{Stk}}$, was approximated using the number weighted mean radius, ⟨r⟩ [cm], rather than the mass weighted mean radius, r_c, $$\overline{\mathrm{S}\mathrm{t}\mathrm{k}}=\frac{\mathrm{\Delta }{v}_{\mathrm{f}}{v}_{\mathrm{f}}\left(r_{\mathrm{c}}\right)}{g⟨r⟩}=\frac{ϵv_{\mathrm{f}}^2\left(r_{\mathrm{c}}\right)}{g⟨r⟩},$$(89)

where the relation between the number and mass weighted radii is given by $$⟨r⟩=r_{\mathrm{c}}\mathrm{\Gamma }\left(\frac{4}{3}\right).$$(90)

$\overline{E}$ was then delineated using the number density averaged Knudsen number, $$\overline{E}=\{\begin{array}{ll}max[0,1-0.42{\overline{\mathrm{S}\mathrm{t}\mathrm{k}}}^{-0.75}]& \overline{{\mathrm{K}\mathrm{n}}_N}<1\\ 1.& \overline{{\mathrm{K}\mathrm{n}}_N}\ge 1\end{array}.$$(91)

Overall, assuming that the collisional efficiency factors are equal, the collisional growth rate from gravitational coalescence assuming an exponential distribution is ≈0.85 times slower than the monodisperse approximation.

4.5 Summary of collisional growth results

Similar to the condensation factor of 0.9, the >1.1 enhancement factor for Brownian motion in the Kn ≪ 1 regime makes intuitive sense for the exponential distribution as the distribution is more skewed towards smaller particles, with an exponential drop-off for particles larger than the mean particle size. The collision rate is therefore enhanced due to the consideration of the smaller particle population in the exponential distribution. For the Kn ≫ 1 regime, the increase by ≈1.37 also makes intuitive sense, as in this regime collisions are guided by the relative cross-section of the particles, which is reduced for smaller particles but more strongly counteracted by the increased thermal velocity of small particles. Similarly, for gravitational coalescence, the collisional rate is dependent on the relative cross-section of the particles, resulting in an overall reduction in the collisional growth rate of around ≈0.85 for the exponential distribution compared to the monodisperse distribution.

5 Settling velocity interpolation function

In the Kn ≪ 1 regime, the Stokes law with an included Reynolds number-dependent fitting factor can be used from Ohno & Okuzumi (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},$$(92)

where g [cm s s⁻²] is the gravitational acceleration, and we have re-included the additional buoyancy term from the original Stokes law derivation. However, in the Kn ≫ 1 regime, which is commonly a regime in which small particles in exoplanet atmospheres are found, the Epstein regime drag law should be used. Here, we can follow Woitke & Helling (2003), who used the Schaaf (1963) drag coefficients in the limiting case for settling velocities well below the atmospheric thermal velocity, v_f ≪ c_T, to derive $$v_{\mathrm{f}}(r)=\frac{\sqrt{\pi }g{\rho }_{\mathrm{d}}r}{2{\rho }_{\mathrm{a}}{c}_{\mathrm{T}}},$$(93)

where $c_{\mathrm{T}}=\sqrt{2k_{\mathrm{b}}T/m_{\mathrm{a}}}\left[\mathrm{c}\mathrm{m}{\mathrm{s}}^{-1}\right]$ is the thermal velocity of the atmosphere.

We can again use the same tanh function from Eq. (40) to smoothly interpolate between the two expressions, $$v_{\mathrm{f}}=f\left({\mathrm{K}\mathrm{n}}^{\mathrm{\prime }}\right)v_{\mathrm{f},\text{\hspace{0.17em}Stokes\hspace{0.17em}}}+[1-f\left({\mathrm{K}\mathrm{n}}^{\mathrm{\prime }}\right)]v_{\mathrm{f},\text{\hspace{0.17em}Epstein\hspace{0.17em}}},$$(94)

where Kn′ = Kn/Kn_cr. Choosing a suitable critical value for the settling velocity transition region is a tricky endeavor; for example, Woitke & Helling (2003) find a critical value of Kn_cr = 1/3 equating limits of their drag force derivations. Through trial and error and visual inspection, we find that Kn_cr = 1 with a scaling factor of a = 2 gives a good balance and smooth transition for the Knudsen number regime limit velocities. Figure 2 shows an example of this interpolation scheme for various Knudsen numbers for typical atmospheric properties.

Fig. 2 Example of the settling velocity tanh smooth interpolation scheme. The blue line shows the Stokes law values (Eq. (92)), the orange line the Epstein law values (Eq. (93)), and the dashed black line the interpolation function (Eq. (94)). The background gas is assumed to be T = 500 K, p = 10−3 bar, $\log \mathrm{g}=3.25,\overline{\mu }=2.33\mathrm{g}{\mathrm{m}\mathrm{o}\mathrm{l}}^{-1}$, and 100% H2.

6 Y-dwarf KCI cloud simple 1D example

In this section, we apply our new equation set to a simple 1D KCl Y-dwarf example. We followed the same set-up as in Gao et al. (2018), using the Y-dwarf temperature-pressure (T-p) profiles, log g, and constant K_zz parameters used in their study. We used a simple, second-order explicit 1D vertical diffusion and advection time stepping scheme, similar to the one used in recent Exo-FMS GCM simulations of brown dwarfs (Lee et al. 2024; Lee & Ohno 2025). For the cloud microphysics, we used the ODE system presented in Lee & Ohno (2025), which includes condensation and evaporation, collisional growth, and homogeneous nucleation, but with modifications for the exponential distribution rates from this study when appropriate. We assumed that the settling velocity for each moment is well defined at the mass-weighted radius (Eq. (8)) derived from the moment solutions at each timestep. We also assumed that the particles are well coupled to the turbulent flow, and so are also vertically diffused in the atmosphere at the given K_zz eddy diffusion coefficient (e.g Woitke et al. 2020). We integrated the 1D column in time for 10⁹ seconds, after which all variables show relative differences of less than 10⁻⁶ between timesteps. We plot the relative difference between the exponential and monodisperse results as $$\delta N_{\mathrm{c}}=\frac{N_{\mathrm{c}}(\exp )-N_{\mathrm{c}}(\text{\hspace{0.17em}mono\hspace{0.17em}})}{N_{\mathrm{c}}(\text{\hspace{0.17em}mono\hspace{0.17em}})},$$(95)

with the equivalent expression for the mean particle size.

Figure 3 shows our results for the 1D column model. In the log g = 3.25 case, we find around an increase of ~15% in the total number density in the exponential distribution compared to the monodisperse distribution, with the mean particle sizes differing by ~1–2%. We suggest that, due to the mean particle size being quite large in this case, ~10 μm, this increase in number density arises from the decrease in the coalescence rate between the monodisperse and exponential distributions. This is reflected in the slightly smaller particle sizes seen for the exponential case.

In the log g = 4.25 case, larger differences are seen between the size distribution assumptions, with around a 25% reduction in cloud particle number density in the exponential case compared to the monodisperse case in the upper atmosphere, but a >200% difference closer to the cloud base. A relative difference of ~15% is seen in the mean particle radius in the upper atmosphere, but this increases to >40% at the cloud base. Due to the small mean particle size of ~0.1 μm, this suggests that the increase in coagulation rate in the Kn ≫ 1 regime is responsible for this difference in the upper atmosphere, while slower condensation occurring near the cloud base accounts for the increased number density and decreased particle radii.

In both size distribution cases, the mass mixing ratio and vapour mixing ratio profile remain highly similar. This suggests that the total mass that is condensed is not markedly altered between each assumption, but the balance in mass distributed between the number density and particle size is changed.

Fig. 3 KCl cloud structures using the log g = 3.25 (left) and log g = 4.25 (right) 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 and the solid lines the exponential size distribution results. The top panel shows the temperature-pressure profiles with the mass mixing ratio of the condensate, qc, vapour, qv, and saturation point, qs. The middle panel shows the mean particle size, rc, and total number density, Nc. The bottom panel shows the relative difference in rc and Nc between the monodisperse and exponential cloud structure.

Fig. 4 Reconstructed exponential particle size distributions (Eq. (1)) from the Y-dwarf KCl case moment results (Left: log g = 3.25, Right: log g = 4.25) in Section 6 when the exponential distribution is assumed. The colour bar shows the atmospheric pressure of the distribution. This shows the limitations of the exponential distribution as limited flexibility is present due to the assumed variance, resulting in a highly similar distribution shape across the cloud layer.

7 Discussion

Several assumptions were made to simplify the coagulation collisional growth integrals analytically. Specifically, in the Kn ≪ 1 regime we assumed a linear Cunningham slip factor, which has around a 10% error near the Kn ~ 1 transition regime. Possibilities to improve accuracy in this regime would be to fit a quadratic equation to the Cunningham slip β function, which would still enable an analytic solution, but increase the number of gamma correction terms in the Kn ≪ 1 coagulation rate expression. In the Kn ≫ 1 regime, beyond fitting for the H factor, no further assumptions were required to produce a solution. For gravitational coalescence, we assumed the collisional efficiency factor was well defined at the mean particle size. However, this assumption may not be entirely accurate, particularly in the Kn < 1 regime, where efficiency factors can vary widely, including reaching zero. Should the efficiency factor vary significantly from the average value across the size distribution, we can expect a significantly altered overall coalescence rate compared to the ones presented here. Accounting for these additional effects analytically or through a numerical or interpolation approach and relaxing the assumptions may be possible in the future, but is beyond the scope of the current study.

Our assumption of an exponential size distribution has many limitations, especially when compared to more flexible size distributions such as the gamma or log-normal distribution. The single parameter limitation of the exponential distribution inherits the assumption that the overall particle size distribution shape does not change in time, with the variance remaining constant at the square of the average particle mass and across the time evolution of the cloud structure. The gamma and log-normal distribution contain an additional parameter compared to the one-parameter exponential distribution, allowing them to capture a much less strict and flexible size distribution shape. A major advantage of the more complex distributions would be allowing the representation of both narrow and broad size distributions through directly calculating the variance using the moment results, rather than the assumed variance of the exponential distribution. To visualise the exponential distribution limitations, in Figure 4 we show the reconstructed exponential distributions for the results in Section 6 that used the exponential equation set from this study. From this, it is clear that the overall shape of the distribution does not change much across the cloud structure, a direct result of the assumed variance in the exponential distribution, meaning representation of a narrower or broader size distribution is limited in the exponential distribution. This can be compared to when the variance can be flexible; for example, Helling et al. (2008) reconstructed a potential exponential distribution that produced both narrow and broad distributions with height. However, it should be noted that Helling et al. (2008) reconstructed the size distribution from a monodisperse microphysics method.

However, with its useful analytical properties and single parameter enabling a two-moment representation, rather than three or more moments for the more complex distributions, the exponential distribution offers a simple way to attempt to take into account polydispersity in the size distribution. With the equations presented here, modellers can be more self-consistent in exploring the effects of polydispersity, rather than reconstructing distributions from moment solutions obtained from the monodisperse equation set (e.g. Helling et al. 2008; Stark et al. 2015; Lee 2023). Since the exponential distribution is a special parameterised form of the gamma distribution, we suspect that designing a three-moment framework with the gamma distribution will take a similar form to that derived here but with additional correction terms to the monodisperse equation set.

Our results in Section 6 used the same Y-dwarf T-p profiles as Gao et al. (2018), which also allows us to compare the moment method results to the 1D CARMA bin-resolving model used in Gao et al. (2018), who also modelled pure KCl cloud particle formation. Despite the exact detailed condensation, nucleation, and collisional growth equations being different in this study from the ones in Gao et al. (2018), this serves as a useful comparison of the two methodologies. In Gao et al. (2018) in the log g = 3.25 case, the results of the presented moment method and CARMA are remarkably similar. Both models show a similar condensed mass fraction structure. The mean particle size also shows a good agreement in this case, hovering around the 10 μm range, but the CARMA results tend to produce slightly smaller particles in the upper atmosphere by a few factors in this case. However, it is difficult to compare one to one, since Gao et al. (2018) present the effective particle radius (weighted by area), whereas our results show the mass-weighted radius (weighted by mass or volume). Despite this, our results are quite consistent as the effective particle radius is always smaller than the mass-weighted radius by some factor in a particle size distribution. In the log g = 4.25 case, our results are much more different from Gao et al. (2018); in particular, Gao et al. (2018) shows a near constant effective particle radius of around 3 μm, whereas we produce particles at around ~0.08 μm. We suggest that this may be a case in which the bin model produces a bimodal distribution, whereby the particle size distribution is split between a small particle regime and a large particle regime, and unable to be captured by the current moment methodology. This skews the effective particle radius to the large particle population in the bin model, even if the number density of the large particles is significantly lower. We suspect that the moment method is capturing the small particle regime in this case but is unable to model the large particle population due to the assumption of a continuous distribution inherent in the method.

8 Conclusions

In this short paper, we have derived condensation, evaporation and collisional growth expressions for a two-moment cloud microphysical scheme applicable to sub-stellar atmospheric clouds. We have gone beyond the monodisperse assumptions commonly used in previous studies and found simple analytical and tractable solutions when assuming an exponential particle size distribution.

Our analysis found that assuming an exponential distribution for particles leads to simple analytical formulae, which, in effect, modify the monodisperse formulations with constant factors for the time-dependent moment equations. These are ≈0.9 for the condensation or evaporation rate (first moment) and >1.1 for Brownian collisional growth (zeroth moment) in the Kn ≪ 1 regime and ≈1.37 for Brownian collisional growth in the Kn ≫ 1 regime and ≈0.85 for the gravitational coalescence collisional rate. However, we note that several simplifying assumptions were made regarding integrating the collision growth equations in order to make the problem more analytically tractable. In addition, we proposed simple interpolation functions for evaluating the condensation, settling velocity, and collisional growth rates in the Kn ~ 1 transition regime.

In a simple 1D column model example, we found significant differences of a factor of >200% in number density and >40% in particle size between the monodisperse and exponential size distribution assumptions, with the differences dependent on the T-p profile and gravity of the simulation. We suggest that, in our particular cases, the total condensed mass in the cloud structure is similar for each assumption, but the balance between number density and particle sizes is greatly affected by the choice of size distribution assumption.

In future papers, we shall extend our two-moment framework to three moments to enable the self-consistent evolution of more complex, three-parameter size distributions, but these, such as the log-normal and gamma distributions, are still somewhat analytically tractable. Due to the exponential distribution case being a special case of the gamma distribution, we find it likely that derivations using the gamma distribution will follow a similar structure to the equations presented in this study.

Acknowledgements

E.K.H. Lee is supported by the CSH through the Bernoulli Fellowship. We thank K. Ohno for initial discussions on the topic. We thank the referee for suggestions that led to improving the consistency of the two-moment framework.

References

All Figures

Fig. 1 Visualisation of the exponential size distribution compared against a monodisperse size distribution in size units (left) and mass units (right). This shows exponential distributions with (solid lines) λ = 10−3 (blue), 10−2 (orange), 0.1 (green), and 1 (red) μm, assuming Nc = 1 cm−3 and ρd = 1 g cm−3, with the equivalent monodisperse distribution (dashed lines). Overall, the exponential distribution enhances the representation of smaller particle sizes relative to the monodisperse distribution. The particle mass distribution is spread over more orders of magnitude than the particle size distribution.

In the text

Fig. 2 Example of the settling velocity tanh smooth interpolation scheme. The blue line shows the Stokes law values (Eq. (92)), the orange line the Epstein law values (Eq. (93)), and the dashed black line the interpolation function (Eq. (94)). The background gas is assumed to be T = 500 K, p = 10−3 bar, $\log \mathrm{g}=3.25,\overline{\mu }=2.33\mathrm{g}{\mathrm{m}\mathrm{o}\mathrm{l}}^{-1}$, and 100% H2.

In the text

Fig. 3 KCl cloud structures using the log g = 3.25 (left) and log g = 4.25 (right) 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 and the solid lines the exponential size distribution results. The top panel shows the temperature-pressure profiles with the mass mixing ratio of the condensate, qc, vapour, qv, and saturation point, qs. The middle panel shows the mean particle size, rc, and total number density, Nc. The bottom panel shows the relative difference in rc and Nc between the monodisperse and exponential cloud structure.

In the text

Fig. 4 Reconstructed exponential particle size distributions (Eq. (1)) from the Y-dwarf KCl case moment results (Left: log g = 3.25, Right: log g = 4.25) in Section 6 when the exponential distribution is assumed. The colour bar shows the atmospheric pressure of the distribution. This shows the limitations of the exponential distribution as limited flexibility is present due to the assumed variance, resulting in a highly similar distribution shape across the cloud layer.

In the text