Effect of Hydrophobic Silicon Substrate on the Ice Formation at Different Temperatures

In this work, we investigate the effect of silicon substrate and temperature on the ice formation through molecular dynamics simulations. It was found that silicon substrate accelerates the ice formation rate. Depending on the temperatures of the systems, the formation rate can be accelerate or decelerated. The obtained results showed that, as the temperature increased from the lowest temperature 100K to the highest temperature 220K, the number of water molecules transformed into ice decreased. We found that the hydrophobic property of silicon is not a barrier for approaching ice formation. The high correlation between Si and the hydrogen atoms of water molecules stimulated the ice formation on the surface of silicon and allowed the formation of amorphous ice directly on the surface of silicon. The temperature effect on ice formation was in distinguishing two different behaviors one from 100K to 180K and the second above 180K. These behaviors were related to the temperature of which amorphous ice is stable. The results of the average number of hydrogen bonds and their lifetime are in parallel with the results of the number density proles that can best explain how silicon affects the ice formation. The coordination number of water molecules increased with the decrease of the temperature for silicon substrate which means that bigger ice clusters were found. The ice molecules, which were formed near the silicon substrate, are more recognized than that of the bulk. To better understand these effects, bulk systems were used as comparative systems.


Introduction
The study of water behavior and ice nucleation has been important for many investigators. Many experimental and theoretical studies were carried on ice nucleation in the presence of different substrates and systems under various conditions. Understanding how ice is formed is worthy since nucleation is becoming an important approach in many elds such as atmospheric studies [1,2], water transport in plants [3], microbiology [4] ,and food preservation [5].
Some of the factors that can provoke or inhibit ice formation are: pressure, the addition of a seed, and electric elds. Investigating the effect of applying pressure on the ice formation phenomenon showed that ice melts under pressure while it reforms when the pressure is lifted [6]. In recent years another method was detected which is the addition of a "seed"; such as ice; to the simulation from which nucleation will start [7,8]. Several conditions were studied to explore their link with ice nucleation whether in inhibiting or provoking. The application of an electric eld in simulations favored the formation of cubic ice structure because electric eld polarizes ice and hence cubic ice can exist in the poles, whereas hexagonal ice clusters do not form under polarizable conditions [9][10][11].
Other factors related to the substrate type and its properties affect the ice nucleation. If we want to explain the hydrophobicity of the substrate, it can be concluded as follow: the ice-nucleating e ciency will be lost if adsorption is too strong due to the high coverage of water molecules. The energy of adsorption will increase with the increase of the hydrophilicity of the substrate [12]. Lupi et al. found that graphite surface can promote ice nucleation, but the structure of the surface can affect such nucleation as well. They also concluded that rough surfaces inhibit nucleation while the at ones enhanced it [13]. Moreover, the presence of graphite favors the layering of water molecules [14]. Layering is established by the ordering of the rst layer and then the second layer will be stacked directly above the rst one [13,15,16]. Another study [17] that depended on silica and graphite substrates was studied. The results manifested that the density of water increased in the presence of these substrates where the solid substrate perturbs the interfacial water layer for only a few molecular diameters away from the surface.
Kaolinite surface is a pretty interesting inorganic surface when studied. The hydroxyl groups play the role of the donors and the acceptors of hydrogen bonds [18]. The hydroxyl group reorients in the presence of aluminum Al to form a hexagonal ice structure. While, from the site of rigid Si of the kaolinite, layers of hexagonal and cubic ice patterns are seen with no need for hydroxyl groups reorientation [19]. Alterminated and OH-terminated alumina surfaces showed that the nature of the surface termination determines the organization of interfacial water molecules yielding well-de ned patterns [20].
Temperature takes a vital role in controlling ice formation or melting. In the case of selecting temperature few degrees below the melting temperature, depending on the hydrophobicity of the substrate, it can lead to the formation of ice [21]. The formation of ice passes through two stages: First, the pre-critical ice nucleus will be formed during a period named the induction period, then, it will keep forming and remelting until it will be stable and becomes a critical nucleus and from this seed, the ice will start forming [22]. When choosing temperature close to the melting point the growth rate of ice becomes slower for the reason that there is a rapid reduction in the translation motion of water then ice-like crystal will be formed and then the ice-like layer will melt. Recent studies attempted to enhance ice nucleation by using Nanogrooves. The work done by Zhang using a silicon substrate with nano-groove indicated that the roughness of the hydrophobic substrate does not affect nucleation. The nucleation rate was found to be sensitive to the size of the groove where ice can be formed inside the groove more than on a at surface [23].
The effect of surface defects was also investigated to interpret its consequences on ice nucleation. Silver iodide (AgI) had captured the attention in these studies due to its ability in ice formation even at 270K [24]. One of the important results that were seen is the presence of both hexagonal and cubic ice structures on AgI [25]. Defects over silver iodide surface depicted that defects did not affect the ice nucleation. However, the surface area around the defects reduced the nucleation rate. This reduction is due to the perturbation of the hydration layer of the surface which will lead to retarding the critical nucleus formation [26].
To describe the ice nucleation there should be a bulk reference to compare the results with and especially when using water on a substrate to know what effect does the substrate has [27]. The work done by Qiu aimed to demonstrate whether ice can nucleate on the graphite surface or not [28]. The results described that the thermodynamic properties of water in the interface layer are similar to that of bulk; thus nucleation can take place. The results of the bulk water system were also used to well-understand the upshot of graphene oxide con ning walls on the phase transition temperature [21]. It was seen that the phase transition temperature was less with 34 K in the presence of con ning walls compared to that of bulk water. Another approach to generalize the behavior of hydrated water in the presence of proteins revealed that this water, depending on radial distribution function (RDF) results, had the same behavior of bulk after 8 Å. Before 8 Å the peaks were located at the same distance but those for hydrated water were sharper [29]. Hence, the aim of this work tends to analyze the effect of hydrophobic substrate on ice formation by comparing the results with bulk data at different temperatures.

Simulation details of bulk system
In order to know the effect of the hydrophobic substrate on the ice formation, the bulk system was used as a reference system. The simulations were performed using LAMMPS [30] and TIP4P/ice model [31], as a model that was designed to cope with the solid-phase transition. TIP4P/ice is considered as a modi ed potential model of TIP4P model. It is an attractive model due to its simplicity, in addition to giving accurate results without increasing the computational time. For this reason, it becomes highly used in both equilibrium and non-equilibrium studies. Predictions of this model for both densities and coexistence curves of ice have been demonstrated to be better than those for other water models. It intends to give great investigation about hexagonal ice, cubic ice, and denser ice forms [31][32][33][34]. In this system, epsilon for oxygen and hydrogen atoms were 0.21084 and 0.0 Kcal/mol respectively. The sigma of the used atoms were 3.1668 Å and 0.0 Å respectively. Partial charges (q) applied for oxygen and hydrogen were − 1.1794 and 0.5897 e, respectively. All the interatomic interactions were modeled by the Lennard-Jones potential. The Nose-hoover thermostat and barostat [35] [36] were applied to x the temperature and pressure respectively. The simulations were carried out using periodic boundary conditions in the three directions. The number of hexagonal ice molecules was 1728 molecule with 1500 water molecules set in a box structure. Simulations were performed at different temperatures with isothermal-isobaric (NPT) ensemble and the pressure equals to 1 atm. The time step was 2 fs conducted for 20 ns.

Simulation details of system with ice and water over silicon substrate
Silicon substrate consisting of 1584 atoms was applied. Using the same dimensions and number of atoms as the bulk system, a box of hexagonal ice and a box of water was set over the silicon substrate. The hexagonal ice box was added to accelerate the ice formation. The typical size of simulation box was 84 x42 x28 Å 3 .The interaction between the water molecules and silicon was described by 8-10 Lennard-Jones potential, with σ Si = 3.387 Kcal/mol and ε Si = 0.585 Å. Periodic boundary conditions were applied in all directions. The velocity Verlet [37] algorithm was used to solve the equations of motions. Shake algorithm to x the bonds were also applied. 10 Å cutoff was applied for all non-bonded interactions.
Several NPT molecular dynamics simulations were carried out at four different temperatures. CHARMM force eld was utilized [38]. Temperature and pressure were stabilized using Nose-Hoover thermostat and Nose-Hoover barostat respectively. The temperatures used in this study were 100K, 140K, 180K, and 220K. The simulations were performed for 20 ns with time step of 2 fs. Figure 1 represents the con guration of the studied systems.

Radial distribution function
The rst peak in the bulk system was found at 2.74 Å with a good agreement with the experimental data [39,40]. By comparing the results between the different silicon systems rst, we found that the peak of the oxygen-oxygen radial distribution function (RDF) is the highest for both Si-100 and bulk system at 100K, in which the intensities of the rst peak reached 7 and 6 respectively. The intensity of the peak decreases as the temperature was increased. This is due to less accumulation of water in the liquid state compared to that of solid state. In other words, the correlation of ice molecules with the neighboring is higher than the correlation of liquid. Figure 2 reveals that at all temperatures; there is a peak exactly at 2.65 Å for O-O, which indicates the closest distance of the accumulation of water molecules relative to each other. For bulk systems, the intensities of the peaks were lower than that of silicon substrate systems. Figure 3 which represents the O-H RDFs shows that the four systems with silicon substrate are the same. The g(r) O-H exhibited a peak at 1.74 Å for the bulk systems which is in correlation with the experimental data [39]. The O-H RDF cannot give comprehensive data about the structure of water whether it is in the liquid or solid state since the peaks for bulk and systems with silicon substrate at the same temperature are almost the same with a small shift only in the second peak of bulk systems. Displaying the H-H RDF in gure 4 shows that the second peak at 100 K is the highest compared to the peak at 140 K, 180 K, 220 K, and for bulk systems. This strong correlation between H-H can provoke the cluster structure. The second peak was increasing with the decrease of the temperature for the systems with silicon substrates. What can be established from g(r) Si-H is that the high correlation is occurring at 100 K while the least one is at 220 K, whereas at 140K and 180 K the correlation is approximately the same as it is expressed in gure 5. The high correlation between Si and H inhibited the water molecules to move apart from silicon substrate and thus ice formation process directly over silicon substrate.
The evolution of the curve at 220 K was remarkably different from the other curves. This remark can be used to classify two behaviors depending on the temperature scale: one from 100 K to 180 K and another above 180 K. The results reported in reference [34] interpreted that for bulk systems as the temperature increases the intensity of the rst peak decreases and this what can be noticed in the g(r) O-O for the bulk systems.
The radial distribution function for the water molecules of Si-100 was compared with g(r) O-O of liquid water to verify the state of the water in the presence of silicon at low temperature (hexagonal ice molecules were excluded from this calculation). The progress of the curves depicted in gure 6, revealed an obvious difference and the proceeding of the water molecules of Si-100 showed an amorphous ice evolution. This result is validated by the experimental study performed by Yongchang Wu et al., in which they detected the presence of amorphous ice layer on the hydrophobic surfaces [41]. The high correlation between silicon and hydrogen means that the water molecules are not getting apart from the substrate and the ice layer is formed directly on the surface.
According to the study published in 2006 [42] explaining about amorphous ice structure, the g(r) O-O of amorphous ice is well-matched with that of amorphous ice structure at low density (similar to this work). The rst peak for the simulation study was seen at 2.6 Å with intensity of 6.4 and the second peak at 4.4 Å [43]. Oetjen stated that amorphous ice is stable from -160˚ till -125˚C and it needs a long time to change from amorphous into another type of ice [44]. The behavior that was seen for Si-H at 220 K can be explained by the less stability of amorphous ice above 180 K. This nding gives more details about the two behaviors: one between 100K and 180 K with stable amorphous ice, and the second above 180 K with unstable amorphous ice structure. In other words, low temperatures do not favor the formation of hexagonal ice. Yet, at low temperatures amorphous ice is formed rst then ordering of the clusters will take place. This ordering took a long time where the stable critical crystal was completely formed after 1µs [22].

Coordination number
By de nition, the coordination number is established from the area under the rst peak in g(r) O-O. The value obtained can give details about the presence of water or ice. When only water is present, this number increases with the increase of the temperature. In the case of hexagonal ice, it should not be above four [45]. For low density amorphous ice (LDA :density below 1 g/cm 3 ) the coordination number is also about four [46]. The LDA was mentioned here because the calculated density of the system was below 1 g/cm 3 . Figure 7 merely exhibited the coordination number. As can be observed, for the bulk systems, the coordination number was lower than that when silicon substrate is present. For the systems with the silicon substrate, the coordination number uctuated around 3.8 while for bulk it was around 3.7. At the lowest temperature in the presence of Si, the coordination number was the highest. This is due to the highest correlation for the rst neighbors (see g(r) O-O). As the correlation decreases, the coordination number decreases. In other words, as the temperature decreases larger cluster of ice is formed. Dowell et al. demonstrated that at temperatures above -160˚C small ice clusters are found [47]. This result is in parallel with our result for having the highest coordination number at the lowest temperature (100K).

The number of ice molecules
In order to understand how many ice molecules were formed, the algorithm developed by Clancy and Ba'ez was used [48]. The number of ice molecules was detected using the following equation: F i is termed the tetrahedral parameter, i is the number of water molecules, n i is the number of water molecules in the rst solvent shell and ijk are the indices of the nearest neighbors. is the angle between two oxygen atoms. Since the hexagonal ice structure is fourfold, cos should be -0.33 to have the tetrahedral structure. As mentioned earlier the number of hexagonal ice molecules at the beginning of the simulation was 1728 for all systems. It can be noted that as the distance from the substrate increases the presence of ice decreases. The same phenomena were seen for bulk systems, as going through bulk layers, the number of ice formed decreases. The effect of the hydrophobic substrate can be seen where Si substrate favored ice formation near it. Lupi et al. concluded that hydrophobicity is not a good predictor of ice nucleation ability [15]. This predication can be used to well-explain the ice formation on silicon surface despite the hydrophobic property of this substrate. The bar diagram plotted in gure 8, designated the total number of ice formed and the number of water molecules transformed into ice. The total number of ice is the summation of the hexagonal ice molecules and the water molecules transformed into ice. To obtain the number of water molecules transformed into ice, equation (1) was applied for the water molecules only and the hexagonal ice molecules were excluded from this calculation. Concerning the effect of temperature, it can be stated that, as the temperature increases, the number of ice formed decreases. Figure S1 demonstrates the effect of the distance on the number of ice formed (see supporting information, gure S 2 shows the position of z-axis). It is in correlation with the result of RDF Si-H where at 100 K the system with silicon substrate clari ed the high formation of ice at z < 5 Å. Moreover, as the distance from the substrate increases the number of ice molecules formed decreases dramatically.
3.4. Survival probability and nucleation rate 3.4.1. Survival probability Ice formation passes through two different phases: induction and growth [17]. The induction phase, in which the critical nucleus is formed. At this phase, the ice melts and reforms until it becomes stable. In other words, the water molecules start to uctuate between ice and liquid phase. This uctuation creates the subcritical nucleus. When the uctuation stops, the subcritical nucleus develops to become a critical nucleus, which is a stable seed. From this seed, ice formation launches. After the formation of the critical nucleus, which is the important seed for starting the nucleation, growth period starts. The critical nucleus develops into crystallite and cause the water in the liquid form to change into ice. According to Molinero et al., the induction phase takes 900 ns where the subcritical nucleus was formed out of 20 water molecules. After the induction phase, the growth phase started at around 1000 ns until it ceases [22]. Figure S2 reveals the beginning of the growth phase for Si-100 (see supporting information). The ice nucleus started to form at the edge of the hexagonal ice cube near the water molecules. At this level, the survival probability is calculated to know the change in the identity whether liquid or ice. The survival probability is de ned by [22]: N (τ;τ+t ) is the number of molecules that remain in phase x during the time interval between τ and τ+t. M(τ) is the number of molecules at the interval τ. I T is the number of intervals averaged over. For calculating the survival probability, only water molecules were used after determining the beginning of the growth phase. Figure 9 explains the survival probability of ice in the growth phase where τ used in this calculation was 10 ns. The probability of nding ice at the beginning of the growth phase was high for the systems in the presence of silicon. This probability reached 0.95 for Si-100, while for bulk systems this probability reached 0.4 at 220 K. The survival probability decreased slightly for Si-100, Si-140, and Si-180. Whereas, for Si-220 the descending was very sharp. Silicon substrate can help in ice formation more than bulk at low temperatures outreaching 180K. As described before the amorphous ice is stable at very low temperatures, yet at 220K this non-crystalline ice type is not stable. This instability is re ected in the curve decay for Si-220. Thus, the probability of nding ice in the growth phase for Si-220 was the lowest. From the based results for calculating the survival probability, the number of ice formed uctuated sometimes increasing and then decreasing which can assure our discussion for melting and reforming of ice during the growth phase till the stable ice nucleus is formed.

Nucleation rate
To calculate the rate of ice formation, the hexagonal ice molecules were precluded and only the 1500 water molecules were used. Figure 10 illustrates the rate of ice formation after determining the time of the launching of the growth phase. The gure re ects that as the temperature decreases, the rate of nucleation increases. Sharp decrease was seen for Si-220 where the rate descended from 8.2x 10 23 m -3 s -1 to 3.5x 10 23 m -3 s -1 . The time needed for the nucleation to start was 2 x 10 -8 fs for Si-100, Si-140, and Si-180 while that for Si-220 was doubled which demonstrates the different behavior of ice formation above 180 K. This behavior is similar to the decay realized in the survival probability, where the probability of nding ice at this temperature was the lowest among the rest of temperatures for systems with silicon substrate. It can be deduced that the survival probability is directly proportional to the rate of ice formation. As the number of ice found in the growth phase increases, the rate of ice formation becomes faster. A clear difference in the rate of nucleation between silicon systems and bulk systems was observed. For all the bulk systems, the time for nucleation was 1x 10 -7 fs. Thus, it can be demonstrated that the silicon substrate provokes the nucleation phenomena.

Number density pro le
It is clear that, the ice formation is starting from the contact region between water molecules and the hexagonal ice molecules. Thus, ice formation is in progress along x-axis. For this reason, three regions were used to well-understand this progress. The regions range between -5 Å and -50 Å, where the hexagonal ice molecules were excluded from this partition. The rst region is -20 < x < -5, the second -35 < x < -20, and the third -50 < x < -35 as gures 11, 12, and 13 depict. Figure S 2 displays the partition regions (see supporting information). At 100 K the system showed cleaved peaks for both bulk and Si-100. The average intensity of the peaks is the same for both systems. The intensity of the peaks increased obviously in the region -35 < x< -20 and the peaks became more compatible between bulk and Si-100. In the -50< x <-35 region, the intensity of the peaks decreased by distance for Si-100 but still the cleaved structure for the ice structure found. The sharp peaks with high intensities can be recognized at -48 Å for the bulk system. This structure is referred to the periodic boundary conditions. The same pattern was seen for systems at 140K and 180K. The cleaved peaks continued to appear by getting apart from the hexagonal ice molecules. The effect of temperature can be concluded by comparing the intensity of the peaks in the -35< x <-20 and -50< x <-35 region, which decreased by distance. The results at 220K expressed short intensity peaks with a noticeable difference between these results and that at different temperatures.

Average hydrogen bond and hydrogen bond life time
In agreement with H-H RDF, it is clear that at 100 K the average number of hydrogen bonds was the highest then it decreased as the temperatures increased for the systems with silicon substrates. This result is similar to the work performed by Bleszynski et al., who showed that hydrogen bonds between the TiO 2 surface and water were more than hydrogen bonds between ice and water [49]. The average hydrogen bond, expressed in gure 14, descended sharper for the bulk systems with the decrease of the temperature than the systems with silicon substrates. Hence, it can be deduced that the presence of silicon favors more hydrogen bond to be formed. The ice nucleation process can be explained by the high correlation between silicon and hydrogen. Besides, silicon substrate provoked the formation of more hydrogen bonds at low temperatures. Figure 15 represents the lifetime for all the systems. The lifetime for hydrogen bonds is also in agreement with the average hydrogen bond where it is high for systems with the highest average bonds (Si-100 and bulk 100). The stability of amorphous ice depends on the time the hydrogen bonds are broken down. The low value of average hydrogen bond for Si-220 had the least lifetime and the least number of water molecules transformed into ice. The fast breakage of hydrogen bond at 220 K caused the instability of the amorphous ice structure, and so, two behaviors were observed below and above 180K.

Mean square displacement
Mean square displacement (MSD) can give an idea about the melting behavior of the solid form. It is known that when materials are in the liquid state their mobility is higher than that of solid. The formation of clusters can inhibit the movement of the particles. Since MSD depends on the position the molecules move according to a reference point, this facilitates the analysis of the ice formation or melting. In the case where ice is melting the curve will be increasing since the molecules are free to move in the medium, while in the case of ice formation the curve will uctuate around a constant value [50]. Figure 16 indicates the mean square displacement result of the bulk systems and the system with Si substrate at different temperatures. The curve at 100K is the lowest and linear, then as the temperature increases, MSD increases. The result of bulk at 220K presented an increasing order curve which is the only the different behavior that can be seen. Thus, the effect of the substrate at low temperatures can be assessed. Ice clusters were formed near silicon substrate and the size of these clusters decreased with the increase of the temperature. As the size of the cluster increases, the movement decreases, which means low values of MSD. Besides, silicon substrate is preventing the movement of the molecules due to the interaction between silicon atoms and hydrogen atoms (explained in the Si-HRDF). Hence, at 220 K for the bulk system, the curve is sharply increasing than that of the Si-220. Moreover, the difference between Si-180 and Si-220 curves was obvious. The reason for this is that in the absence of the substrate the molecules exhibits freedom to move as the temperature increases.

Conclusion
From all the above results, it can be concluded that silicon accelerates ice formation by increasing the Si-H correlation. The high correlation between silicon and hydrogen means that the hydrophobic property of silicon does not repel the water. Water forms amorphous ice clusters at very low temperature on hydrophobic silicon substrate. The presence of silicon increased the rate of ice formation and the number of water molecules that transformed from the liquid state to the solid state as it increased the correlation between silicon and hydrogen in comparison to that of O-H. H-H RDF, coordination number, number ice, average hydrogen bond, and hydrogen bond lifetime signi ed decreasing values with the increase of temperature in the presence of silicon. It can be noticed that the silicon substrate favors hydrogen bonding as temperature decrease to provoke the cluster formation. Another effect of temperature was indicated in the simulations executed. Two different behaviors were detected one before 180 K and the other above this temperature. These two behaviors were related to the stability of amorphous ice whereby amorphous ice became less stable at temperature above 180 K. This stability depends on the correlation between silicon and hydrogen and the lifetime of the formed hydrogen bonds. Declarations Funding Figures Figure 1 Initial con guration of the system in the presence of silicon substrate The red, cyan and blue colors represent oxygen, hydrogen, and silicon atoms respectively. The initial con guration of the bulk system is the same but with the absence of the substrate.