{'fraction_non_alphanumeric': 0.03660908362887542, 'fraction_numerical': 0.021469702169850894, 'mean_word_length': 4.399670544630907, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Collision fluctuations of lucky droplets with superdroplets

It was previously shown that the superdroplet algorithm for modeling the collision-coalescence process can faithfully represent mean droplet growth in turbulent clouds. But an open question is how accurately the superdroplet algorithm accounts for fluctuations in the collisional aggregation process. Such fluctuations are particularly important in dilute suspensions. Even in the absence of turbulence, Poisson fluctuations of collision times in dilute suspensions may result in substantial variations in the growth process, resulting in a broad distribution of growth times to reach a certain droplet size. We quantify the accuracy of the superdroplet algorithm in describing the fluctuating growth history of a larger droplet that settles under the effect of gravity in a quiescent fluid and collides with a dilute suspension of smaller droplets that were initially randomly distributed in space ('lucky droplet model'). We assess the effect of fluctuations upon the growth history of the lucky droplet and compute the distribution of cumulative collision times. The latter is shown to be sensitive enough to detect the subtle increase of fluctuations associated with collisions between multiple lucky droplets. The superdroplet algorithm incorporates fluctuations in two distinct ways: through the random spatial distribution of superdroplets and through the Monte Carlo collision algorithm involved. Using specifically designed numerical experiments, we show that both on their own give an accurate representation of fluctuations. We conclude that the superdroplet algorithm can faithfully represent fluctuations in the coagulation of droplets driven by gravity.


Introduction
Direct numerical simulations (DNS) have become an essential tool to investigate collisional growth of droplets in turbulence (Onishi et al. 2015;Saito and Gotoh 2018). Here, DNS refers to the realistic modeling of all relevant processes, which involves not only the use of a realistic viscosity, but also a realistic modeling of collisions of droplet pairs in phase space. The most natural and physical way to analyze collisional growth is to track individual droplets and to record their collisions, one by one. However, DNS of the collision-coalescence process are very challenging, not only when a large number of droplets must be tracked, but also because the flow must be resolved over a large range of time and length scales.
Over the past few decades, an alternative way of modeling aerosols has gained popularity. Zannetti (1984) introduced the concept of "superparticles, i.e., simulation particles representing a cloud of physical particles having similar characteristics." This concept was also used by Paoli et al. (2004) in the context of condensation problems. The application to coagulation problems was pioneered by Zsom and Dullemond (2008) and Shima et al. (2009), who also developed a computationally efficient algorithm.
Less is known about how the superdroplet algorithm represents fluctuations in the collisional aggregation process. Dziekan and Pawlowska (2017) compared the results of the superdroplet algorithm with the predictions of the stochastic coagulation equation of Gillespie (1972) in the context of coalescence of droplets settling in a quiescent fluid. Dziekan and Pawlowska (2017) concluded that the results of the superdroplet algorithm qualitatively agree with what Kostinski and Shaw (2005) called the lucky droplet model (LDM). To assess the importance of fluctuations, Dziekan and Pawlowska (2017) computed the time 10% , after which 10% of the droplets have reached a radius of 40 m. In agreement with earlier Lagrangian simulations of Onishi et al. (2015), which did not employ the superdroplet algorithm, they found that the difference in 10% between their superdroplet simulations and the stochastic model of (Gillespie 1972) decreases with the square root of the number of droplets, provided that there are no more than about nine droplets per superdroplet. The number of droplets in each superdroplet is called the multiplicity ( ). When this number is larger than 9, they found that a residual error remains. We return to this question in the discussion of the present study, where we tentatively associate their findings with the occurrence of several large (lucky) droplets that grew from the finite tail of their initial droplet distribution.
The role of fluctuations is particularly important in dilute systems, where rare extreme events may substantially broaden the droplet-size distribution. This is well captured by the LDM, which was first proposed by Telford (1955) and later numerically addressed by Twomey (1964), and more recently quantitatively analyzed by Kostinski and Shaw (2005). The model describes one droplet of 12.6 m radius settling through a dilute suspension of background droplets with 10 m radius. The collision times between the larger ("lucky") droplet and the smaller ones are exponentially distributed, leading to substantial fluctuations in the growth history of the lucky droplet. Wilkinson (2016) derived analytic expressions for the cumulative distribution times using largedeviation theory. Madival (2018) extended the theory of Kostinski and Shaw (2005) by considering a more general form of the droplet-size distribution than just the Poisson distribution.
The goal of the present study is to investigate how accurately the superdroplet algorithm represents fluctuations in the collisional growth history of settling droplets in a quiescent fluid. Unlike the work of Dziekan and Pawlowska (2017), who focused on the calculation of 10% , we compare here with the distribution of cumulative collision times, which is the key diagnostics of the LDM. We record growth histories of the larger droplet in an ensemble of different realizations of identical smaller droplets that were initially randomly distributed in a quiescent fluid. We show that the superdroplet algorithm accurately describes the fluctuations of growth histories of the lucky droplet in an ensemble of simulations. In its simplest form, the LDM assumes that the lucky droplet is large compared to the background droplets, so that the radius of those smaller droplets can be neglected in the geometrical collision cross section and velocities of colliding droplets; see Eqs.
Since fluctuations early on in the growth history are most important (Kostinski and Shaw 2005;Wilkinson 2016), this can make a certain difference in the distribution of the time it takes for the lucky droplet to grow to a certain size. As the small droplets are initially randomly distributed, their local number density fluctuates. Consequently, lucky droplets can grow most quickly where the local number density of small droplets happens to be large.
The remainder of this study is organized as follows. In section 2 we describe the superdroplet algorithm and highlight differences between different implementations used in the literature (Shima et al. 2009;Johansen et al. 2012;Li et al. 2017). Section 3 summarizes the LDM, the setup of our superdroplet simulations, and how we measure fluctuations of growth histories. Section 4 summarizes the results of our superdroplet simulations. We conclude in section 6.

a. Superdroplet algorithm
Superdroplet algorithms represent several physical droplets by one superdroplet. All droplets in superdroplet are assumed to have the same material density d , the same radius , the same velocity v , and reside in a volume around the same position x . The index labeling the superdroplets ranges from 1 to s ( 0 ) ( Table 1), where 0 denotes the initial time.
The equation of motion for the position x and velocity v of superdroplet reads:  Here g is the gravitational acceleration, and the hydrodynamic force is modeled using Stokes law, so that is the droplet response (or Stokes) time attributed to the superdroplet, = 10 −5 m 2 s −1 is the viscosity of air, and is the mass density of the airflow. Droplets are only subject to gravity and no turbulent airflow is simulated. Droplet collisions are represented by collisions of superdroplets (Shima et al. 2009;Johansen et al. 2012;Li et al. 2017), as mentioned above. Superdroplets and (collision partners) residing inside a grid cell collide with probability = , where is the integration time step and is their collision rate.
A collision happens when < , where 0 ≤ ≤ 1 is a uniformly distributed random number. To avoid a probability larger than unity, we limit the integration step through the condition ≪ 1/ (Johansen et al. 2012;Li et al. 2017). The collision rate is given by where is the collision efficiency, max = max( , ) is the larger one of the multiplicities and of superdroplets or (Table 1), and 3 is the volume of the grid cell closest to the superdroplet. The number density of physical droplets in superdroplet is then = / 3 . Note that Eq. (4) implies that droplets having the same velocity (v = v ) never collide. This also implies that no collisions are possible between physical particles within a single superdroplet. For the purpose of the present study, it suffices to limit ourselves to the simplest, albeit unrealistic assumption of = 1, but we also consider in one case a slightly more realistic quadratic dependence on the radius of the larger droplet. To assess the effects of this assumption, we compare with results where the efficiency increases with droplet radius (Lamb and Verlinde 2011). Following Kostinski and Shaw (2005) and Wilkinson (2016), we adopt a simple power law prescription for the dependence of the efficiency on the droplet radius.
What happens when two superdroplets collide? The collision scheme suggested by Shima et al. (2009) amounts to the following rules; see also Fig. 1 for an illustration. To ensure mass conservation between superdroplets and , when > , which is the case illustrated in Fig. 1(b), droplet numbers and masses are updated such that where and are the droplet masses. When < , which is the case shown in Fig. 1(a), the update rule is also given by Eq. (5), but with indices and exchanged. In other words, the number of droplets in the smaller superdroplet remains unchanged (and their masses are increased), while that in the larger one is reduced by the amount of droplets that have collided with all the droplets of the smaller superdroplet (and their masses remain unchanged).
To ensure momentum conservation during the collision, the momenta of droplets in the two superdroplets are up- after a collision of superdroplets. Finally, when = , which is the case described in Fig. 1(c), droplet numbers and masses are updated as and it is then assumed that, when two superdroplets, each with one or less than one physical droplet, collide, the superdroplet containing the smaller physical droplet is collected by the more massive one; it is thus removed from the computational domain after the collision, still conserving mass and momentum. We emphasize that Eq. (5) does not require to be an integer. Since we usually specify the initial number density of physical particles, can be fractional from the beginning. This is different from the integer treatment of in Shima et al. (2009). The superdroplet simulations are performed by using the particle modules of the P C (Pencil Code Collaboration et al. 2021). The fluid dynamics modules of the code are not utilized here. To reduce the computational cost and make it linear in the number of superdroplets per mesh point, s ( ), Shima et al. (2009) supposed that each superdroplet interacts with only one randomly selected superdroplet per time step rather than allowing collisions with all the other superdroplets in a grid cell (they still allow multiple coalescence for randomly generated, non-overlapping candidate pairs in one time step, which is what they referred to as random permutation technique. This technique was also adopted by Dziekan and Pawlowska (2017) and Unterstrasser et al. (2020). However, this is not used in the P C . Instead, we allow each superdroplet to collide with all other superdroplets within one grid cell to maximize the statistical accuracy of the results. This leads to a computational cost of O( 2 s ( )), which does not significantly increase the computational cost because s ( ) is relatively small for cloud-droplet collision simulations. In the P C , collisions between particles residing within a given grid cell are evaluated by the same processor which is also evaluating the equations of that grid cell. Due to this, together with the domain decomposition used in the code, the particle collisions are automatically efficiently parallelized as long as the particles are more or less uniformly distributed over the domain.

b. Numerical setup
In our superdroplet simulations, we consider droplets of radius 10 m, randomly distributed in space, together with one droplet of twice the mass, so that the radius is 2 1/3 × 10 m = 12.6 m. The larger droplet has a higher settling speed than the 10 m droplets and sweeps them up through collision and coalescence. For each simulation, we track the growth history of the larger droplet until it reaches 50 m in radius and record the time it takes to grow to that size.
In the superdroplet algorithm, one usually takes ( 0 ) ≫ 1, which implies that the actual number of lucky droplets is also more than one. This was not intended in the original formulation of the lucky droplet model (Telford 1955;Kostinski and Shaw 2005;Wilkinson 2016) and could allow the number of superdroplets with heavier (lucky) droplets, (luck) s , to become larger than unity. This would manifest itself in the growth history of the lucky droplets through an increase by more than the mass of a background droplet. We refer to this as "jumps". Let us therefore now discuss the conditions under which this would happen and denote the values of ( 0 ) for the lucky and background droplets by luck and back , respectively. First, for luck = back , the masses of both lucky and background superdroplets can increase, provided their values of ( 0 ) are above unity; see Fig. 1(c). Second, even if luck < back initially, new lucky superdroplets could in principle emerge when the same two superdroplets collide with each other multiple times. This can happen for two reasons. First, the use of periodic boundary conditions for the superdroplets (i.e., in the vertical direction in our laminar model with gravity). Second, two superdroplets can remain at the same location (corresponding to the same mesh point of the Eulerian grid for the fluid) during subsequent time steps. The simulation time step must be less than both the time for a superdroplet to cross one grid spacing and the mean collision time, i.e., the inverse collision rate given by Eq. (4). Looking at Fig. 1, we see that back can then decrease after each collision and potentially become equal to or drop below the value of luck . This becomes exceedingly unlikely if initially back ≫ luck , but it is not completely impossible, unless luck is chosen initially to be unity.
The initial value of back can in principle also be chosen to be unity. Although such a case will indeed be considered here, it would defeat the purpose and computational advantage of the superdroplet algorithm. Therefore, we also consider the case back ≫ luck . As already mentioned, jumps are impossible if luck is unity. For orientation, we note that the speed of the lucky droplet prior to the first collision is about 3.5 cm s −1 , the average time to the first collision is 490 s, and thus, it falls over a distance of about 17 m before it collides.
The superdroplet algorithm is usually applied to threedimensional (3-D) simulations. If there is no horizontal mixing, one can consider one-dimensional (1-D) simulations. Moreover, we are only interested in the column in which the lucky droplet resides. In 3-D, however, the number density of the 10 m droplets beneath the lucky one is in general not the same as the mean number density of the whole domain. This leads to yet another element of randomness: fluctuations of the number density between columns.
Equation (1) is solved with periodic boundary conditions using the P C (Pencil Code Collaboration et al. 2021), which employs a third-order Runge-Kutta time stepping scheme. The superdroplet algorithm is implemented in the P C , which is used to solve equations (3)-(7). For the 1-D superdroplet simulations, we employ an initial number density of background droplets of 0 ≈ 3 × 10 8 m −3 within a volume = × × with = = 0.002 m, = 0.214 m, and s ( 0 ) = 256 such that the multiplicity is luck ( 0 ) = back ( 0 ) = 1. For each simulation, 7,686,000 time steps are integrated with an adaptive time step with a mean value of = 2.942×10 −4 s. For a superdroplet with an initial radius of 12.6 m to grow to 50 m, 123 collisions are required. For the purpose of the present study, we designed a parallel technique to run thousands of 1-D superdroplet simulations simultaneously (see details in appendix A1).

a. Basic idea
The LDM describes the collisional growth of a larger droplet that settles through a quiescent fluid and collides with smaller monodisperse droplets, that were initially randomly distributed in space. This corresponds to the setup described in the previous section. We begin by recalling the main conclusions of Kostinski and Shaw (2005). Initially, the lucky droplet has a radius corresponding to a volume twice that of the background droplets, whose radius was assumed to be 1 = 10 m. Therefore, its initial radius is 2 = 2 1/3 1 = 12.6 m. After the ( − 1)th collision step with smaller droplets, it increases as Fluctuations in the length of the time intervals between collision −1 and give rise to fluctuating growth histories of the larger droplet. These fluctuations are quantified by the distribution of the cumulative time corresponding to 123 collisions needed for the lucky droplet to grow from 12.6 m to 50.0 m (note that Kostinski and Shaw (2005) used one more collision, so their final radius was actually 50.1 m). The time intervals between successive collisions are drawn from an exponential distribution with a probability ( ) = exp(− ). The rates depend on the differential settling velocity |v − v 1 | between the colliding droplets through Eqs. (3) and (4). Here, however, the background droplets have always the radius 1 , so the collision rate at the ( − 1)th collision of the lucky droplet with radius obeys where = ( , 1 ), and v and v 1 are approximated by their terminal velocities.
While the LDM is well suited for addressing theoretical questions regarding the significance of rare events, it should be emphasized that it is at the same time highly idealized. Furthermore, while it is well known that ≪ 1 (Pruppacher and Klett 1997), it is instructive to assume, as an idealization, = 1 for all , so the collision rate (10) can be approximated as ∼ 4 (Kostinski and Shaw 2005), which is permissible when ≫ 1 . It follows that, in terms of the collision index , the collision frequency is where * = (2 /9)( d / )( / ) 4 1 , and is the number density of the 10 m background droplets. This is essentially the model of Kostinski and Shaw (2005) and Wilkinson (2016), except that they also assumed ≠ 1. They pointed out that, early on, i.e., for small , is small and therefore the mean collision time −1 is long. We note that the variance of the mean collision time is −2 , which is large for small . The actual time until the first collision can be very long, but it can also be very short, depending on fluctuations. Therefore, at early times, fluctuations have a large impact on the cumulative collision time. Note that for droplets with ≥ 30 m, the linear Stokes drag is not valid (Pruppacher and Klett 1997).

b. Relaxing the power law approximation
We now discuss the significance of the various approximations being employed in the mathematical formulation of the LDM of Kostinski and Shaw (2005). To relax the approximations made in Eq. (11), we now write it in the form would correspond to the expression Eq. (10) used in the superdroplet algorithm. In Eq. (11), however, it was assumed that A = B = . To distinguish this approximation from the form used in Eq. (12), we denote that case by writing symbolically " A ≠ ≠ B "; see Fig. 2. In Eq. (13), we have introduced A and B to study the effect of relaxing the assumption A = B = , made in simplifying implementations of the LDM. Both of these assumptions are justified at late times when the lucky droplet has become large compared to the smaller ones, but not early on, when the size difference is moderate. By comparison, if fluctuations are ignored, the collision times are given by = −1 . This is what we refer to as mean-field theory (MFT). In Fig. 3 we demonstrate the effect of the contributions from A and B on the mean cumulative collision time in the corresponding MFT, where are the inverse of the mean collision rates. We see that, while the contribution from A shortens the mean collision time, that of B enhances it. In Fig. 2, we also see that the F . 4. Comparison of ( ) in a double-logarithmic representation for the LDM appropriate to our benchmark (black solid line) with various approximations where A = B = (red dashed line) along with a case where only B = is assumed (blue dash-dotted line). Here we used 10 10 realizations. Note that we plot the distribution of the cumulative times versus the normalized time, / , as was done in the work of Kostinski and Shaw (2005). Normalizing by allows us to see changes in the shape of ( / ), thus allows a more direct comparison of the subtle differences in the shapes of the different curves and ensures that the peaks of all curves are at approximately the same position.
contributions to the two correction factors 2 / 2 A and 2 / 2 B have opposite trends, which leads to partial cancelation in their product.
In Fig. 4 we show a comparison of the distribution of cumulative collision times for various representations of . Those are computed numerically using 10 10 realizations of sequences of random collision times . We refer to appendix A1 for details of performing this many realizations.
The physically correct model is where A ≠ ≠ B (black line in Fig. 4). To demonstrate the sensitivity of ( ) to changes in the representation of , we show the result for the approximations A = = B (red line) and A ≠ = B (blue line). The ( ) curve is also sensitive to changes in the collision efficiency late in the evolution. To demonstrate this, we assume ∝ 2 when exceeds a certain arbitrarily chosen value * between 10 and 40 m, and = const below * (Lamb and Verlinde 2011). To ensure that ≤ 1, we take with * = ( * /50 m) 2 . However, the normalized ( ) curves are independent of the choice of the value of * . In Fig. 5, we show the results for A ≠ ≠ B using * = 40 m and 30 m (red and blue lines, respectively) and compare with the case = const. The more extreme cases with * = 20 m and 10 m are shown as gray lines. The latter is similar to the case ∼ 6 considered by Kostinski and Shaw (2005) and Wilkinson (2016). To quantify the shape of ( ), we give in Table 2 the average of ≡ ln( / ), its standard deviation = 2 1/2 , where ≡ − , its skewness skew = 3 / 3 , and its kurtosis kurt = 4 / 4 − 3. We recall that, for a perfectly lognormal distribution, skew = kurt = 0. The largest departure from zero is seen in the skewness, which is positive, indicating that the distribution broadens for large . The kurtosis is rather small, however.
The main conclusion that can be drawn form the investigation mentioned above is that, as far as the shapes of the different curves are concerned, it does not result in any significant error to assume ≫ 1 . The value of is only about 10% smaller if A = = B is used (compare the red dashed and black solid lines in Fig. 4). This is because the two inaccuracies introduced by A and B almost cancel each other. When * = 40 m or 30 m, for example, the values of increase by 3% and 15%, respectively; see Table 3, where we also list the corresponding values of MFT 124 .
On the other hand, the actual averages such as ≈ MFT 124 vary by almost 50%. A straightforward extension of the LDM is to take horizontal variations in the local column density into account. Those are always present for any random initial conditions, but could be larger for turbulent systems, regardless of the droplet speeds. In 3-D superdroplet simulations, large droplets can fall in different vertical columns that contain different numbers of small droplets, a consequence of the fact that the small droplets are initially randomly distributed. To quantify the effect of varying droplet number densities in space, it is necessary to solve for an ensemble of columns with different number densities of the 10 m background droplets and compute the distribution of cumulative collision times. These variations lead to a broadening of ( ), but it is a priori not evident how important this effect is. A quantitative analysis is given in appendix A3.

c. Relation to the superdroplet algorithm
To understand the nature of the superdroplet algorithm, and why it captures the lucky droplet problem accurately, it is important to realize that the superdroplet algorithm is actually a combination of two separate approaches to solving the LDM, each of which turns out to be able to reproduce the lucky droplet problem to high precision. In principle, we can distinguish four different approaches (Table 3) to obtaining the collision time interval . In approach I, was taken from an exponential distribution of random numbers. Another approach is to use a randomly distributed set of 10 m background droplets in space and then determine the distance to the next droplet within a vertical cylinder of possible collision partners to find the collision time (approach II). A third approach is to use the mean collision rate to compute the probability of a collision within a fixed time interval. We then use a random number between zero and one (referred to as Monte Carlo method; see, e.g., Sokal 1997) to decide whether at any time there is a collision or not (approach III). This is actually what is done within each grid cell in the superdroplet algorithm; see Eqs. (3) and (4). The fourth approach is the superdroplet algorithm discussed extensively in section 2.a (approach IV). It is essentially a combination of approaches II and III. We have compared all four approaches and found that they all give very similar results. In the following, we describe approaches II and III in more detail, before focussing on approach IV in section 4.

d. Solving for the collisions explicitly
A more realistic method (approach II; see Table 3) is to compute random realizations of droplet positions in a tall box of size 2 h × , where h and are the horizontal and vertical extents, respectively. We position the lucky droplet in the middle of the top plane of the box. Collisions  are only possible within a vertical cylinder of radius + 1 below the lucky droplet. Next, we calculate the distance Δ to the first collision partner within the cylinder. We assume that both droplets reach their terminal velocity well before the collision. This is an excellent approximation for dilute systems such as clouds, because the droplet response time of Eq. (2) is much shorter than the mean collision time. Here we use the subscript to represent the time until the ( − 1)th collision, which is equivalent to the th droplet. We can then assume the relative velocity between the two as given by the difference of their terminal velocities as The time until the first collision is then given by 2 = Δ /Δ 2 . This collision results in the lucky droplet having increased its volume by that of the 10 m droplet. Correspondingly, the radius of the vertical cylinder of collision partners is also increased. We then search for the next collision partner beneath the position of the first collision, using still the original realization of 10 m droplets. We continue this procedure until the lucky droplet reaches a radius of 50 m. Approach II is an explicit method compared to other approaches listed in Table 3.

e. The Monte Carlo method to compute
In the Monte Carlo method (approach III; see Table 3) we choose a time step and step forward in time. As in the superdroplet algorithm, the probability of a collision is given by = ; see Eq.
(3). We continue until a radius of 50 m is reached. We note that in this approach, is kept constant, i.e., no background droplet is being removed after a collision.
Approach III also allows us to study the effects of jumps in the droplet size by allowing for several lucky droplets at the same time and specifying their collision probability appropriately. These will then be able to interact not only with the 10 m background droplets, but they can also collide among themselves, which causes the jumps. We will include this effect in solutions of the LDM using approach III and compare with the results of the superdroplet algorithm.

a. Accuracy of the superdroplet algorithm
We now want to determine to what extent the fluctuations are correctly represented by the superdroplet algorithm. For this purpose, we now demonstrate the degree of quantitative agreement between approaches I-III and the corresponding solution with the superdroplet algorithm (approach IV; see Table 3). This is done by tracking the growth history of each lucky droplet. As the first few collisions determine the course of the formation of larger droplets, we also use the distribution ( ) of cumulative collision times . We perform real superdroplet simulations with different random seeds using ( 0 ) = 1.
We begin by looking at growth histories for many individual realizations obtained from the superdroplet simulation. Fig. 6 shows an ensemble of growth histories (thin gray lines) obtained from real ≈ 10 3 independent simulations, as described above. The times between collisions are random, leading to a distribution of cumulative growth times to reach 50 m. Also shown is the mean growth curve (thick black line), obtained by averaging the time at fixed radii . This figure demonstrates that the fluctuations are substantial. We also see that large fluctuations relative to the average time are rare.
To quantify the effect of fluctuations from all realizations, we now consider the corresponding ( ) in Fig. 7. We recall that ( 0 ) = 1 for our superdroplet simulation in Fig. 7. However, a simulation with ( 0 ) = 50 yields almost the same result; see appendix A2.
The comparison of the results for the LDM using approach I and the superdroplet algorithm shows small dif- ferences. The width of the ( ) curve is slightly larger for approach I than for the superdroplet simulations. This suggests that the fluctuations, which are at the heart of the LDM, are slightly underrepresented in the superdroplet algorithm. However, this shortcoming may also be a consequence of our choice of having used only 256 superdroplets, i.e., one lucky and 255 background superdroplets. Given that the multiplicities of lucky and background droplets was unity, each collision removed one background droplet. Thus, after 123 collisions, almost 50% of the background droplets were removed by the time the lucky droplet reached 50 m. Nevertheless, as we will see below, this has only a small effect.
An important question is to what extent our results depend on the number density of background droplets and the size of the computational domain. To examine this with the superdroplet algorithm (approach IV), we consider three values of the initial number density: 0 = 3 × 10 8 m −3 , T 4. Runs of Fig. 8(a), where and ,luck are in units of 0 = 7.5 × 10 7 m −3 , is in units of 0 = 3 × 10 8 m −3 , luck,0 is in units of 1.2 × 10 5 m −3 , using 3 = 2.6 × 10 −8 m 3 for all runs, except for the last one (gray), where it is a factor 2 smaller. color in Fig. 8 F . 8. (a): ( ) for 0 (red), 10 0 (blue), and 100 0 (black) with 0 = 3 × 10 8 m −3 and = 0.214 m. In the last case, 0.5% of the background droplets were removed; the orange symbols denote a case with 100 times larger value of luck , where 50% of the background droplets were removed. The gray symbols denote the case of Fig. 7 with ten times more physical lucky droplets (10 luck,0 ). (b): ( / ) for = 0 (red), 2 0 (green), 10 0 (blue), 100 0 with 10 luck (black), 1000 luck (orange), and 0 with 10 luck,0 (gray). (c): ( / ) for , 2 , 8 , and 64 with 100 0 , obtained using the superdroplet algorithm (approach IV). The red dash-dotted line in (b) represents the LDM (approach I) with A ≠ ≠ B and 0 = 3 × 10 8 m −3 , which is the same simulation as the one in Fig. 7. The green dots in (b) is for 8192 realizations, while all the other simulations are for 1024 realizations. 10 0 , and 100 0 , while the initial number density of the lucky droplet is luck,0 ≡ 1.2×10 5 m −3 , 10 luck,0 , and again 10 luck,0 , respectively; see Table 4 for a summary. Thus, even though the lucky droplet has to collide 123 times to reach 50 m, it only removes 123 luck / 0 = 5%, 5%, and 0.5% of the droplets, respectively. Fig. 8 shows ( ) for these three cases using first the cumulative time [ Fig. 8(a)] and then the normalized time / [ Fig. 8(b)]. We see that the positions of the peaks in ( ) change linearly with the initial number density 0 , but ( / ) are very similar to each other. This is related to the fact that, after normalization, 0 drops out from the expression for / in the LDM (approach I); see Eq. (9). At small values of / , however, all curves show a similar slight underrepresentation of the fluctuations as already seen in Fig. 7. In all these simulations, we used 1024 realizations, except possibly for one case where we used 8192 realizations; see the green symbols in Fig. 8(b). The distribution of cumulative growth times is obviously much smoother in the latter case, but the overall shape is rather similar.
In the above, the number density of the lucky droplets has been much smaller than the number density of the background droplets. This means that for each collision the physical number of background droplets changed by only a small amount (5% or 0.5%). To see how sensitive our results for ( ) are to this number, we now perform an extra experiment where 50% of the background droplets are removed by the time the lucky droplet reaches 50 m. This is also shown in Fig. 8(a) and (b); see the orange symbols, where luck = back = 200. We see that even for 50% removal the results are essentially unchanged.
In our superdroplet simulations (approach IV; see Table 3), the vertical extent of the simulation domain is only = 0.214 m. This is permissible given that we use periodic boundary conditions for the particles. Nevertheless, the accuracy of our results may suffer from poor statistics. To investigate this in more detail, we now perform 1-D simulations with 2 , 8 , and 64 . At the same time, we increased the number of mesh points and the number of superdroplets by the same factors. Since the shape of ( /¯ ) is almost independent of 0 , as shown in Fig. 8(b), we use 0 = 3 × 10 10 m −3 instead of 0 = 3 × 10 8 m −3 to reduce the computational cost. As shown in Fig. 8(c), ( /¯ ) is insensitive to the domain size. Therefore, our results with = 0.214 m can be considered as accurate with respect to ( /¯ ).
In the following, we discuss how our conclusions relate to those of earlier work. We then discuss a number of additional factors that can modify the results. Those additional factors can also be taken into account in the LDM. Even in those cases, it turns out that the differences between the LDM and the superdroplet algorithm are small.

b. The occurrence of jumps
One of the pronounced features in our superdroplet simulations with ( 0 ) > 1 is the possibility of jumps. We see examples in Fig. 9 where luck = back = 2 and the jumps are visualized by the red vertical lines. Those jumps are caused by the coagulation of the lucky droplet with droplets of radii larger than 10 m that were the result of other lucky droplets in the simulations. What is the effect of these jumps? Could they be responsible for the behavior found by Dziekan and Pawlowska (2017) that the difference in their 10% between the numerical and theoretical F . 9. Same as Fig. 6 but with initial condition ( 0 ) = 2 using s ( 0 ) = 128, corresponding to the same number of physical droplets as in Fig. 6, where ( 0 ) = 1. Note the occurrence of jumps, indicated in red.
calculation decreases with the square root of the number of physical droplets, as we discussed in section 1?
It is clear that those jumps occur mainly during the last few steps of a lucky droplet growing to 50 m (see Fig. 9) when there has been enough time to grow several more lucky droplets. Because the collision times are so short at late times, the jumps are expected to be almost insignificant. To quantify this, it is convenient to use approach III, where we choose (luck) s = 3 superdroplets simultaneously. (As always in approach III, the background particles are still represented by only one superdroplet, and is kept constant.) We also choose luck = 1, and therefore The lucky droplets can grow through collisions with the 10 m background droplets and through mutual collisions between lucky droplets. The collision rate between lucky droplets and is, analogously to Eq. (12), given by where luck is the number density of physical droplets in the superdroplet representing the lucky droplet. To obtain an expression for luck in terms of the volume of a grid cell 3 , we write luck = luck / 3 . The ratio of the physical number of lucky droplets, To investigate the effect of jumps on ( ) in the full superdroplet model studied above (see Figs. 6 and 9), F . 10. Comparison of models with = 0 (no jumps), 0.004 (the value expected for the simulations), 0.02, and 0.05 using approach III; see Table 3. we first consider the case depicted in Fig. 6, where luck = back ≡ ( 0 ) = 1. Here, we used s = 256 superdroplets, of which one contained the lucky droplet, so (luck) s = 1, and the other 255 superdroplets contained a 10 m background droplet each. In our superdroplet solution, the ratio (19) was therefore ≈ 1/255 = 0.004. Using approach III, enters simply as an extra factor in the collision probability between different lucky droplets. (In approach III, all quantities in Eq. (19) are kept constant.) The effect on ( ) is shown in Fig. 10, where we present the cumulative collision times for models with three values of using approach III. We see that for small values of , the cumulative distribution function is independent of , and the effect of jumps is therefore negligible (compare the black solid and the red dashed lines of Fig. 10). More significant departures due to jumps can be seen when = 0.02 and larger.
Let us now compare with the case in which we found jumps using the full superdroplet approach (approach IV). The jumps in the growth histories cause the droplets to grow faster than without jumps. However, jumps do not have a noticeable effect upon ( ) in the superdroplet simulations we conducted; see Fig. 11. By comparing ( ) for back = 40 (blue crosses in Fig. 11) with that for back = 2 (black circles), while keeping luck = 2 in both cases, hardly any jumps occur and the lucky droplet result remains equally accurate.
For larger values of , jumps occur much earlier, as can be seen from Fig. 12, where we show 30 growth curves for the cases = 0.004, which is relevant to the simulations of Fig. 7, as well as = 0.02, and 0.05. We also see that for large values of , the width in the distribution of arrival times is broader and that both shorter and longer times are possible. This suggests that the reason for the finite residual error in the values of 10% found by Dziekan and Pawlowska (2017) for (  ( / ) of simulations in Fig. 9 (black circles) and the ones with initially back = 40 (blue crosses). luck = 2 in both cases. The red line denotes the LDM (approach I) with A ≠ ≠ B , which is the same simulation as the one in Fig. 7. T

5.
Comparison of the moments of = ln( / ) for approaches I-III. be due to jumps. In our superdroplet simulations, by contrast, jumps cannot occur when ( 0 ) = 1 or back ≫ luck .

c. The two aspects of randomness
Let us now quantify the departure that is caused by the use of the Monte Carlo collision scheme. To do this, we need to assess the effects of randomness introduced through Eqs. (3) and (4) on the one hand and the random distribution of the 10 m background droplets on the other. Both aspects enter in the superdroplet algorithm.
We recall that in approach II, fluctuations originate solely from the random distribution of the 10 m background droplets. In approach III, on the other hand, fluctuations originate solely from the Monte Carlo collision scheme. By contrast, approach I is different from either of the two, because it just uses the exponential distribution of the collision time intervals, which is indirectly reproduced by the random initial droplet distribution in approach II and by the Monte Carlo scheme in approach III.
In Fig. 13, we compare approaches I, II, and III. For our solution using approach II, we use a nonperiodic domain of size 10 −4 × 10 −4 × 700 m 3 , thus containing on average 2100 droplets. This was tall enough for the lucky droplet to reach 50 m for all the 10 7 realizations in this experiment. The F . 13. Comparison of ( ) for approaches I, II, and III. differences between them are very minor, and also the first few moments are essentially the same; see Table 5. We thus see good agreement between the different approaches. This suggests that the fluctuations introduced through random droplet positions is not crucial and that it can be substituted by the fluctuations of the Monte Carlo scheme alone.
It is worth noting that we were able to perform 10 7 and 10 6 realizations with approaches II and III, respectively, and 10 10 realizations with approach I, while in the superdroplet algorithm (approach IV), we could only run 10 3 -10 4 realizations due to the limitation of the computational power. This may be the reason why fluctuations appear to be slightly underrepresented in the superdroplet algorithm; see Fig. 7 and the discussion in section 4.a. Nevertheless, the agreement between the LDM and the superdroplet simulations demonstrates that the superdroplet algorithm is able to represent fluctuations during collisions and does not contain mean-field elements. This can be further evidenced by the fact that the results of approaches II and III agree perfectly with those of approach I, and the superdroplet algorithm is just the combination of approaches II and III.

Discussion
Fluctuations play a central role in the LDM. We have therefore used it as a benchmark for our simulation. It turns out that the superdroplet algorithm is able to reproduce the growth histories qualitatively and the distribution of cumulative collision times quantitatively. The role of fluctuations was also investigated by Dziekan and Pawlowska (2017), whose approach to assessing the fluctuations is different from ours. Instead of analyzing the distribution of cumulative collision times, as we do here, their primary diagnostics is the time 10% , after which 10% of the mass of cloud droplets has reached a radius of 40 m. In the LDM, such a time would be infinite, because there is only one droplet that is allowed to grow. They then determined the accuracy with which the value of 10% is determined. The accuracy increases with the square root of the number of physical droplets, provided that the ratio ( 0 ) is kept below a limiting value of about 9. For ( 0 ) > 9, they found that there is always a residual error in the value of 10% that no longer diminishes as they increase the number of physical droplets. We have demonstrated that, when ( 0 ) > 1, jumps in the growth history tend to occur. Those jumps can lead to shorter cumulative collision times, which could be the source of the residual error they find.
For a given fraction of droplets that first reach a size of 40 m, they also determined their average cumulative collision time. They found a significant dependence on the number of physical droplets. This is very different in our case where we just have to make sure that the number of superdroplets is large enough to keep finding collision partners in the simulations. However, as the authors point out, this is a consequence of choosing an initial distribution of droplet sizes that has a finite width. This implies that for a larger number of droplets, there is a larger chance that there could be a droplet that is more lucky than for a model with a smaller number of droplets. In our case, by contrast, we always have a well-known number of superdroplets of exactly 12.6 m, which avoids the sensitivity on the number of droplets.
The ( 0 ) = 9 limit of Dziekan and Pawlowska (2017) does not hold in this investigation. In this context we need to recall that their criterion for acceptable quality concerned the relative error of the time in which 10% of the total water has been converted to 40 m droplets. In our case, we have focussed on the shape of the ( ) curve, especially for small .

Conclusions
We investigated the growth histories of droplets settling in quiescent air using superdroplet simulations. The goal was to determine how accurately these simulations represent the fluctuations of the growth histories. This is important because the observed formation time of drizzle-sized droplets is much shorter than the one predicted based on the mean collisional cross section. The works of Telford (1955), Kostinski and Shaw (2005), and Wilkinson (2016) have shown that this discrepancy can be explained by the presence of stochastic fluctuations in the time intervals between droplet collisions. By comparing with the lucky droplet model (LDM) quantitatively, we have shown that the superdroplet simulations capture the effect of fluctuations.
A tool to quantify the significance of fluctuations on the growth history of droplets is the distribution of cumulative collision times. Our results show that the superdroplet algorithm reproduces the distribution of cumulative collision times that is theoretically expected based on the LDM. However, the approximation of representing the dependence of the mean collision rate on the droplet radius by a power law is not accurate and must be relaxed for a useful benchmark experiment.
In summary, the superdroplet algorithm appears to take fluctuations fully into account, at least for the problem of coagulation due to gravitational settling in quiescent air. Computing the distribution of cumulative collision times in the context of turbulent coagulation would be rather expensive, because one would need to perform many hundreds of fully resolved 3-D simulations. Our study suggests that fluctuations are correctly described for collisions between droplets settling in quiescent fluid, but we do not know whether this conclusion carries over to the turbulent case.
Acknowledgments. This work was supported through the FRINATEK grant 231444 under the Research Council of Norway, SeRC, the Swedish Research Council grants 2012-5797, 2013-03992, and 2017-03865, Formas grant 2014-585, by the University of Colorado through its support of the George Ellery Hale visiting faculty appointment, and by the grant "Bottlenecks for particle growth in turbulent aerosols" from the Knut and Alice Wallenberg Foundation, Dnr. KAW 2014.0048. The simulations were performed using resources provided by the Swedish National Infrastructure for Computing (SNIC) at the Royal Institute of Technology in Stockholm and Chalmers Centre for Computational Science and Engineering (C3SE). This work also benefited from computer resources made available through the Norwegian NOTUR program, under award NN9405K.
Data availability statement. The source code used for the simulations of this study, the P C , is freely available on https://github.com/pencil-code/. Datasets for "Collision fluctuations of lucky droplets with superdroplets" (v2021.05.07) are available under https://doi.org/10.5281/zenodo.4742786; see also http://www.nordita.org/~brandenb/projects/lucky/ for easier access. The plotting and analysis scripts are also included. Some of the data is stored in the proprietary save file format.

A1. Numerical treatment of approach I
In section b, we noted that solutions to approach I have been obtained with the P C (Pencil Code Collaboration et al. 2021). This might seem somewhat surprising, given that this code is primarily designed for solving partial differential equations. It should be realized, however, that this code also provides a flexible framework for using the message passing interface, data analysis such as the computation of probability density distributions, and input/output.
To compute the probability distribution of with approach I, we need to sum up sequences of random numbers for many independent realizations of drawn from an exponential distribution. We use the special/lucky_droplet module provided with the code. Each point in the computational domain corresponds to an independent realization, so each point is initialized with a different random seed. The domain is divided into 1024 smaller domains, allowing the computational tasks to be performed simultaneously on 1024 processors, which takes about 4 min on a Cray XC40. A2. Dependence on initial s / grid and d / s In this appendix, we first test the statistical convergence of ( ) for the initial number of superdroplets per grid cell, s ( 0 )/ grid . As discussed in section 2.b, we set s ( 0 )/ grid = 4 for 1-D simulations. Using the same numerical setup, we examine the statistical convergence of ( ) for different values of s ( 0 )/ grid . As shown in Fig. A1(a), ( ) converges even at s ( 0 )/ grid = 1. This is important because one can use as few superdroplets as possible once grid is fixed, without suffering from the statistical fluctuations.
The most practical application of the superdroplet algorithm is the case when ≥ 1. Thus, we investigate how affects fluctuations by performing the same 1-D simulation as described in section 2.b with different values of ( 0 ). Fig. A1(b) shows that ( ) is insensitive to ( 0 ), which suggests that the superdroplet algorithm can capture the effects of fluctuations regardless of the value of ( 0 ). This is different from Dziekan and Pawlowska (2017), who found that the approach can represent fluctuations only if d ( 0 )/ s ( 0 ) ≤ 9.

A3. Horizontal variations of droplet densities
In this appendix, we analyze in more detail the effect of horizontal variations of droplet densities discussed section b. This is relevant for computing the 3-D distribution function from a 1-D distribution function. The LDM applies to a given value of the number density. Other columns have somewhat different number densities and therefore also different mean cumulative collision times. The LDM with approaches I-III can be extended to include this effect by computing cases with different number densities and then combining ( ) and normalizing by the for the combined ( ). This can be formulated by introducing the column density as where 1 and 2 denote the vertical slab in which the first collision occurs, and using this Σ( , ) as a weighting factor for the 1-D distribution functions 1D ( ) to compute the 3-D distribution functions as (A2) Since the first collision matters the most, we choose 2 = max (where the lucky droplet is released) and 1 = max − 2 / 2 (where it has its first collision). Our reference model had a number density of 0 = 10 8 m −3 . We now consider compositions of models with different values, where we include the densities (i) 0.9 × 10 8 m −3 and 1.1 × 10 8 m −3 , as well as (ii) 0.8 × 10 8 m −3 and 1.2 × 10 8 m −3 , and finally also (iii) 0.7 × 10 8 m −3 and 1.3 × 10 8 m −3 . All these compositions have the same mean droplet number density but different distributions around the mean. We first average the distribution function and then normalize with respect to the mean collision time for the ensemble over all columns. The parameters of the resulting distributions are listed in Table A1 for three compositions with different density dispersions. We see that, as we move from composition (i) to compositions (ii) and