{'fraction_non_alphanumeric': 0.037339523154516274, 'fraction_numerical': 0.01114740944520862, 'mean_word_length': 4.618579938818225, 'pattern_counts': {'":': 1, '<': 8, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, '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}}
Premelting layer during ice growth: role of clusters

The premelting plays an important role in ice growth, but there is a significant gap in our knowledge between the atomistic premelting surface structure and the macroscopic growth mechanism. In this work, using large-scale molecular dynamics simulation, we reveal the existence of clusters on premelting surface, as an intermediate feature bridging the gap. We show the spontaneous formation and evolution of clusters, and they form a stable distribution determined by the growth rate. We demonstrate how this stable distribution is related to the growth mode of ice, connected by the growth of clusters. We come to a bilayer-by-bilayer growth mode under simulation-reachable high growth rates, but another mechanism, namely"cluster stacking", is speculated to exist at lower growth rates. Our work builds a connection between the microscopic structure of premelting layer and the macroscopic growth of ice, making a step forward to the full understanding of premelting and ice growth.

A direct conclusion of such an evolutionary process is that small clusters are unstable. Because all remaining clusters are large, small clusters (if they cannot grow to large ones) have to disappear. This is the result of thermal motion, which tends to break up clusters. This explains why large models are required for observing clusters in simulation (see Supplementary Information Section 2).
It is also noteworthy that, though single molecules can move in QLL, clusters tend to stay still once they are formed. This can be seen by comparing the last two images (100 ns and 250 ns) in Fig. 1(c): clusters with the same labels (numbers and geometries) stay at almost same locations. The physics behind here, as we believe, is that the clusters in layer 1 affect the arrangement of molecules in layer 0. Though liquid water and solid ice co-exist in both layer 0 and clusters, the region in layer 0 covered by clusters is almost exclusively filled with ice (see Supplementary   Information Section 3). The solid-like structure under clusters provides stable support for the clusters, preventing them from moving freely. The stability of clusters plays an important role in the description of QLL, as we will see through the paper.
It should be pointed out that, however, such stability doesn't prevent clusters from disappearing. This is shown in the 3 rd image of Fig. 1(c), where some clusters (labeled with red circles) disappear later (4 th image) in the simulation.
In a thermodynamic view, we believe this reflects that a lot of small clusters are not as stable as fewer but larger ones due to higher "edge" energy (as an analogy of surface energy, proportional to the total circumference of clusters). Therefore, we can expect clusters to further merge after 4 th image if there is sufficient time. However, as we will show in the next section, this is not true for the QLL during ice growth, and such difference is what makes QLL during growth distinct in perspective of clusters.

QLL during ice growth
Now we shall move on to the QLL during ice growth. Inspired from the last section, here we mainly focus on the clusters in QLL. To simulate ice growth, water molecules are added randomly towards the ice surface at controlled rates. In addition, some water molecules on the surface are deleted beforehand from initial structure, to guarantee no clusters form at the beginning of simulation (see Methods section). This allows us to observe the complete process of cluster formation and growth. Fig. 2(a) illustrates the horizontal distribution of molecules in layer 1 in a typical simulation of ice growth. As we expected, no clusters are found in layer 1 at the beginning (112 ns, 1 st from left). With water molecules continuously coming to surface, the layer 1 becomes denser and after some time clusters appear (225 ns, 2 nd from left). Similar to the equilibrium case, a lot of clusters exist at this time but most of them are small. Clusters further merge thereafter (375 ns, 3 rd from left), resulting in fewer and larger clusters, distributed roughly uniformly in layer 1.
Clusters keep growing after that, but the number of clusters and their location do not change anymore (450 ns, 4 th from left).  The fact that the number and location of clusters stop changing after a certain time (like the 3 rd image in Fig. 2(a)) is remarkable. The location of clusters getting fixed also occurs under equilibrium, but the number of clusters stopping changing is a unique behavior in ice growth scenario and does not happen under equilibrium. This behavior is better illustrated by the variation of cluster size distribution over time, see Fig. 2(b). In the first half of the simulation (blue to green lines), the cluster size curve goes upwards in a smooth way. After a certain time, small clusters start to shrink while large clusters keep growing, resulting in a "cliff" in the cluster size curve. As time goes further (orange to yellow lines), clusters on the "right" side of the cliff continue shrinking until disappear. On the contrary, clusters on the "left" side of the cliff (we will call them "stable clusters") continue growing, making the cliff steeper. This trend continues even after all small clusters have disappeared. Due to the continuous supplement of water molecules during ice growth, remaining stable clusters can grow simultaneously without any of them dissipating. As a result, the number of clusters does not change anymore.
Further we studied how the growth rate of ice affects the number density of stable clusters. In this paper, we use the inverse growth rate (IGR), the average time of adding a water molecule in unit area, as the measure of how fast ice grows (note that higher IGR means slower growth). We found that the inverse number density of stable clusters (= total area divided by the number of stable clusters) is positively linear dependent on IGR (see Fig. 2(c)). This relation can be qualitatively understood as faster growth increase the possibility of cluster formation and suppressing cluster shrinking. The linear dependence can be further explained by a simple model assuming free molecules in layer 1 are moving like random walks (see Supplementary Information Section 4).
It is worth noticing the significance of these stable clusters to real systems. Under realistic oversaturation, the growth rate of ice is much lower than those simulated in this work 36 . That means fewer stable clusters in a given area, or large "white space" between stable clusters. This resulting in the stable cluster distribution dominates the process of ice growth (see Supplementary Information Section 5), and determines its growth mechanism. This is the foundation of our discussion in the following sections.
Growth mechanism of ice Based on the picture of clusters, now we can discuss the growth mode of ice. Since cluster distribution is stable (yellow frame in Fig. 3), the growth mode of ice fully depends on the growth scheme of clusters. Theoretically, there are two possible schemes for cluster growth. In the first scheme (green frame), the clusters extend horizontally and keep the height of a single bilayer (1 st sketch). The clusters continuously extend until merge together (2 nd sketch) and form a full bilayer (3 rd sketch), then the next bilayer of ice grows from here. In this paper, we name this the "bilayer-by-bilayer" mode of ice growth. In the second scheme (orange frame), the cluster extends both horizontally and vertically at the same time. Since ice is divided by bilayers in the vertical direction, the vertical growth of clusters means that new clusters (which have bilayer structure) form on top of existing ones (1 st -2 nd sketch). Because clusters do not move after forming, new clusters will remain at their original positions after lower clusters merge into full bilayers (3 rd -4 th sketch). In this paper we name this as the "cluster stacking" mode of ice growth. It is worth noticing that in this mode the clusters have similar behavior as droplets observed in experiments (exist during ice growth, have similar macroscopic appearance, look stable from surface, etc.).
It might be noticed that the two growth modes of ice here are somewhat similar to the ones in classical crystal growth theory 36,41 . It is yet unknown how to establish such theories for premelting surfaces though, due to complex structure and dynamics of QLL. By utilizing the concept of clusters and large scale molecular simulations, we can directly model the growth of ice without relying on empirical assumptions or external thermodynamic data. To determine which growth mode actually happens in a given condition, we shall focus on the core difference between them: whether new clusters can form on top of existing ones. If new stable clusters can form on existing ones, the cluster stacking mode is preferred; otherwise the ice grows in the bilayer-by-bilayer mode. To this end, let's have a closer look at the growth process of clusters, see Fig. 4: Generally, new water molecules are continuously falling on the cluster during ice growth (vertical arrows in Fig. 4).
Roughly speaking, since clusters already have bilayer structure (see Supplementary Information Section 1), new molecules are unlikely to insert directly into the cluster (layer 1, colored red) but instead they accumulate above the cluster (layer 2, colored yellow). In the meantime, molecules may leave layer 2 by diffusing to the edge of the cluster (horizontal arrows) and jump into layer 1 (curved arrows). These two processes compete with each other and lead to a steady molecule number density in layer 2. In this paper the molecule number density in layer 2 is characterized by the layer concentration (LC) defined below: represents a cluster, ( ) is the area of the cluster, 2 ( ) is the number of molecules in layer 2 on the cluster (yellow molecules in Fig. 4). is a normalization factor equal to the molecule number density of perfect ice crystal (11.78 / nm 2 , see Fig. 1(b)), so perfect crystal has 100% LC and other systems are scaled accordingly.
Intuitively higher LC should favor cluster formation in layer 2, thus favoring cluster stacking. Specifically, there should exist a critical LC (CLC), where new clusters will form in layer 2 after sufficient time if LC > CLC, and impossible to form otherwise. We found that the CLC is negatively related to the region size ( ( ) here), but tends to become a constant when the region is large enough to reach thermodynamic limit (see Supplementary Information Section 6). Furthermore, since in real scenarios stable clusters are far away from each other, stable clusters can grow very large so we can safely assume they have reached thermodynamic limit. That means we only need to consider CLC under thermodynamic limit when comparing it with LC.  One quantity we are interested in is , the maximum IGR (thus minimum growth rate) for cluster stacking under a base cluster size. This is given by the boundary IGR between LC increase and decrease (remind generally lower IGR leads to higher LC). For example, from Fig. 5(a) we find ≈ 66 with base cluster size ~30 nm.
However, the accuracy of determining decreases when the base cluster is large, due to the LC boundary becoming fuzzy (e.g. in Fig. 5(d)). This is caused by the longer time required for LC to reach steady, reflecting the challenge in studying QLL dynamics. Nevertheless, we can still get a rough estimate on , which is not so quantitative but enough for our discussion below.

Discussion
In the last section, we determined under given base cluster sizes. We may restate the result as we determined , the minimum base cluster size required for cluster stacking under given IGR (since both low IGR and large base cluster size favor cluster stacking). Meanwhile, the maximum base cluster size in a real system is determined by the number density of stable clusters, which (as we know) is also determined by IGR. That points out the way to determine the feasibility of cluster stacking, namely comparing "minimum and maximum" under the same IGR: (1) , minimum base cluster size required for cluster stacking; (2) , maximum base cluster size, derived from the number density of stable clusters.
To form cluster stacks, the base cluster size must be larger than ; meanwhile, the base cluster size on ice surface is limited by size . Therefore, > means cluster stacking is impossible and the ice grows in bilayer-by-bilayer mode; and < means cluster stacking can occur and the ice grows in cluster stacking mode. With the ⁄ ratio going higher, more layers can exist in cluster stacks. Fig. 6 shows both sizes under IGR range investigated in this paper. It's clear that is larger here, indicating a bilayer-by-bilayer growth mode in this IGR range (note that the "crossover" in ultra-low IGR region is not physically meaningful; cluster stacking would not occur there, because base clusters in that region are too small to form cluster stacks). The , maximum base cluster size, is equal to the half of average distance between stable clusters (since we characterize cluster size by radius). The , minimum cluster size for stacking, is equal to the final base cluster size of corresponding simulations (34,44,53, 63 nm, respectively), so it's slightly larger than the initial sizes used in Section 3.3 (30, 40, 50, 60 nm). Such difference is minor and doesn't affect our discussion.
Unfortunately we cannot directly determine the growth mode in realistic scenarios, because the IGR is way larger than those investigated here and unreachable for today's MD simulations. However, we may touch this by how molecules in layer 2 diffuse on clusters, since it's directly related to . For full bilayers, molecules above the bilayer diffuse like random walks (see Supplementary Information Section 4), resulting in a ( + ) 0.5 behavior for (blue curve in Fig. 6). We may fit with a similar model ( + ) , resulting in c = 0.21 (orange curve in Fig. 6). The smaller exponent indicates that the diffusion of molecules on clusters is more convoluted than random walk, which is probably caused by the barrier from edges. Furthermore, we speculate such effect might persist to the thermodynamic limit, because the barrier always exists and the "diffusion" in a liquid-like system is global. Therefore, the small exponent may persist to large IGR scenarios, resulting in < and cluster stacking happening in real systems. Confirming or denying such speculation requires a systematic investigation of how water molecules transport in clusters, which will significantly enhance our understanding of QLL.
Finally, it is worth mentioning that some other factors may also affect the growth mode of ice. An obvious one is temperature, which is suggested to cause structural transition in QLL 21,29 . Such transition may lead to an abrupt change in the diffusion behavior of molecules, as well as stable cluster distribution. Another possible factor is defects, which may induce or restrain the formation of clusters, and impact the diffusion. Some experiments found that the distribution of defects and droplets are correlated on ice surface 42 , strengthening such speculations.

Conclusion
In this paper, we report a systematic molecular-level study of QLL during ice growth, under large spatial scales (70~280 nm) using mW water model. We discovered the spontaneous formation of clusters, and simulated their distribution and evolution. Specifically, a stable cluster distribution exists during ice growth, with the number density of stable clusters linearly dependent on IGR. Based on this picture we further purpose two possible growth modes of ice, namely bilayer-by-bilayer and cluster stacking. The actual growth mode is determined by calculating LC evolution on clusters and comparing " and ", and we found the bilayer-by-bilayer mode being real under simulated high growth rates. However, the growth mode under realistic low growth rates is yet unrevealed and we suggest potential researches in the future. Overall, the concept of clusters and their distribution points out a way of connecting molecular-level dynamics and macroscopic growth of ice, making a step forward to the full understanding of QLL in ice growth. QLL during ice growth. First we constructed a perfect Ih ice crystal with sufficient bilayers in z direction and (0001) surface exposed on both sides. The model is then relaxed by NVT simulation of 1 ns. After that, we extracted 4.5 bilayers (0.5 means a half-bilayer) from the top of the model as the initial structure of our simulation. We shall point out that, we believe the procedure above is not necessary, and anyone who wants to reproduce may simply use perfect ice crystals as their initial structures (like we did in the simulation under equilibrium). However, such a procedure (though overcomplicated) is unlikely to affect our conclusion, since the number of stable clusters (an integer) is relatively robust and our conclusion doesn't rely on very accurate numbers either.

Model construction
Anyway, after the construction above, we further removed some molecules from the surface to make sure no clusters may form at the beginning of simulation. Specifically, the initial structure should have LC << CLC in layer 1, and remaining molecules in layer 1 should be as scattered as possible. In simulations the bottom 1.5 bilayers are fixed and the top 3 layers are movable. ~6 nm of vacuum is preserved above the surface, but without a reflection layer.
Growth mechanism of ice. Model construction also begins from perfect Ih ice crystals same as the NVT simulation, Growth mechanism of ice. Similar to Section 2.2, but molecules are added at 3~4.5 nm above the surface (such difference is arbitrary, and we believe it's not important because mW force-field is short-ranged), and a snapshot is made every 0.05 ns.

Data analysis
Dividing layers. The lower boundary of layer 0 is the bisector of first two bilayers in the initial structure. The layer 0 is one bilayer thick (0.36 nm), and layer 1 generally represents the whole region above layer 0. In the context of requiring more than two layers (Fig. 3, Fig. 4), each layer is one bilayer thick, and the region of each layer is determined successively from layer 0. The region corresponding to each layer is fixed in the whole simulation.
Cluster analysis. Cluster analysis is performed with distance-based neighbor criterion: two molecules with distance < 0.45 nm belong to the same cluster. This is done with the "cluster analysis" modifier in Ovito. Only molecules in layer 1 are considered.
Cluster tracing (Fig. 1(d)). To get Fig. 1(d) we need to know where a cluster in a certain snapshot comes from (e.g. "cluster B in snapshot 498 comes from cluster A in snapshot 497"). This is done in the following procedure: (1) Calculate the center of each cluster in snapshot [n] (we will call them A1, A2,…) and snapshot [n+1] (we will call them B1, B2,…). Only clusters with >50 molecules are considered.
(3) Select the smallest element in the distance matrix. If it's larger than 10 nm, go to (6); else go to (4).
(4) Find the corresponding Ax and By for this element, and we will claim "cluster By in snapshot [n+1] comes from cluster Ax in snapshot [n]".
(5) Set this element to +infinity, and go back to (3).
(6) If a cluster in snapshot [n+1] doesn't come from any cluster in [n], we will claim it's newly evolved at this snapshot.

(7) If a cluster in snapshot [n] doesn't have any corresponding cluster in [n+1], we will claim it disappears in [n+1].
The choice of cutoff ("10 nm" in (3) and "50 molecules" in (1)) is theoretically arbitrary but should be reasonable to get meaningful results. Information Section 3). The method described in this article 50 (called "identify diamond structure" in Ovito) is used to determine if a molecule belongs to water or ice (hexagonal / cubic) structure. This method is originally designed for Si and Ge systems, but we found it also works well for mW water, including identifying structures in QLL. Note that the 1 st and 2 nd nearest neighbors of ice molecules are also categorized as ice molecules in this work (themselves are "on the lattice site" but their neighbors may not).

Distribution of water and ice (Supplementary
Stable cluster distribution. The stable cluster distribution is considered established when all of the following criteria are met: (1) every stable cluster has >400 molecules; (2) every other cluster has <100 molecules; (3) cluster distribution doesn't change after this time (except merging due to cluster growth). Due to thermal fluctuations, the time when cluster distribution reaches stable (used in Fig. S4) has an uncertainty of ~10 ns.

LC on clusters.
LC is calculated in the layer above the base cluster (layer 2 in Fig. 4). The area occupied by the base cluster is calculated by making a boundary polygon for molecules in the base cluster. To calculate the boundary polygon, we project the molecules to XY plane and construct alpha shapes from them. All distinct connected alpha shapes are sorted by their alpha value, and the shape with alpha value at the lower quartile is used to extract the boundary polygon. The whole process is done with the function "boundary" in Matlab.

Vertical distribution of molecules
In the "QLL under equilibrium" section of main text we studied the horizontal distribution of molecules in layer 1, and showed the existence of clusters. To further clarify the structure of these clusters, we calculated the vertical distribution of molecules (also in this simulation), at two typical time points: (1) 1 ns, where no clusters have formed; (2) 250 ns, where several large clusters exist. The results are shown in Fig. S1. Fig. S1(a) shows the vertical molecule distribution profile at two time points, along with the profile of initial structure (perfect crystal) as a comparison. It can be seen that the molecule distribution has almost reached equilibrium after 1 ns, and the shape of the profile agrees with previous studies 1,2 . However, if we look closely into the region of layer 1 (the inset of Fig. 1(a)), we can notice a significant difference between the profile of two time points (1 ns vs. 250 ns). Namely, the molecule count decays monotonically toward the surface at t = 1 ns (orange line), but a peak at ~1.2 nm exists at t = 250 ns (yellow line). Further analysis shows that this is related to the existence of clusters. Fig. S1(b) explicitly illustrated the vertical profile of layer 1, and decomposed the contribution from clusters and "free" molecules (molecules don't belong to any cluster) for t = 250 ns (all molecules are free at t = 1 ns). The profile of free molecules is always monotonically decreasing, but a peak exists for the profile of clusters (the "t = 250 ns" graph of Fig. S1(b), orange line). Furthermore, we noticed that the profile of clusters and the equilibrated first bilayer (layer 0) are very similar if we plot them together (bottom graph). This similarity indicates that the clusters actually have the structure of bilayers. Meanwhile, this also means that the clusters are actually "located" in layer 1, and relatively isolated from layer 0 (this can be confirmed by the vertical profile of clusters, which is ~0 at the boundary of two layers). Therefore, both the concept of clusters, and the division between two layers, are well-behaved. In the "QLL under equilibrium" section of main text we show that small clusters are unstable. In small simulation models, there are not enough molecules in layer 1, prohibiting large clusters to form. The remaining small clusters tend to dissipate quickly, thus are unable to be observed consistently. This is probably the reason why clusters are not reported before.
This argument can be confirmed directly by an MD simulation with a much smaller model (~20 x 20 nm horizontal, ~28000 molecules). All other simulation settings are kept same as the "QLL under equilibrium" section of main text. Results are shown in Fig. S2. Fig. S2(a) shows the horizontal distribution of molecules in layer 1 at three different times. It is possible to form clusters in such small models, as shown in the middle image (t = 150 ns). However, the cluster breaks up later (see the right image, t = 225 ns), meaning that such clusters are unstable and cannot be consistently observed. Fig. S2(b) shows the change of largest cluster size during simulation. Some indication of forming clusters can be seen (peaks in the curve), but they are too small (<200 molecules) thus cannot exist for a long time.
It is worth mentioning that, it is still possible to observe stable clusters in small models. As demonstrated in the main text later (and Supplementary Information Section 6), layer concentration (LC) affects cluster formation, and increasing LC will help the formation of clusters in such small models. This can be achieved by adding extra molecules on the surface of initial structure, so more molecules will stay in layer 1 at equilibrium. However, the number of molecules added must be carefully controlled, otherwise clusters may contact each other (or with the periodic mirror of itself).   S3 demonstrates the distribution of water (more accurately, molecules with liquid-like structure) and ice (molecules with solid-like structure) in QLL. Fig. S3(a) shows the water-ice distribution of several clusters. They are composed of both water and ice, with stacking faults occurring occasionally. Fig. S3(b) shows the same thing but for the whole layer 0, which shows similar behavior as that of clusters. Such a mixed distribution agrees with results in previous study 3 . However, this is not true for the region in layer 0 covered by clusters. The shadowed area in Fig. 3(b) shows the location of these clusters, and most of these area is composed of ice (colored yellow/blue) in layer 0. That is to say, clusters are mostly "sat" on solid ice rather than liquid water.

Figure S3. Distribution of water and ice in QLL (clusters & layer 0). Hexagonal ice, cubic ice, and water molecules are colored yellow, blue, and grey, respectively (see Methods section of main text for how they're categorized). The structure is taken from the end of simulation in the "QLL under equilibrium" section of main text (t = 250 ns). (a) four typical clusters. (b) layer 0 as a whole; shadows indicate region covered by clusters (note that clusters themselves are not shown in (b), since they're not in layer 0).
The linear dependence between number density of stable clusters and inverse growth rate (IGR) can be explained by a simple model. First note that molecules in layer 1 can be divided into two types: (1) molecules in clusters, and (2) free molecules (molecules not in any cluster). We will make the following assumptions about these two types of molecules:  Free molecules move as random walking in layer 1 until they touch any cluster;  After touching a certain cluster, free molecules will become molecules in that cluster;  Molecules in clusters always stay in the same cluster.
We will further make the following simplification:  Area occupied by clusters in layer 1 is small compared to the system size (this is usually true, see Supplementary Information Section 5), so (almost) all molecules newly come to the surface are free molecules;


The growth of clusters is slow (compared to the interval of molecules coming to surface), so the area of "free region" (region not covered by clusters) can be considered as a constant.
Since the amount of clusters does not change, the LC in free region should roughly be a constant (if LC increases to too high, new clusters will form in free region; if LC decreases to too low, existing clusters will dissipate). Since the area of free region is considered constant, the number of free molecules has to be constant. That is to say, the following two rates are equal: (1) Rate of new free molecules added; (2) Rate of free molecules touching clusters.
(1) is (by simplification) equal to the rate of new molecules coming to the surface, or growth rate. So the inverse of (1) is IGR. The inverse of (2) is the average time required for a free molecule to touch a cluster. Since free molecules are assumed to do random walks, this average time is proportional to d 2 , where d is the average distance between a random point in free region and the nearest cluster. Since clusters are (roughly) uniformly located, d is further proportional to the average distance between clusters. Then 1/d 2 is proportional to the number of clusters in unit area, so d 2 itself is proportional to inverse number density of clusters. In other words, the inverse of (2) is proportional to the inverse number density of clusters. Since (1) is equal to (2), the inverse number density of clusters is proportional to IGR. The intercept in Fig. 3(c) of main text represents the finite size of clusters, so the number density of clusters cannot be infinite. However, such effect is negligible in realistic scenarios (IGR is very large).
It's necessary to point out that the real world is more complicated. The "mean free path" of free molecules in layer 1 is actually very short, due to frequent exchanges of molecules between layer 0 and layer 1. So the free molecules are not really doing random walks. More likely, the "random walk" of free molecules in layer 1 represents the diffusion of molecules in layer 0, which arguably shows t~d 2 behavior if traditional diffusion theory is applicable.
In the "QLL during ice growth" section of main text, we have known about the stable cluster distribution during ice growth. In this section, we will discuss how long does the "stable distribution" hold in the whole process of growth.
Here we assume that the ice grows in bilayer-by-bilayer mode (see Fig. 3 of main text; for cluster stacking mode the stable distribution arguably holds at all times). In this mode, the stable distribution establishes when the number of clusters stops changing, and lasts until clusters start to contact each other. Fig. S4(a) shows the horizontal map of molecules in layer 1, right after the establishment of stable distribution. It can be seen that, with IGR increases, the number of molecules in layer 1 at this moment decreases. Fig. S4(b) demonstrates this quantitively by showing LC corresponding to (a). Therefore, under high IGR in realistic scenarios, the stable distribution starts holding when LC is very low (<10%, at most). On the other hand, since clusters are (roughly) uniformly located, the LC when clusters start contacting is mostly determined by how they arrange in the 2D plane.
In a typical example (simple square arrangement), the LC at this time is π/4 ≈ 78.5%. Therefore, during the period of growing a full bilayer, the stable distribution holds for ~70% of total time. So we can say the stable cluster distribution dominates the whole process of ice growth. Figure S4. In the "growth mechanism of ice" section of main text we claim that high LC prefers cluster formation, and the critical LC (CLC) is a constant unrelated to region size in the thermodynamic limit. However, we don't provide direct evidence there. In this section, we will study the relation between CLC and region size directly by MD simulations.
To this end, we shall first slightly extend the definition of LC. In the main text the region LC defined on is the base cluster. However, it's hard to accurately control the LC on a cluster, and the cluster size continuously changes during ice growth, making quantitative study difficult. Instead, it's more feasible to study the relation between LC and region size in a full-bilayer model under NVT simulations. In such context, we can define LC as the following: where represents a full-bilayer model, 1 ( ) is the number of molecules in layer 1 (red molecules in Fig. 1(a) Fig. S5(a)), the maximum cluster size increases consistently over time, implying that stable clusters can form under this LC (and model size). In the simulation with lower LC (3.52%, dashed lines in Fig.   S5(a)), however, the maximum cluster size goes up only for a while and then goes down again, implying that current LC (and model size) cannot support forming stable clusters. Such observation indicates that high LC favors cluster formation. (1) average LC, calculated over last 80% of the simulation; (2) final cluster size, averaged over the maximum cluster size of the last 10% simulation. This is the measure of whether stable clusters may form.
Due to the stochastic nature of such simulations (technically "forming clusters" are somewhat rare events), some abnormal behaviors are observed in Fig. S5(b) (e.g. simulation with same model size and higher LC shows lower final cluster size). However, we can still tell that high LC and large system sizes generally favor the formation of clusters.
Now we shall turn to the relation between CLC and region size. In the main text we claim that theoretically the CLC should be a constant when region size reaches thermodynamic limit. More specifically, when the region is sufficiently large (so numerous molecules exist in it), forming a cluster in the layer would only have a marginal effect on the LC of other locations in the layer. That means forming clusters effectively becomes a local behavior (have no effect elsewhere), thus its feasibility only depends on local LC and is irrelevant to the region size. An intuitive analogy for this is the crystallization from solutions, whose feasibility only depends on the concentration of solution and is irrelevant to its volume (as long as it's large enough so crystallization itself has negligible effect on the concentration).
This conclusion is also supported by the MD simulation above. The CLC under each model size can be estimated from results in Fig. S5(b), if we ignore abnormal data points (see footnotes of Fig. S5). Fig. S5(b) shows the estimated CLC for all four model sizes, using 300 molecules as the criterion of forming stable clusters. We found that the CLC is negatively correlated with system size in 70~210 nm, showing that large models can form clusters under lower LC. However, the difference of CLC between 210 nm and 280 nm is small (and the order is "reversed": the larger 280 nm model has a higher estimated CLC). The reversed order is likely unphysical, but it also implies that the effect of system size on CLC would be small in large systems (so stochastic factors may dominate), thus supporting our conclusion. This also suggests that the thermodynamic limit is achieved when region size reaches ~210~280 nm, and suggesting a CLC of ~3.3%. In Section 2.3 of main text we used a slightly lower value (3.1~3.3%) for CLC under thermodynamic limit, since such simulations may fail to catch some cluster formation events thus tend to overestimate CLC.