{'fraction_non_alphanumeric': 0.02955167333193012, 'fraction_numerical': 0.00862976215533572, 'mean_word_length': 4.588332157139496, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Directional freezing of binary colloidal suspensions: A model for size fractionation of graphene oxide

The performance of graphene oxide(GO)-based materials strongly depends on the lateral size and size distribution of GO nanosheets. Various methods are employed to prepare GO nanosheets with a narrow size distribution. One of the promising method was proposed recently by directionally freezing of a GO aqueous dispersion at a controlled growth rate of the freezing front. We develop a theoretic model of binary colloids suspension, incorporating both the moving freezing boundary and the preferential adsorption of colloidal particles to the ice phase. We numerically solve the coupled diffusion equations and present state diagrams of size fractionation. Selective trapping of colloids according to their size can be achieved by a suitable choice of the experimental parameters, such as the adsorption rates and the freezing speed.


I. INTRODUCTION
The size of graphene oxide (GO) nanosheets plays an extraordinarily crucial role in controlling the microstructures and properties of GO-based materials. This size tunability of the material properties offers great potential for applications in material science and nanotechnology. The distinct sizes of GO can demonstrate their individual merits. Large GO sheets are desirable for fabricating conductive thin films for optoelectronic devices. 1,2 On the contrary, small GO sheets are favorable for sensing and biological applications. 3,4 Furthermore, the sizes of GO sheets play an important role in controlling the electrical and thermal conductivities of GO-based materials. 5,6 In practice, however, solutions of GO nanosheets usually have a broad size distribution due to particular preparation processes. Therefore, a tremendous variety of chemical and physical size fractionation methods have been proposed and developed to achieve the GO sheets with the relatively narrow distribution in size. Shi et al. designed a size separation method based on the pH-dependent amphiphilicity of GO sheets via adjusting the pH value of dispersion. 7 Taking into account of technical difficulties of precisely adjusting the pH values of GO solutions, they further developed a method based on the exclusive diffusion of GO sheets from its aqueous dispersion through precisely size-defined pores of track-etched membranes. 8 By virtue of the peculiar dispersibility and stability of GO sheets in polar solvent, a polar solvent-selective sedimentation method was developed by Zhang et al. to obtain series of products with size homogeneity. ? Sun et al. conducted a density gradient ultracentrifugation rate separation method to sort GO sheets by taking advantage of their differences in sedimentation rate. 9 Later, Lee et al. developed a facile size selection method based on the spontaneous accumulation of large-size GO flakes within the liquid crystalline phase. 10 Nevertheless, these methods mentioned above have their individual drawbacks in practice, such as requirement of special equipment, sophisticated operating procedures to remove the chemical additives, time or energy consuming cost.
Recently, we developed a simple and efficient method to achieve the size fractionation of GO nanosheets by directionally freezing a GO aqueous dispersion at a controlled growth rate of the freezing front. 11 The stratification of GO sheets in the solution leads to a preferential adsorption of small GO sheets in the ice phase, and the final results show a narrow size distribution of GO sheets with the size ranging from several nanometers to tens of mi-crometers. This typical non-equilibrium phenomenon is simultaneously controlled by several factors, such as the size distribution of GO sheets, the freezing rate of ice front, and the preferential adsorption of GO sheets to the ice phase. However, the underlying mechanism for the size fractionation of GO sheets is still not clear.
Here, we shall provide a detailed numerical study of the directional freezing of binary colloidal suspensions. We make a crude simplification and model the GO sheet as a hard sphere of similar size. The spherical approximation should be valid for nanometer-size GO sheets due to the rapid rotational motion in comparison to their translational Brownian motion, but is questionable for GO sheets of micrometer size. Here we proceed with this approximation and show later that a qualitative agreement with the experiments can be achieved. In the present work, we primarily focus on the freezing front moving with a low velocity when the instability of the freezing front does not happen. 12 In this case, the ice front can be regarded as a planar-like interface, and the evolution of the system be reduced to be a one-dimensional problem. Furthermore, the slow freezing rate allows us to consider the trapping of colloids in the ice phase as a quasi-equilibrium process, i.e., the preferential trapping of colloids is governed by the equilibrium of chemical potentials between ice phase and liquid phase. More complete treatment involves the intricate competition from many types of forces. [13][14][15][16] In addition, the dilute limit we concern here ensures the validity of the two-body interaction approximation describing the cross-section interaction which has been demonstrated as the crucial factor in controlling the stratification of binary colloids in drying films. 17,18 The remainder of this article is organized as follows: In Section II, we derive the dynamical equations that describe the evolution of binary colloids with a moving freezing front.
We then numerically solve these equations, and plot state diagrams of size fractionation by varying the system parameters, such as the size ratio, the freezing rate, and the difference in the chemical potentials of colloids in ice/solution phases (Section III). We then undertake the comparison of our predictions with experimental observations and demonstrate the mechanism of fractionation. Finally, we conclude in Section IV with a brief summary.

II. THEORETICAL MODEL
In this section, we set up a diffusion model that describes the distribution of colloidal particles in the ice and in the solution, as the freezing front moves from top to bottom.
We focus on the case of slow freezing rates, where the advancing ice front keeps a planar profile 12,19 and some colloids will be trapped in the ice due to adsorption. However, at faster freezing rates, a non-planar front arises 20,21 due to interfacial instability during the solidification of colloidal suspensions, particularly for the constitutional supercooling where the temperature of freezing front is below its equilibrium freezing temperature. The system out of equilibrium with non-planar interface is extremely complicated, related to the intricate interaction between the freezing front and colloids, and is beyond the scope of our current study. Here, we only consider the fractionation of the planar freezing front for colloidal mixtures. In addition, the lateral flow can also be ignored because of the relatively low concentration we consider here. Thus, our system can be regarded to freeze one dimensionally.
The effect of the gravity is also neglected in our model.
We consider a colloidal solution that is composed of two types of particles with different sizes, as schematically shown in Fig. 1. The solution is frozen from the top surface, i.e., z = h 0 , and the freezing front moves downward along the −z direction with a constant speed V f r . The evolution of the solution height satisfies h(t) = h 0 − V f r t. The colloids are considered as hard spheres with the radius r S (small) and r L (large), and their volumes are v i = 4πr 3 i /3(i = S, L). We define the size ratio of colloids by α = r L /r S . Because of the interaction between the freezing front and the colloids, some particles will be trapped inside the ice as the freezing front moves. The time-dependent volume fraction can be denoted by φ(t, z ∈ [0, h]) in the solution and φ Ci (t, z ∈ [h, h 0 ]) in the ice. The number density can be related to volume fraction by φ i = n i v i . Initially, the volume fractions for both species in the solution are assumed to be homogeneous, i.e. φ i (t = 0, z) = φ 0i .
We shall derive the evolution equations for the colloidal volume fractions using the Onsager variational principle. This principle essentially states that the time evolution of the system is determined by the balance of two forces: the thermodynamic force that drives the system to the state of equilibrium, and the viscous frictional force that resists this change. 22 The evolution equations can be determined by the minimization of Rayleighian R =Ȧ + Φ, whereȦ is the rate of change of free energy A, and Φ is the dissipation function accounting for the frictional force exerted by the fluid when particles move.

The free energy can be expressed by
where the free energy density for the hard sphere mixture reads In the free energy expression (2), we have used the local density approximation (so no gradient terms come into the free energy) and kept up to the second-order virial expansion. 23 These approximations should be valid for dilute solutions. In Eq. (2), k B and T denote the Boltzmann constant and absolute temperature, respectively. In addition, the second-order virial coefficient is given by The number of particles is conserved, leading to the conservation equatioṅ where υ i is the averaged particle velocity. The time derivative of free energy can be calculated using the integral by partṡ where the chemical potential µ i is given by The dissipation function is a quadratic function of υ i , where ξ ij is the friction matrix. In the present work, we ignore the hydrodynamic interaction and assume the friction matrix is diagonal: where ζ i denotes the friction constant per particle.
Minimization of the Rayleighian R =Ȧ + Φ with respect to υ i leads to the expression for the average velocities υ i Here, the Einstein relation D i = k B T /ζ i has been used.
The evolution equations for the colloidal concentrations can be obtained by the substi- Equations (9) and (10) are coupled diffusion equations which should satisfy a set of boundary conditions. At the bottom, i.e. z = 0, the boundary condition is υ i = 0. However, at the top surface of the solution i.e. z = h, the situations ranging from trapping and rejection are complicated due to the intricate interactions between dispersed colloids and the freezing front, such as thermodynamic forces and hydrodynamic forces. 13,14,16 Nevertheless, regardless of the elaborate competition among these range-dependent forces, here, we take the simplest assumption that the quasi-static equilibrium condition is satisfied as the ice front moves at a slow speed. Then, we can determine this boundary condition by the equivalence of chemical potentials at interface between the ice phase and the suspension phase for the identical The chemical potential in the ice is and its counterpart in the solution is Note here µ 0 Ci represents the change of the chemical potential for one colloids previous in the solution phase enters in the ice phase. Therefore, this quantity is related to the preference of colloids to either the ice or the solution.
Substituting Eqs. (12) and (13) into Eq. (11), we can obtain the boundary condition at Here k i is the Boltzmann factor of the quantity µ 0 Ci . In the very dilute limit if one neglect the second-order virial contribution, k i = φ Ci (h)/φ i (h) gives the ratio of the colloidal concentrations in the ice phase and in the solution phase at the freezing front. When µ 0 Ci > 0, k i < 1 and φ Ci (h) < φ i (h), the colloids prefer staying in the solution phase. When µ 0 , the colloids are adsorbed into the ice phase.
Associating with the condition that particle numbers are conserved across the interface

the boundary condition at the freezing interface reads
Here the boundary conditions at the freezing front are coupled as well.
The coupled diffusion Eqs. (9) and (10)  The dimensionless diffusion equations are This procedure leads to two Peclet numbers which quantify the competition between the freezing of the solution and the diffusion of the colloids. For Pe < 1, the system is diffusion-dominant and tends to exhibit a uniform distribution. On the contrary, for Pe > 1, the colloidal accumulation will lead to the creation of a concentration gradient.
In solving the coupled Eqs. (9) and (10), normal numerical routines are unstable when the concentrations occasionally become less than zero (also nonphysical). In order to avoid this problem, we perform the following substitution of variables In this way, the positiveness of φ S and φ L are guaranteed. Then the dimensionless diffusion equations become These coupled differential equations were solved using Matlab routine pdepe.

III. RESULTS AND DISCUSSION
It was reported that the larger α or strong concentration gradient of small colloids ∂φ S /∂z drives the large colloids toward the bottom, causing the higher concentration of small colloids than large colloids near the freezing front. 17 When the freezing front moves ahead, due to the adsorption of the interface between the solution and the ice, some colloidal particles would be trapped in the ice. The detail of the adsorption mechanism depends on many factors, such as the size of the colloids, functional groups on the surface, etc. To be general, here we consider two possibilities: (I) k S < k L and (II) k S > k L . We study different values of α and Pe S , and numerically solve the coupled diffusion equations (20) and (21). The total concentrations of the colloids in the ice are then calculated by integration.
A. Case (I): k S < k L At first, we analyze the case k S (= 0.04) < k L (= 0.08), i.e., the freezing front has a weaker adsorption effect on small colloids than large colloids. We select three representative sets of α and Pe S . The concentration profiles are shown in Fig. 2. The left and middle columns show the concentration profiles in the solution for small colloids and large colloids, respectively.
Each curve corresponds to a different time. The right column shows the concentration profile in the ice, and we stop the calculation when half of the initial solution is frozen. The initial concentrations of two type of colloids are small, φ 0S = φ 0L = 0.004.
In Fig. 2(a), we show the case of Pe S = 5 and α = 3. In this case, the Peclet numbers for both types are great than unity, indicating that the diffusion is slow in comparison to the movement of the freezing front. Therefore, both types of colloids are accumulated near the ice front. The concentration of the large colloids near the front is higher than that of small colloids because Pe L = 3Pe S > Pe S . Since the large colloids also has a higher adsorption rate, k L > k S , more large colloids are trapped inside the ice phase, shown by the in-ice concentration profiles in the right panel of Fig. 2(a).
When Pe S = 15 and α = 5 [ Fig. 2(b)], the small-on-top structure appears in the distribution of colloidal particles in the solution, similar to the case of evaporation-induced stratification. In the ice phase, at the early times, most of colloids gather at the interface, and k S < k L makes more large colloids to be stuck in the ice. At later times, the concentration gradient of small colloids becomes large, driving big colloids to the bottom eventually.
The number of large colloids decreases quickly at the interface, so the concentration of large colloids in the ice will also be reduced. This leads to a broad peak in the concentration profile in the ice for large colloids.
When the Peclet number and the size ratio increase further: Pe S = 25 and α = 8 [ Fig. 2(c)], the small-on-top structure fully develops in the solution due to large values of α and Pe S . In the ice phase, the concentration peak of large colloids becomes sharp and shifts In addition, the average concentration of two types of colloids in the ice can be obtained by integrating volume fraction profiles (the right column in Fig. 2) at time τ = 0.5 This value is a good indicator for the efficiency of fractionation. Comparing the average concentration of small colloids with that of large colloids, a state diagram (Fig. 3) in the Pe S -α plane is acquired. The circles are for systems that the average concentration of large colloids is higher than small colloids, while squares show the case that large colloids is less than small colloids in the ice. The color scale is the numerical value of log(φ SC /φ LC ). The separation of the large/small colloids is more efficient in two parts of the state diagram: In the range when Pe S is small while the α > 5, most particles in the ice are large colloids [ Fig. 2(a)]. When the parameters Pe S and α are both large, small colloids predominate in the ice phase. The later phenomenon corresponds to the Fig. 2(c).
We consider the other condition: k S (= 0.1) > k L (= 0.01), i.e., the adsorption of small colloids at the freezing front is stronger than that of large colloids. We also select three typical examples via substitution of representative sets (α, Pe S ). Figure 4 shows concentration profiles of small colloids and large colloids at various times in the solution and in the ice.
The initial concentrations of two type of colloids are φ 0S = φ 0L = 0.003.
When Pe S = 5 and α = 3 [ Fig. 4(a)], the concentration gradients for both colloids are small due to the relatively small values of the parameters. There is a small rise of the concentration near the freezing front on account of Pe S ¿1, meaning that the freezing front moves faster than the diffusion. Since the adsorption rate for small colloids is larger than that of large colloids, the in-ice concentrations of small colloids is also larger than that of large colloids.
A different scenario appears when Pe S = 15 and α = 6 [ Fig. 4(b)]. At the beginning, the large colloids gather at the interface because Pe L > Pe S . In later times, the concentration gradient of the small colloids is large enough, pushing large colloids to the bottom. Similar to the intermediate state in Sect. III A, a broad peak in the concentration distribution of large colloids in the ice appears. In the ice phase, small colloids are trapped more than large colloids initially due to k S > k L . Later, due to the accumulation of larger colloids in the solution near the interface, the concentration of large colloids in the ice gradually increases.
Finally, the concentration gradient of small colloids becomes stronger, pushing the large colloids away, and the large colloids trapped in the ice are reduced.
When we increase the parameters even further, Pe S = 25 and α = 10 [ Fig. 4(c)], the small-on-top structure is formed at very early times. Due to the depletion of the large colloids near the ice front, the peak of large colloids in the ice shifts to right and decreases sharply, and most of the colloids trapped are small ones.
Similarly, we integrate volume fraction profiles in the ice at time τ = 0.5 to obtain the where the small colloids are trapped more in the ice than the large colloids (indicated by the red regions). The cases in Fig. 5 when Pe S and α are both large is similar to that in Fig. 3; the small colloids are trapped predominately due to the small-on-top structure developed in the solution. When both Pe S and α are small, there are also more small colloids in the ice than the large colloids. However, the mechanism for this case is mainly that the adsorption rate of small colloids is large. Our theoretical model demonstrates a delicate interplay between the adsorption rates and freezing speed. The size fractionation could be realized by careful control of the adsorption rates (by surface modification of GO sheets) and the directional freezing speed. When the case that more large colloids in the ice is desired, one should choose k S < k L while keep the freezing slow. In the opposite case of more small colloids in the ice, one can resort to a large freezing speed, but not too fast to be in the regime of constitutional supercooling.
There are a few limitations of our theoretical model and here we discuss possible extensions and improvements. We have employed a simple diffusion model and kept up to the second-order virial contribution to the free energy. These assumptions become questionable when the solution becomes concentrated. In general, more sophisticate free energy models can be used at high colloidal concentrations [24][25][26] . For realistic GO nanosheets, the planar geometry is important and one has to consider the orientation degrees of freedom. There are also free energy models for disks and plates [27][28][29] . In practice, size fractionation needs to be performed for commercial GO solutions, which tend to be in the high concentration regime. In experiment, however, the range of the freezing rate for effective fractionation narrows at high concentrations 11 , possibly due to the inhibition of the diffusion. Therefore, we have limited our studies to dilute solutions. More complicated separation methods may be preferred for concentrated solutions 7? -10 .
Another important effect is missing from our model is the hydrodynamic interaction. We are using an implicit solvent model and the solvent is treated as a viscous background. One can incorporate the hydrodynamic interactions by including the concentration-dependence in the diffusion constant. There are many recent studies on the effect of hydrodynamics on the non-equilibrium colloidal dynamics, mostly by evaporation [30][31][32][33][34][35] . One study on polymer solutions by Statt et al. 32 showed that stratification happens in implicit solvent simulations, but not in explicit solvent simulations when hydrodynamic effects are included. One recent study on colloidal solutions by Tang et al. 35 , however, observed the stratification for both explicit and implicit solvent simulations. Clearly there are more to be done to clarify the effect of hydrodynamic interactions.
In our model, the adsorption of the colloids to the ice phase is assumed to be in a quasiequilibrium. This simplification thus allows us to use one simple parameter, k i in Eq. (14), to characterize the adsorption rate. The real situation is more complicated. There are many forces involved when a solidification front approaches a colloidal particle: (a) The disjoining force, originated from the thin layer of fluid in between the colloid's surface and the ice front.
This force is the result of the temperature gradient and is of the van der Waals type 13,14,16 . For large particles and fast freezing speed, the dynamics is governed by the balance of the disjoining force and the hydrodynamic force. This situation has been studied extensively 15,21,36,37 , and the focus is on using the fast freezing to create complex composite materials 38,39 . Here we purposely work in the regime of low freezing speed to avoid the instability of the freezing front. In this case, especially for small colloids, the dynamics is determined by the balance of the specific interaction, entropic force, and the viscous friction.
The specific interaction can be the results of hydrophobic interaction, the hydrogen-bonding interaction, etc., mostly of the short-range type. A recent experimental work on the direc-tional freezing of emulsions demonstrated that long-range interactions are also important 40 .
In our model, the effect of these interactions are grouped into one phenomenological parameter k i . More works have to be done to determine the connection between the k i parameter and the detailed interaction between the colloids/nanosheets and the freezing front.