Freezing point depression in model Lennard-Jones solutions

Crystallisation of liquid solutions is of uttermost importance in a wide variety of processes in materials, atmospheric and food science. Depending on the type and concentration of solutes the freezing point shifts, thus allowing control on the thermodynamics of complex fluids. Here we investigate the basic principles of solute-induced freezing point depression by computing the melting temperature of a Lennard-Jones fluid with low concentrations of solutes, by means of equilibrium molecular dynamics simulations. The effect of solvophilic and weakly solvophobic solutes at low concentrations is analysed, scanning systematically the size and the concentration. We identify the range of parameters that produce deviations from the linear dependence of the freezing point on the molal concentration of solutes, expected for ideal solutions. Our simulations allow us also to link the shifts in coexistence temperature to the microscopic structure of the solutions.


Introduction
Thermodynamics and kinetics of solid-liquid phase transitions rule several processes that have a massive impact on life on Earth. Freezing and melting of water and aqueous solutions deeply influence climate, atmospheric and geological processes, materials processing and food processing and conservation. At a given pressure, crystallisation and melting rates are controlled in first instance by the temperature at which the liquid and the solid phase coexist in equilibrium, i.e. the coexistence temperature (T c ). Even though T c is an intrinsic property of materials, it can be affected by external factors, for example by the high curvature of capillary walls [1] under micro or nano-confinement [2], by the confinement length [3] or by the presence of impurities, even at very low concentration [4,5].
In solutions, T c is usually reduced by an amount proportional to the solute concentration. This is commonly experienced, for example, when ice is molten at temperatures lower than 0 • C by adding salt. In the ideal case of a binary system, in which the solute is fully diluted, the depression of the coexistence temperature, T c = T solution c − T pure solvent c , depends linearly on the molality m of the solute, with a proportionality constant k. In aqueous solutions, k is a dubbed cryoscopic constant. In the ideal case, k does not depend on the type of solute but only on its concentration, so that The basic rationale behind such ideal behaviour is that the addition of solute to a pure solvent increases the entropy of the solution by an amount proportional only to the number of added solute particles, independently on the chemical nature of the solute [6]. The excess entropy shifts the phase balance between liquid and solid in favour of the liquid, thus shifting the T c . The simplicity and universality of this relation is exploited to measure the vapour pressure of solutions of electrolytes and/or their molar number [7,8]. However, Equation (1) is only justified in the case of an ideal dilution, that is at very low solute concentrations, and it was observed to break down already for 0.5 molal aqueous solutions [9]. Deviations from ideal behavior stem from specific molecular interactions, namely solute-solute, solvent-solute interactions, and solute-induced variations in solvent-solvent interactions [10]. At low concentration, solvent-solute interactions are expected to play the most important role, and empirical corrections to the ideal freezing point depression law have been devised in terms of solute-solvent interaction or activity coefficients [10][11][12]. Equation (1) can be generalised to take into account deviations from the linear behaviour due to the partial dissociation of solutes, by including the Van't Hoff dissociation factor i, so that Even if empirical models may result successful for specific cases a microscopic theory linking the thermodynamics of solid-liquid phase transitions of binary solutions to their structural properties islacking.
In this work, we verify to what extent Equation (1) holds for a simple system, and provide the missing link between the shift in the coexistence temperature and the structural changes in the solvent, induced by the presence of different types of solutes. We especially want to target the conditions often occurring in biological environments or in food processing, in which the presence of few large macromolecules produces significant shifts in the coexistence temperature of water. To this scope, we employ molecular dynamics (MD) simulations on model systems, consisting of binary Lennard-Jones (LJ) mixtures, in which one of the two components, made of larger particles, is kept at low concentration (<2.5% in number of particles) and acts as solute, and the other acts as solvent.
LJ mixtures are valuable models to understand the general features of crystallisation in binary solutions, from nucleation [13] to growth [14][15][16][17]. More complex LJ mixtures with suitably chosen interaction parameters were also used to model how additives modulate the nucleation and crystal growth of solutes in supersaturated solutions [18,19]. The phase diagram of binary LJ mixtures was extensively characterised in former works [20][21][22]. Nevertheless, these works take into account mixtures of particles with rather similar size (diameter ratio up to ∼1.2), and do not discuss the effect of large solute particles at small concentrations on the solid-liquid transition temperature of a solvent made of smaller LJ particles, which, instead, is the focus of this work.
Two types of solute-solvent interaction and a broad range of sizes and concentrations of solute particles are tested. We show that T c deviates significantly from the ideal behaviour in Equation (1) for solutes smaller than a certain size even at very low concentrations, while the effect of larger solutes on T c is correctly predicted. We then propose a link between T c and how solutes modify the short-and medium-range order of the solvent, probed by computing the static structure factor S(k) and the excess configurational entropy.

Models and numerical methods
We compute the melting temperature of a binary LJ mixture by MD simulations. Solid-liquid phase transitions in singlecomponent LJ are well characterised: the solid-liquid phase coexistence line was computed by Agrawal and Kofke [23], and several studies investigated nucleation from the melt [24][25][26][27] and crystal growth [28,29].
In our simulations of LJ binary mixtures, the interaction between two particles of type a and b at distance r is described by the truncated and shifted LJ potential [30], where v ab (r) = 4ε ab σ ab r 12 − σ ab r 6 . (4) We define type-1 particles as solvent and type-2 particles as solute (at low concentration), and σ 11 = 1 and ε 11 = 1 as length and energy units. All particles have unitary mass and all the other quantities are expressed in reduced units derived from σ 11 and ε 11 . For example, time is expressed in τ = σ 11 √ m/ε 11 and pressure in ε 11 /σ 3 11 . The interaction parameters between mixed species ε 12 and σ 12 , are defined as geometric averages ε 12 = √ ε 11 ε 22 and σ 12 = √ σ 11 σ 22 . Interactions are truncated at a cut-off r c = 3.7 and long-range corrections to potential energy and pressure are applied to compensate the effect of truncation to the LJ equation of state [31]. A so large r c is necessary to account correctly for the interactions, when σ 22 is much larger than 1.
The direct calculation of the coexistence temperature T c by either heating a crystal until it melts or cooling a liquid until it freezes would be affected by large hysteresis, due to the finite size of the simulation cell and the presence of the nucleation barrier. To circumvent this issue, we adopt the two-phase method [32], which allows an accurate and direct evaluation of T c at a given pressure in a single simulation. This approach consists of simulating a crystal in contact with its melt in the same simulation cell with periodic boundary conditions ( Figure 1). After equilibration at a chosen temperature, the system with a fixed number of particles N is allowed to evolve at constant pressure and constant enthalpy (NPH). In the NPH ensemble, one of the two phases, depending on the initial temperature, grows at the expenses of the other producing or adsorbing latent heat and the system evolves spontaneously towards the coexistence temperature. Provided that the size of the system is sufficiently large to avoid interactions between solid-liquid interfaces, the T c obtained by this approach is in excellent agreement with that computed by the thermodynamic integration [33].
The starting configuration, with dimension 25 × 25 × 85, containing about 45 · 10 3 particles, was equilibrated using the Berendsen thermostat and barostat [34] at the melting temperature of the single component system, T pure c = 0.694(1), at a pressure of 0.002387 ε 11 /σ 3 11 . The NPH ensemble is enforced by using the Parrinello-Rahman barostat [35], an extension to non-isotropic cell deformation of the Andersen barostat [36]. This algorithm is based on an extended Hamiltonian, in which the cell parameters are additional dynamical variables, and its conserved quantity is the total enthalpy. The pressure is controlled by using a barostat coupling parameter τ p = 30 and the default compressibility, κ T = 0.021, and the angles of the simulation box are constrained. The systems consist of a solid bock made of type-1 (solvent) particles, and a liquid phase containing also type-2 solute particles at the MD simulations are performed using GROMACS [37] and the equations of motion are integrated with a timestep t = 0.005. Production runs in which thermodynamic averages are calculated are 4 × 10 5 steps long, resulting in 2000τ long trajectories. Five independent runs for each parameter pair (size, initial concentration) were carried out. During the NPH simulation, the number of particles in the liquid state changed due to the phase transition. In order to obtain the final concentration of the system, the number of particles in the liquid state was averaged over the last 100 frames (starting at t = 1500 τ in steps of 5) and the concentration computed to c = (number of solutes/average number of solvent particles). Taking the mean value (along with the standard deviation) of all five runs leads to the final result for a given solute parameter pair. The melting temperature was obtained by averaging over the last 500 τ of each independent trajectory.
To determine which solvent particles are in either the liquid or the solid phase, we use Steinhardt's order parameter [38] modified by Wolde [39]. We calculate for each solvent particle i, the complex vector q 6m (i) as follows: for each neighbour j = 1. . .N b (i) of particle i, where N b (i) is the total number of neighbours around i within a chosen cut-off range (we chose the first minimum in the radial distribution function), calculatẽ wherer ij is the normalised distance vector between i and j and Y 6m (r ij ) are spherical harmonics of angular momentum l = 6. The sum of the normalised scalar products leads to q 6 (i), and quantifies the symmetry around particle i. If this value is greater than 0.7, the particle is identified as crystalline. The local order parameter q 6 has proven to work particularly well for the expected face-centred-cubic symmetry in LJ systems [25]. We also performed microcanoncial MD simulations to analyse the structure and thermodynamics of liquid solutions in a simulation box containing about 22 · 10 3 particles. Specifically, we computed the radial distribution function g 11 (r) of the solvent and the static structure factor S 11 (k). Hereon for the sake of simplicity, we refer S 11 (k) and g 11 (r) as S(k) and g(r). S(k) is calculated as the Fourier transform of g(r): We calculated the conformational entropy (S) of the solvent using the excess entropy formalism [40,41]. The total thermodynamic entropy is the sum of the ideal gas entropy S id and the excess entropy S e , The entropy of an ideal gas reads where N is the number of particles of the system, V the volume, k B the Boltzmann constant, β the inverse temperature, h the Planck constant and λ the de Broglie wavelength, The excess entropy is expressed as a multiparticle correlation expansion [40,41], with S n as the n-particle correlation contribution. S 2 , the contribution due to pair correlations between atoms of the same type, can be calculated from the pair distribution function as follows: where the integral is computed from 0 to a radius R max at which g(R max ) 1. Our expansion of the excess entropy is truncated after the pair wise term S 2 . The calculation of S is carried out taking into account the solvent only, therefore the quantities N, g(r) and ρ refer solely to type-1 particles.
Solutes in the binary LJ mixture start to form clusters with increasing solute size and concentration depending on the solute-solvent and solute-solute interactions. The amount of clusters, or the affinity to form clusters, is computed by using a clustering order parameter , defined as = number of clusters number of solutes .
The clusters in the liquid mixture are identified and counted using a nearest-neighbour criterion. For each solute particle, all neighbouring solute particles within a cut-off corresponding to the first minimum in the radial distribution function g 22 (r) are recursively marked as belonging to the same cluster.

Freezing-point depression
The proportionality constant in Equation (1) for an ideal solution is given by [42]: where H is the enthalpy difference between the solid and the liquid phase. The melting enthalpy per particle for our system at P = 0.002387 ε 11 /σ 3 11 is H = 1.12467 ε 11 , in  (1). Considering the differences in some simulation parameters, e.g. the interaction cut-off and the size of the system, this result agrees well with the result obtained from the thermodynamic integration by Agrawal and Kofke [23], T c = 0.687(4). agreement with the previous calculations at similar conditions [31,43]. At T c = 0.694(1) with the gas constant in reduced units, R = 8.3473 · 10 −3 /K, k for an ideal solution in reduced units at the thermodynamic conditions considered here is 3.575(8) · 10 −3 .
Using the two-phase approach, we compute the melting temperature for the LJ binary mixtures described in the previous section. We found that the expected linear dependence of T c on the number, or equivalently molal, concentration in Equation (1) is retrieved for all the considered solvophilic solutes. Examples are shown in Figure 2, in which we plot T c as a function of concentration of solvophilic solutes of two different sizes. The figure clearly shows that the slope of the linear fit is not independent on the size of the solute particles, as predicted by the ideal solvation model, but varies from 0.0003(1) to 0.0039 (2). We can already argue that these systems deviate from the ideal behaviour, with the consequence that k is not solely a colligative property of the solvent. Hence, we will call the proportionality factor in Equation (1) α(σ 22 ) to avoid confusion with the theory of fully diluted systems. It is then worth computing systematically the proportionality factor α(σ 22 ) as a function of the size and the interaction type of solutes, as shown in Figure 3. For solvophilic solutes (Figure 3(a)), for sizes σ 22 ≥ 1.3, the proportionality factor α in T c (c) does not depend on σ 22 within error bars and is fairly close to the k predicted by theory (solid line in Figure 3). On the other hand, for σ 22 ≤ 1.25, α(σ 22 ), is much smaller. α(σ 22 ) increases with σ 22 and saturates at a plateau for 1.3 ≤ σ 22 ≤ 2.0 at α ≈ 0.00334(4), 7% below the theoretical value for the cryoscopic constant calculated from Equation (13).
The linear dependence of T c on the entire concentration range of c up to 2.5% is observed only for solutes with sizes σ 22 ≤ 1.3. For larger solutes, T c is proportional to c only up to a threshold concentration c max . The reason is that at c ≥ c max , the solutes start to cluster and one should correct the definition of number concentration taking into account such effect. Thus, for weakly solvophobic solutes with size σ 22 between 1.3 and 1.6, the slope α(σ 22 ) was obtained from a linear fit up to c max . Larger solutes, σ 22 ≥ 1.7, form clusters even at very low concentrations, making it impossible to estimate α(σ 22 ). The maximum concentration c max was defined by a threshold value of the clustering order parameter of Equation (12). With increasing solute size and concentration, weakly solvophobic solutes form more clusters and decreases. At ≈ 80%, the linear dependence of T c on c breaks down and, depending on the solute size, determines the maximum concentration c max . The maximum freezing point depression obtained for weakly solvophobic solutes (1.03%) in the range of concentrations considered here is significantly lower than that calculated for solvophilic solutes, with σ 22 = 2.5, at the concentration of 2.3%, thus indicating that large solvophilic particles can be used more effectively to control the thermodynamics of phase transitions of solutions. In fact, weakly solvophobic solutes depress T c as efficiently as solvophilic solutes only in the low concentration range, when clustering is avoided. On the other hand, when solutes begin to cluster, the freezing temperature of the solution stops decreasing and  becomes independent on the concentration of solute particles, as shown in Figure 4, where T c (c) is plotted together with the clustering order parameter .
In general, our results show that in contrast with the theory for ideal solutions, size and interaction type of the solute particles affect the freezing-point depression. Deviations from the ideal behaviour expressed in Equation (1) need to be analysed in detail.

Entropy corrections for solid state solution
For 'small' solute particles, we observe a significant reduction of proportionality factors α(σ 22 ) with respect to the ideal cryoscopic constant k. The origin of such deviations from the ideal behaviour is that in deriving Equation (13), one ignores the possibility that incorporating solutes in the crystal to form a solid solution may be thermodynamically favourable. In fact, when solutes are sufficiently small, we have to consider also how the Gibbs free energy of the solid phase changes upon incorporation of solute particles. We calculate the change of the free energy G when changing single component crystal of N LJ particles to a solid solution made of N 1 = N − N 2 solvent particles and N 2 solute particles, at constant temperature and pressure. Here we consider solvophilic solutes, both for the sake of simplicity and because incorporation is more effective in this case. The difference in the free energy of a pure solid solvent and a solid solution is G = U − T S + P V , where U is the internal energy and S the entropy of the system. We neglect the relatively small volume term P V, and compute U in a 24,576 particles system at coexistence. To improve the statistics, a set of different number of solutes, {122, 147, 171, 195, 219, 243, 303, 363, 482} was probed. For solutes with a large σ 22 , a smaller subset of {N 2 } was used because of equilibration issues for large N 2 . While the energy difference is computed by MD, the difference in entropy, S = S ss − S sc is derived as follows. To get the configurational entropy of the solid solution S ss , we consider N crystal lattice sites and replace N 2 of them with solute particles. The number of configurations is With S = k B ln( ), we find for the entropy S ss where we used Stirling's approximation and x 2 = N 2 /N. Since the configurational entropy of the single component crystal is zero, S = S ss . The entropy gain in solid solution S counteracts the difference in internal energy U, making the solid solution thermodynamically more favourable than the single component crystal for solutes smaller than ∼1.18, as it is shown in Figure 5, which displays G as a function of solute size.
This shift in the reference free energy of the solid phase leads to deviations of the proportionality constant α from the ideal behaviour for σ 22 < 1.3 for solvophilic solutes. Also for solute sizes σ 22 ≥ 1.18, although less likely, incorporation in the solid phase is thermodynamically competitive. Figure 6. Clustering of the solution of weakly solvophobic solutes. The cluster ratio is shown as a function of solute size σ 22 and concentration c. = 100% is obtained for a fully dispersed system, = 0% if the system is completely phase-separated. See Equation (12) for the definition of . These results are determined by MD simulations in the microcanonical NVE ensemble of a liquid system of ≈22 · 10 3 particles in a periodic supercell.

Clustering of weakly solvophobic solute
particles Weakly solvophobic solutes tend to cluster with increasing solute size and concentration. Clustering produces large deviations from the ideal freezing point depression behaviour.
Here we analyse clustering of solute particles quantitatively. In Figure 6, the clustering order parameter , computed from Equation (12), is plotted as a function of solute size σ 22 and concentration c. The figure shows that weakly solvophobic solutes stay dispersed if their size is close to the solvent particle size (σ 22 = 1.0), but start to cluster at σ 22 = 1.5 for c > 1.0%, often resulting in a completely phase-separated system. Solvophilic solutes stay dispersed over the entire concentration range considered. The two consequences of clustering are (1) saturation T c as a function of concentration, which hinders the calculation of the proportionality factors α(σ 22 ) for solutes with large σ 22 , (2) less effective depression of T c .
Comparing the clustering order parameter and the dependence of T c on the concentration of solvophobic solutes (Figure 4), for all solute sizes we find that a breakdown of the linear regime and saturation of T c occurs for ∼ 80%. Therefore, we can conclude that ∼ 80% dictates the threshold concentration at which a the maximum T c is attained for a given solute size. Remarkably, rescaling the concentration on the x-axis toc based on either the Van't Hoff factor as in Equation (2) or the number of clusters as obtained from the clustering analysis does not restore a linear dependence of T c onc.
We need to point out that our simulations have been performed on finite periodic systems, and those with the highest degree of clustering tend to phase separate. However, phase separation does not occur over the simulated time scales, therefore we are actually probing out-of-equilibrium cases, which often occur in materials growth and food processing, as well as in the atmosphere.

Static structure factor and configurational entropy
After discussing the deviations from the ideal behaviour of the freezing point depression, we now aim at establishing a direct link between the thermodynamic of the liquid-solid transition in solutions and the changes in the microscopic structure of the solvent upon insertion of solutes. In a seminal paper for the whole field of MD simulations, Hansen and Verlet demonstrated that a LJ liquid crystallises when the first peak of the structure factor S(k) is about 2.85 [31]. It is then reasonable to expect a correlation between the melting temperature of solutions and the amplitude of the first peak of the S(k) of the solvent. In our simulations, the amplitude of the first peak of the S(k) of the pure solvent at T c is ∼2.81, in agreement with more recent estimates for liquids with various interaction potentials and for colloidal particles [44,45]. To probe such correlation, we computed the S(k) of binary LJ solutions in the microcanonical (NVE) ensemble, after equilibrating the systems at the melting temperature of the pure solvent. We considered both solvophilic and weakly solvophobic solutes of four different sizes (1.1, 1.5, 2 and 2.5 σ ) over the whole range of concentrations up to 2.5%. Figure 7 shows that the amplitude of the first peak of the structure factor and the melting temperature of the solutions are indeed unmistakingly correlated, confirming that the mid-range structural modifications induced by the presence of solutes dictate the thermodynamics of phase transitions in solutions. The plot, however, does not include systems that tend to cluster and phase separate ( < 80%), for which the calculation of the S(k) is made problematic by size effects and by the out-of-equilibrium state of these systems.
While this scenario may be compatible with the basic assumptions of the thermodynamics of freezing point depression in ideal solutions, we find that also systems with non-ideal freezing-point depression proportionality factor fall well within this general trend, and obey the Hansen-Verlet criterion. In fact, the changes in the amplitude of the first peak of the S(k) do not depend solely on the concentration but also on the size and the type of interaction of the solute particles. This fact can be appreciated in Figure 8  ability of different solutes to inhibit freezing, as formerly reported in Figure 3. No further significant changes appear in either the radial distribution function or the S(k), except for a peak arising at very small wavevector. The latter peak is an indicator of a long-range structure in the solute, possibly affected by the limited size of the simulated systems.
The amplitude of the first peak in the static structure factor S(k) is sensitive to medium-range structural disorder, and is a qualitative indicator of the configurational entropy of a liquid. The entropy S of the solvent is the dominating factor among the thermodynamic contributions to the free energy and is the leading influence on T c . A quantitative relation between structure and configurational entropy is given by the definition of the excess entropy S e , as the integral of the radial distribution function [40,41,46]. S e computed according to Equation (11), is a pairwise addition to the entropy of the ideal gas, which depends on the density of the liquid as in Equation (9). The truncation to the pairwise term leads to a systematic overestimate of S, which is relatively small at temperatures near the melting point. For example, this error was estimated about 2% for liquid sodium [40]. The entropy thus computed also correlates linearly with T c , especially for solvophilic solutes, as shown in Figure 9. Finite size and equilibration issues, related to the tendency to phase separate of weakly solvophobic solutes with σ 22 > 1.5, hinder the calculation of S e , yielding too large uncertainties.

Conclusions
In summary, utilising MD simulations of two-phase systems, we have investigated the freezing-point depression in model LJ solutions with fairly low concentration of solutes of varying size and two different types of interactions with the solvent, solvophilic and weakly solvophobic. Our results confirm that T c depends linearly on the molal concentration of the solute, as expected for ideal solutions at the limit of infinite dilution, as long as the solute does not tend to cluster and phase separate. We found that when the dispersion of the solute gets smaller than 80% increasing the molal concentration does not affect the freezing temperature of the solution. For this reasons, solvophilic solutes are more efficient to shift T c . In addition, in the range of sizes and concentrations of solutes in which T c is proportional to the molal/number concentration, the proportionality constant is different from the ideal cryoscopic constant, and exhibits a dependence on the size of the solutes and the type of interaction of solutes-solvent interaction. Our analysis supports the classical picture based on the application of Raoult's law on ideal solutions, in which the presence of solutes shifts the free energy of the solvent by increasing its entropy. Deviations from the ideal cryoscopic constant for solutes of small size stem from the possibility of being incorporated in the solid phase. Even though the formation of solid solutions may shift the free energy of the solid by a G one order of magnitude lower than the free energy changes occurring in the solvent, this effect is sufficient to produce anomalous proportionality factors between T c and c.
We also show that T c is strictly linked to the structural changes of the solvent induced by the presence of solutes. The observed correlations between either the mediumrange structure or the entropy of the solvent and T c then suggest that solutes affect the structure of solvents, thus increasing their entropy and stabilising the liquid phase.