Classical Nucleation Theory for the Crystallization Kinetics in Sheared Liquids

TL;DR

This paper extends classical nucleation theory to sheared liquids, demonstrating through simulations that shear influences crystallization rates via elastic and interfacial work, with size-dependent effects on nucleation energetics.

Contribution

It introduces a modified CNT framework incorporating elastic and interfacial work under shear, validated by numerical simulations of hard spheres.

Findings

Shear impacts nucleation rate through elastic work proportional to droplet volume.

Local shear modulus inside droplets is about half of the bulk value.

Interfacial work change is significant for small nuclei but negligible for larger ones.

Abstract

While statistical mechanics provides a comprehensive framework for the understanding of equilibrium phase behavior, predicting the kinetics of phase transformations remains a challenge. Classical nucleation theory (CNT) provides a thermodynamic framework to relate the nucleation rate to thermodynamic quantities such as pressure difference and interfacial tension through the nucleation work necessary to spawn critical nuclei. However, it remains unclear whether such an approach can be extended to the crystallization of driven melts that are subjected to mechanical stresses and flows. Here, we demonstrate numerically for hard spheres that the impact of simple shear on the crystallization rate can be rationalized within the CNT framework by an additional elastic work proportional to the droplet volume. We extract the local stress and strain inside solid droplets, which yield size-dependent values for the shear modulus that are about half of the bulk value. Finally, we show that for a complete description one also has to take into account the change of interfacial work between the strained droplet and the sheared liquid. From scaling reasons, we expect this extra contribution to dominate the work formation of small nuclei but become negligible compared to the elastic work for droplets composed of a few hundreds particles.

Figures

Equations

$$\Delta F_0(V_s)=-\Delta P{V}_s+\Phi (V_s).$$

$$W_e(V_s)=\frac{{\sigma }_s^2}{2G(V_s)}{V}_s$$

$$\Delta P-\frac{{\sigma }_s^2}{2G_{\text{eff}}}=\frac{2{\Gamma }_{\ast }}{R_{\ast }}$$

$$G_{\text{eff}}=G_{\ast }{\left[1-\frac{R_{\ast }}{3G_{\ast }}{\frac{\partial G}{\partial R}|}_{R_{\ast }}\right]}^{-1}.$$

$$W({\sigma }_s;R_{\ast })=\frac{16\pi {\Gamma }_{\ast }^3}{3(\Delta P-\frac{{\sigma }_s^2}{2G_{\text{eff}}}{)}^2}$$

$$a_W^{\text{CNT}}=\frac{1}{G_{\text{eff}}\Delta P}.$$

$$a_N^{\text{CNT}}=\frac{3}{2}{a}_W^{\text{CNT}},$$

$$a_Z^{\text{CNT}}=-a_W^{\text{CNT}},$$

$$a_k^{\text{CNT}}=-\Delta F_0{a}_W^{\text{CNT}}.$$

$${\dot{\mathbf{r}}}_i=-\frac{D_0}{k_{B}T}{\nabla }_{i}U+\dot{\gamma }{y}_i{\mathbf{e}}_x+\sqrt{2D_0}{\xi }_i$$

$$q_{l,m}(i)=\frac{1}{N_n(i)}\sum _{j=1}^{N_n(i)}{Y}_{l,m}({\theta }_{i,j},{\phi }_{i,j}),$$

$$d(i,j)=\frac{{\sum }_{m=-l}^l{q}_{l,m}(i)q_{l,m}^{\ast }(j)}{({\sum }_{m=-l}^l|q_{l,m}(i){|}^2{)}^{1/2}({\sum }_{m=-l}^l|q_{l,m}(j){|}^2{)}^{1/2}}$$

$${\mathbf{r}}_i^{\gamma }={\mathbf{r}}_i+\gamma y_i{\mathbf{e}}_x$$

$$P_B(n)=\frac{1}{2}(1+\mathrm{erf}[\sqrt{\pi }{z}_c(n-n_c)]),$$

$${\hat{\sigma }}_i=\frac{1}{2v_i}\sum _{j\ne i}\frac{x_{ij}{y}_{ij}}{r_{ij}}{u}^{\prime }(r_{ij}),$$

$${\sigma }_{\text{sub}}=\frac{1}{{\sum }_{i\in C}{v}_i}\sum _{i\in C}{v}_i{\hat{\sigma }}_i,$$

$$\rho (r)=\frac{{\rho }_l+{\rho }_s}{2}+\frac{{\rho }_l-{\rho }_s}{2}\tanh \left(\frac{r-R_0}{w}\right),$$

$$D_{\text{min}}^2(i,t_0\to t)=\sum _{j=1}^n[({\mathbf{r}}_j(t)-{\mathbf{r}}_i(t))-\mathbf{D}\times ({\mathbf{r}}_j(t_0)-{\mathbf{r}}_i(t_0)){]}^2.$$

$$\Delta F(n;\dot{\gamma })=F(n;\dot{\gamma })-F_0(n)=W_e(n;\dot{\gamma })+\Delta W$$

$$W_e(n;\dot{\gamma })=\frac{1}{2}{\sigma }_s(n){\gamma }_s(n)V_s(n)$$

$$\Delta F(n;\dot{\gamma })=F_0(n)+W_e(n;\dot{\gamma })+\Delta W(n;\dot{\gamma }).$$

$$D_0=\frac{3}{8\sqrt{\pi }}{\left(\frac{k_{B}T}{m}\right)}^{1/2}\left(\frac{1}{\rho d^2}\right)$$

$$\exp (-\Delta \mathcal{F})=⟨\exp (-\beta w_{\tau })⟩,$$

$$w_{\tau }=\sum _{t=0}^{\tau -1}[H_{{\lambda }_{t+1}}({\omega }_t)-H_{{\lambda }_t}({\omega }_t)],$$

Full text

I Introduction

Classical nucleation theory (CNT) is an establish thermodynamic framework that helps understanding phase formation and interpreting experiments and numerical simulations. Notable examples are the modeling of ice nucleation rates (impacting the understanding of our climate Zobrist et al. (2007); Hoose et al. (2010); Laksmono et al. (2015)) and the estimation of interfacial tension Knott et al. (2012); Sanz et al. (2013); Espinosa et al. (2016a, b). Originally developed to study systems that are prepared in thermal equilibrium, there have been several attempts to extend CNT to systems driven into a non-equilibrium steady state Desre and Yavari (1990); Kashchiev (1972); Reguera and Rubi (2003a, b); Mokshin et al. (2013); Mura and Zaccone (2016). Of particular interest is crystallization in the presence of mechanical stresses and flows Reguera and Rubi (2003b); Mokshin et al. (2013); Mura and Zaccone (2016). How such driving forces can control not only the nucleation kinetics but also the structure of the newborn solid phase remains poorly understood.

Already for relatively simple models, such as liquids and colloidal suspensions in which particles interact through soft or hard-core repulsions, the effect of flow on the nucleation kinetics is far from trivial and heavily depends on the shear strength as well as the degree of supersaturation (or cooling) Ackerson and Clark (1981); Haw et al. (1998); Blaak et al. (2004a); Holmqvist et al. (2005); Cerdà et al. (2008); Wu et al. (2009); Lander et al. (2013); Shao et al. (2015); Mao et al. (2015); Ruiz-Franco et al. (2018). Aside from shear-induced order (particle layering) at high strain rates Ackerson and Pusey (1988); Xue and Grest (1990); Panine et al. (2002); Nikoubashman et al. (2011); Besseling et al. (2012) one finds that supercooled liquids crystallize basically via the same activated nucelation process as in the quiescent regime: The system remains in the disordered melt until a rare fluctuation leads to a sufficiently large critical nucleus that grows spontaneously. This suggests to write the crystallization rate as $k=\kappa e^{-W}$ with work $W$ (in units of the thermal energy $k_{B}T$) to escape from the metastable state and kinetic prefactor $\kappa$. The theoretical challenge is to find expressions for $W$ (and $\kappa$) as a function of the relevant system parameters.

In a pioneering numerical work Blaak et al. (2004a, b) studying the crystallization of colloidal particles it has been shown that shear flow can suppress crystallization. This has been rationalized through a significant increase of $W$ as a function of the strain rate. However, other experimental and numerical studies have reported the opposite behavior with shear induces order Holmqvist et al. (2005); Mokshin and Barrat (2010); Mokshin et al. (2013); Wu et al. (2009); Shao et al. (2015) and even an optimal strain rate for which the crystallization kinetics is fastest Cerdà et al. (2008); Mokshin et al. (2013); Lander et al. (2013). More recently, we have probed the crystallization of hard spheres as a function of the strain rate and the packing fraction Richard and Speck (2015), the latter serving as a control parameter for the degree of supersaturation. We found a crossover from shear-induced suppression to shear-induced enhancement: At low packing fractions, the rate is dominated by the nucleation barrier tending to increase with the strength of the shear flow. As for soft spheres, small clusters are more likely to dissolve and one finds larger critical nucleus sizes. In contrast, for dense suspensions the (quiescent) nucleation barrier effectively vanishes and the rate is controlled by the kinetic prefactor $\kappa$. The latter strongly follows the particle dynamics and drops close to the glass transition due to caging effects. In this regime external flow enhances the particle diffusion and allows for a better exploration of configuration space, accelerating the formation and growth of solid nuclei. Continuing to increase the strain rate, clusters start to break-up and the rate again is controlled by the nucleation work of critical nuclei.

In this paper, we numerically study the extension of CNT to sheared liquids. As a well-studied testbed, we choose monodisperse hard spheres. Recently, Mura and Zaccone have put forward the idea that the quiescent nucleation work is to be augmented by a reversible elastic work to stress the critical nucleus Mura and Zaccone (2016). Such an extension has already been proposed to model the steady coexistence of a strained solid with its sheared melt Butler and Harrowell (2002, 2003), but was found inconclusive with respect to the existence of a chemical potential out of equilibrium. Here, we numerically test this scenario and demonstrate that the change of nucleation kinetics of sheared hard spheres can indeed be captured through an additional elastic work. Crucially, for quantitative predictions one needs to take into account the diminished density and the elevated shear stress inside finite solid droplets.

The manuscript is organized as follows: We first discuss the theoretical basis for extending CNT in Sec. II before providing details of the model system and the simulations in Sec. III. In Sec. IV, we first show that employing bulk values for the shear modulus and equating the solid stress with the hydrodynamic liquid stress does not yield reasonable predictions. We then present a method to access the local shear stress and shear strain of finite solid droplets, which gives access to the shear modulus and shows that the stress inside the droplet is elevated. Finally, we reconstruct the elastic work as a function of the droplet size and highlight the presence of an extra surface work for small droplets.

II Theory

II.1 Classical nucleation theory under shear

We aim to predict the isothermal nucleation rate $k$ per volume with which a critical solid droplet spontaneously forms within the sheared melt (Fig. 1). The two main control parameters are the liquid density $\rho_{l}$ and the imposed strain rate $\dot{\gamma}$. The latter determines the shear stress $\sigma_{l}=\eta\dot{\gamma}$ in the liquid acting upon the droplet, with $\eta(\dot{\gamma},\rho_{l})$ denoting the shear viscosity of the liquid at a particular density. Note that we explicitly take into account a dependence of $\eta$ on $\dot{\gamma}$ since dense liquids might undergo shear thinning. The solid droplet is characterized by its volume $V_{s}$, its density $\rho_{s}(V_{s})$, and its shear modulus $G(V_{s})$, both of which depend on the droplet size. Throughout, we follow Gibbs idea of a dividing surface separating the solid droplet from the liquid, both of which are modeled as homogeneous systems.

Without shear flow ($\dot{\gamma}=0$), the nucleation rate is dominated by the reversible nucleation work to reach the transition state, i.e., a critical solid droplet of volume $V_{s}$. This nucleation work is given by the free energy difference

[TABLE]

The first term is proportional to the volume and captures the free energy gained through creating space for the droplet. The second term $\Phi(V_{s})=\Gamma A$ is the excess interfacial free energy given as the product of the droplet surface $A\sim V_{s}^{2/3}$ and the interfacial tension $\Gamma$. In principle, $\Gamma$ again depends on the droplet size Richard and Speck (2018a). The thermodynamic driving force for nucleation is the difference $\Delta P=P_{s}-P_{l}$ of pressure $P_{s}$ between the inside of the solid droplet and the ambient liquid pressure $P_{l}$ at the same liquid chemical potential, $\mu_{s}=\mu_{l}$. It is assumed that the surrounding liquid is stress-free and therefore there is no elastic work on the nucleus entering the nucleation work.

Turning on the shear flow with $\dot{\gamma}\neq 0$, the system is steadily driven away from equilibrium into a non-equilibrium steady state and thus constantly dissipates heat. Strictly speaking, there is no thermodynamic potential anymore that determines the behavior of the system. To proceed, we make the three following assumptions: First, in our modeling we neglect the dissipation due to shearing the liquid. We treat the droplet as an inclusion in a stressed medium (the surrounding liquid) but ignore the “housekeeping” work that needs to be spend to keep the medium at a given shear stress $\sigma_{l}$ Speck and Seifert (2009); Gerloff and Klapp (2018). We do, however, consider the excess work $W$ required to form the critical solid droplet. Second, we assume that the nucleation rate is still determined by this excess work, $k=\kappa e^{-W}$, i.e., the droplet still emerges due to a spontaneous thermal fluctuation. Put differently, the liquid acts as a heat reservoir, the fluctuations of which are still characterized by a (possibly effective) temperature. In contrast to the quiescent liquid, however, there is an additional elastic work

[TABLE]

to create the droplet with shear stress $\sigma_{s}$ Landau and Lifshitz (1959); Butler and Harrowell (2003); Mura and Zaccone (2016). Note that this work is not compensated by a reduction in free energy of the liquid since the external work to maintain $\sigma_{l}$ is immediately dissipated. Hence, the total nucleation work now reads $W=\Delta F_{0}+W_{e}$. Third, we assume that the droplet undergoes a pure shear transformation with strain $\gamma_{s}=\sigma_{s}/G$.

We now restrict ourselves to spherical droplets with radius $R$, volume $V_{s}=\frac{4\pi}{3}R^{3}$, and area $A=4\pi R^{2}$. The mechanical equilibrium condition ($\partial W/\partial R=0$) yields

[TABLE]

at the surface of tension $R_{\ast}$ defined through $\partial\Gamma(R)/\partial R|_{R_{\ast}}=0$, where $\Gamma_{\ast}=\Gamma(R_{\ast})$ and the shear modulus $G_{\ast}=G(R_{\ast})$ are evaluated at the critical droplet radius $R_{\ast}$. Here we have introduced the effective shear modulus

[TABLE]

Eliminating volume and area of the critical droplet, the nucleation work then takes the customary CNT form

[TABLE]

known from quiescent nucleation but with $\Delta P$ replaced by the effective driving force $\Delta P-\frac{\sigma_{s}^{2}}{2G_{\text{eff}}}$, which is reduced due to the additional elastic work required to deform the solid nucleus.

For small shear stress $\sigma_{s}$ in the linear response regime, we expand the nucleation work in the form $W(\sigma_{s})\approx\Delta F_{0}(1+a_{W}\sigma_{s}^{2})$, where the CNT prediction for the response coefficient $a_{W}$ depends on the liquid and solid properties through $\Delta P$ and $G_{\text{eff}}$ as

[TABLE]

The same expansion can be performed for the critical nucleus size $N_{c}=\tfrac{4\pi}{3}\rho_{s}R_{\ast}^{3}\approx N_{0}(1+a_{N}\sigma_{s}^{2})$ and the Zeldovich factor $Z_{c}=\sqrt{\frac{|F^{\prime\prime}(q_{c})|}{2\pi}}\approx Z_{0}(1+a_{Z}\sigma_{s}^{2})$, where $N_{0}$ and $Z_{0}$ are the critical nucleus size and Zeldovich factor in the quiescent limit $\dot{\gamma}=0$, respectively. Moreover, for the rate we find $\ln k\approx\ln k_{0}+a_{k}\sigma_{s}^{2}$ with $k_{0}$ the nucleation rate at vanishing stress. From the CNT expression for the work Eq. (5), we thus find the predictions

[TABLE]

For notational convenience, we drop the subscript “eff” for the shear modulus in the following.

III Methods

III.1 Model

We perform non-equilibrium Brownian dynamics (BD) simulations of a mono-component hard-sphere fluid. Particles interact through the pairwise Weeks-Chandler-Anderson potential $u(r)=4\epsilon[(\alpha/r)^{12}-(\alpha/r)^{6}+1/4]$ for $r<2^{1/6}\alpha$. Simulations are done in the canonical ensemble (NVT) with a fixed number of particles $N$, volume $V$, and temperature $T$. The system is composed of $N=5000$ particles if not mentioned otherwise.

As indicated in Fig. 1, the direction of the shear flow is set along $\mathbf{e}_{x}$, the flow gradient and vorticity along $\mathbf{e}_{y}$ and $\mathbf{e}_{z}$, respectively. The coupled equations of motion read

[TABLE]

with $D_{0}$ the bare translational diffusion coefficient, $\mathbf{\xi}_{i}$ Gaussian white noise, and $-\nabla_{i}U$ is the conservative force acting on particle $i$. Additionally, we employ Lees and Edwards periodic boundary conditions Allen and Tildesley (2017). The potential strength $\epsilon$ is set to $40k_{B}T$, with $k_{B}$ the Boltzmann constant. Throughout all our simulations, we scale lengths by $\alpha$, times by $\alpha^{2}/D_{0}$, and energies by $k_{B}T$. The equations of motion are integrated with time step $10^{-5}$. All results are then expressed in hard-sphere units where lengths are measured in units of the effective diameter $d=1.097\alpha$ and the Brownian time is defined as $\tau_{B}=d^{2}/D_{0}$. The packing fraction $\phi$ of the system is given by $\phi=\pi Nd^{3}/(6V)$. Further details concerning this mapping can be found in Refs. Filion et al. (2011); Richard and Speck (2018b).

III.2 Bond-orientational order

We monitor the degree of crystallinity in our system through the local bond-orientational order parameter Steinhardt et al. (1983); Lechner and Dellago (2008)

[TABLE]

which is evaluated for particle $i$, where $Y_{l,m}(\theta,\varphi)$ are spherical harmonics and $N_{n}$ is the number of neighbors within distance $r_{ij}<1.5\alpha$. We construct a bond network through the scalar product

[TABLE]

using $l=6$ with $d(i,j)>0.7$ defining a bond. Finally, a particle is defined as “solid-like” if the number of bonds $\xi\geq 9$, and clusters are constructed from mutually bonded solid-like particles.

III.3 Shear viscosity and shear modulus

In a homogeneous liquid, we monitor the fluctuating shear stress $\hat{\sigma}$ at constant strain rate $\dot{\gamma}$ and deduce the shear viscosity $\eta=\langle\hat{\sigma}\rangle/\dot{\gamma}$, where the average $\langle\cdots\rangle$ here involves only configurations that have less than 5% “solid-like” particles. Since dense suspensions undergo shear thinning at finite strain rates, we first evaluate $\eta$ as a function of $\dot{\gamma}$ and fit the flow curve using a Carreau model Bruus (2008), giving us access to the zero shear rate viscosity $\eta_{0}$. Beyond $\rho_{l}\simeq 1.06$ (and for $N=5000$), monodisperse hard spheres crystallize within $<10\tau_{B}$, preventing the correct estimation of $\eta$. To circumvent this limitation, we have performed additional simulations with a $\Delta=5\%$ Gaussian polydispersity. Results are shown in Fig. 2(a) for $\dot{\gamma}=0.0,0.036$, and $0.084$. At vanishing stress $\dot{\gamma}\to 0$, the viscosity diverges approaching $\rho_{l}\simeq 1.1$ ($\phi\simeq 0.576$) marking the onset of dynamical arrest, which is well modeled by a Vogel-Fulcher-Tammann (VFT) form Debenedetti and Stillinger (2001). For $\dot{\gamma}>0$, we observe a drop of the viscosity indicating shear thinning.

To investigate the role of an elastic work, we require access to the stress $\sigma_{s}$ and the shear modulus $G$ inside critical droplets as a function of droplet size $R_{\ast}$. As reference, we first determine the bulk shear modulus $G_{\infty}=G(R_{\ast}\to\infty)$ as a function of the solid density $\rho_{s}$. To this end, we compute the bulk shear modulus of a face-centered-cubic (fcc) crystal closely following Ref. 46. In these simulations, particle positions are subjected to the affine transformation

[TABLE]

with strain $\gamma$. The shear modulus is computed from the initial slope of the average stress $\langle\hat{\sigma}\rangle_{\gamma}$ plotted against the imposed strain $\gamma$. In Fig. 2(b), we compare our shear modulus with values of true hard spheres. We find a very good agreement between the two systems, supporting the validity of the mapping procedure.

IV Results

IV.1 Extracting nucleation rate and work

We employ the same framework as presented in Ref. 40. We prepare initial configurations at packing fraction $\phi$ without any solid particles using the algorithm developed by Clarke and Wiley Clarke and Wiley (1987), where the non-overlapping distance between particles is chosen to be equal to the effective diameter $d_{\text{eff}}$. We then harvest $300$ independent trajectories. As order parameter we employ the size $n$ of the largest cluster, and trajectories are terminated when the system crosses the barrier and reaches an absorbing boundary at $n_{b}$. Specifically, we choose $n_{b}=400$, which is four times larger than the largest critical size $n_{c}$ found in this work. This construction enables to sample a non-equilibrium steady-state with a net current $j$ of droplets flowing from the metastable melt towards $n_{b}$. As shown in Refs. Manuel et al. (2015); Richard and Speck (2018b) and below, it allows to consistently extract nucleation barriers $F(n)$ through linking the non-equilibrium distribution $P^{+}(n)$ and the splitting probability $P_{B}(n)$ that a configuration at $n$ will commit to the solid phase and reach $n_{b}$. For large barriers, a quadratic expansion of the free energy yields the expression

[TABLE]

which is used to extract critical nuclei sizes $n_{c}$ and Zeldovich factors $z_{c}$. Note that we explicitly distinguish the variables $N_{c}$ and $Z_{c}$ appearing in the CNT expressions from $n_{c}$ and $z_{c}$. The latter are computed using the bond-order parameter described in Sec. III.2 and thus depend on the set of parameters employed to construct the bond network.

Fig. 3(a) shows how $P_{B}$ changes with the imposed strain rate $\dot{\gamma}$ at $\phi\simeq 0.542$, which is close to the melting point located at $\phi_{m}\simeq 0.543$. Progressively increasing the strain rate $\dot{\gamma}$, we first observe a shift of $P_{B}$ towards larger $n$, indicating that small clusters are likely to dissolve under shear. Second, we find a systematic broadening of the splitting probability, implying a flattening of the barrier at $n\simeq n_{c}$, i.e., a smaller Zeldovich factor. Both behaviors are in qualitative agreement with Eq. (7). In Fig. 3(b), we show at $\dot{\gamma}\simeq 0.06$ consecutive unsuccessful nucleation events, whereby nuclei larger than the quiescent critical nucleus size $n_{c}(\dot{\gamma}=0)$ fully dissolved. Furthermore, we have computed the mean first passage time (MFPT) $\tau^{+}(n)$ to reach a given size $n$ starting from the metastable liquid. The MFPT is inversely proportional to $P_{B}(n)$ scaled by the nucleation time $\tau_{x}=1/j$. In Fig. 3(b), we confirm this connection by comparing $\tau^{+}(n)$ with $P_{B}(n)$ scaled by $\tau_{x}$, which is the value at which $\tau^{+}(n)$ plateaus.

Having collected a set of trajectories, we can now compute the stationary probability distribution $P^{+}(n)$ to observe a configuration with a droplet of size $n$. As shown in Refs. Manuel et al. (2015); Richard and Speck (2018b), we can reconstruct the actual distribution $P(n)=P^{+}(n)/[1-P_{B}(n)]$. In our simulations, we also compute the average number of clusters of size $n$, which allows to correct $P(n)$ for small clusters Richard and Speck (2018b). We interpret $F(n)\sim\ln P(n)$ as an effective free energy governing the nucleation kinetics, from which we extract the nucleation work $W=\Delta F$ as the height of the barrier. Many more details of this procedure can be found in Refs. Lundrigan and Saika-Voivod (2009); Leitold and Dellago (2016); Richard and Speck (2018b). In Fig. 3(d), we plot the profile $F(n)$ for several strain rates $\dot{\gamma}$ at the same packing fraction. We observe that the top of the barrier moves to larger values $W$ and larger critical sizes $n_{c}$, and that $F(n\simeq n_{c})$ flattens in agreement with a decrease of the Zeldovich factor seen from the broadening of $P_{B}$ and $\tau^{+}$. For the largest strain rate, we find a barrier increase of about $3k_{B}T$, which is consistent with the two orders of magnitude increase in the nucleation time seen in Fig. 3(c).

IV.2 Response coefficients based on liquid shear stress and CNT predictions

We can now extract the nucleation work $W$ and nucleation rates $k=1/(\tau_{x}V)$ as well as critical sizes $n_{c}$ and Zeldovich factors $z_{c}$. For moderate shear rates, plotting these quantities as a function of the square of the solid shear stress $\sigma_{s}^{2}$ would disclose the derived response coefficients in Eq. (6) and Eq. (7). Unfortunately, defining and measuring a mechanical shear stress within finite droplets is non-trivial. As a first step to assess the impact of shear, we will plot results as a function of the ambient liquid stress $\sigma_{l}=\eta\dot{\gamma}$. The relation $\sigma_{s}=\sigma_{l}$ has been proposed in Ref. 13 but we will show below that it does not hold for finite droplets.

As already located in our previous work Richard and Speck (2015), we expect a change from crystallization suppression to enhancement around $\phi\simeq 0.56$. Here, we focus on the regime where shear suppresses crystallization and have applied our methodology for various packing fractions ranging from $\phi\simeq 0.539$ to $\phi\simeq 0.553$. In Fig. 4, we plot $k$, $W$, $n_{c}$, and $z_{c}$ against $\sigma_{l}^{2}$. As we increase $\sigma_{l}$, both the rate and Zeldovich factor decrease, whereas the nucleation work and critical nucleus size increase, qualitatively in agreement with predictions. At large supersaturations, the impact of the flow on the nucleation kinetics is minimal, consistent with our previous study Richard and Speck (2015). In appendix A, we compare Brownian dynamics simulations with our earlier work employing molecular dynamics and show that both dynamics can be mapped onto each other. It shows that the local dynamical rule has little effect on the nucleation kinetics (at least in the linear response regime). It is worth mentioning that the same observation has been made for the crystallization of soft spheres interacting via the Yukawa potential Blaak et al. (2004a). Continuing to increase the shear stress inside the liquid, we observe a deviation from the linear scaling.

Linear regressions of the data in Fig. 4 gives us access to the response coefficients $\{a_{k},a_{W},a_{n},a_{z}\}$ as a function of the packing fraction $\phi$. Again, we make explicit the difference between the $a_{N}$ and $a_{Z}$ derived from the thermodynamic modeling and the computed $a_{n}$ and $a_{z}$ based on the local bond order. In Fig. 5, we plot $-a_{k}/\Delta F_{0}$, $a_{W}$, $2a_{n}/3$, and $-a_{z}$, which should collapse onto a single curve following the CNT prediction in Eq. (7). We first observe that $-a_{k}/\Delta F_{0}$ does indeed fall onto $a_{W}$, indicating that the nucleation rate can be reasonably modeled by the increase of the nucleation work. Interestingly, we find $2a_{n}/3$ and $-a_{z}$ collapsing on top of each other, but with values about two times larger than $a_{W}$. A discrepancy between $a_{W}$ and $2a_{n}/3$ must stem from a non-linear relation between $N_{c}$ and $n_{c}$ since a simple rescaling would leave $a_{n}$ invariant. Such a non-linear relation can be attributed to the fact that we are considering critical droplets composed of hundred particles and less, which are thus mainly composed by particles located at the interface where the identification of “solid-like” particles is somewhat ambiguous.

Finally, we provide in the inset of Fig. 5 a comparison between $a_{W}$ and $a_{W}^{\infty}=1/(G_{\infty}\Delta P_{\infty})$. The latter assumes that critical droplets have the same density $\rho_{s}^{\infty}$ than a bulk solid crystal at the same ambient liquid chemical potential $\mu_{l}$. We take the pressure difference $\Delta P_{\infty}$ from the bulk equations of state and $G_{\infty}$ from our parametrization of the shear modulus shown in Fig. 2(b). We find that both $a_{W}$ and $a_{W}^{\infty}$ decrease exponentially as a function of $\phi$, which is consistent with the sharp increase of the shear modulus beyond the melting point [cf. Fig. 2(b)]. The main observation, however, is that $a_{W}$ is about two orders of magnitude larger than the bulk prediction $a_{W}^{\infty}$. Although it is known that properties of small droplets deviate from bulk quantities, it is unlikely that such a gap can be explained by a hundred times smaller effective shear modulus. In the following, we will present a method to directly access the stress $\sigma_{s}$ and strain $\gamma_{s}$ of critical droplets and demonstrate that the corrected response coefficient $a_{W}$ agrees with the CNT prediction.

IV.3 Seeding of droplets under shear

To generate configurations of large critical droplets we now turn to a different type of simulations, namely the seeding of droplets Sanz et al. (2013); Espinosa et al. (2016a); Richard and Speck (2018a). Starting configurations are prepared in the same way as described in our previous work Richard and Speck (2018a). To avoid finite size effects, we now study larger systems with $N=40,000$ particles. For each packing fraction, strain rate, and seed size, we generate $20$ fleeting trajectories that are terminated when either $n<n_{a}=10$ or $n>n_{b}=1000$ is fulfilled, see Fig. 6(a). We then compute the probability $P_{B}$ that a run crosses $n_{b}$ without coming back to $n_{a}$, which gives us an estimate for the critical nucleus for which $P_{B}\simeq 1/2$ (cf. Eq. (14)). As a consistency check, we show in Fig. 6(b) the evolution of the total shear stress in the system. We find $\sigma$ reaching a steady-state value for $t\simeq 0.5-1\tau_{B}$, which is a negligible relaxation time compared with the typical fleeting time. Hence, our seeding preparation does not alter the extraction of $n_{c}$. In Fig. 6(c), we show our determination of the critical nucleus size $n_{c}$ as a function of the square of ambient shear stress $\sigma_{l}^{2}$. As found in Sec. IV.2, we observe a linear behavior at small driving and we recover for the largest packing fraction ($\phi=0.542$) our previous results from the direct (unseeded) simulations. Linear regressions allow us to extract new estimates for $a_{n}$ until $\phi\simeq 0.52$ and to compare them with $a_{W}$ and $a_{W}^{\infty}$. Clearly, no change of behavior can be seen although we are now probing critical droplets composed of several hundred particles.

IV.4 Extracting local shear stress and strain

We are now proposing a methodology to access the solid density $\rho_{s}$, shear stress $\sigma_{s}$, and strain $\gamma_{s}$ of finite droplets, which ultimately will give us an estimate for the shear modulus $G$. To gather the necessary statistics, in this subsection we consider not only critical droplets but all droplets of a given size $n$.

IV.4.1 Local density and shear stress

Using configurations generated via the seeding method, we are now able to extract density and stress profiles. To do so, we define for each particle $i$ its microscopic local shear stress $\sigma_{i}$ by Egami et al. (1980)

[TABLE]

where $v_{i}$ is the Voronoi volume of the particle $i$. We can evaluate the average shear stress inside a given sub-volume $C$ via a weighted sum as

[TABLE]

which we use to compute the radial stress profile $\sigma(r)$ with bin width $\Delta r=1.5d$ by measuring the particle positions with respect to the center of mass of the droplet. The radial density profile $\rho(r)$ follows simply from counting the number of particles in each bin. Profiles at a particular size $n$ are computed from a minimum of $400$ configurations and are taken with the criteria $n_{i}\in[n-5\%,n+5\%]$, $n_{i}$ being the size of the largest droplet in the system. In Fig. 7(a), we show such density profiles for $\phi\simeq 0.542$ and nucleus sizes $n=80,250,500,$ and $1000$. Note that results for $n=80$ are extracted from configurations generated by direct simulations. As already discussed extensively in our previous work Richard and Speck (2018a), we observe a gradual increase of the density inside the solid droplet as its size grows. Profiles can be well modeled by the mean-field expression

[TABLE]

where $\rho_{s}=\rho(0)$ is the density at the center of the solid droplet, $\rho_{l}$ is the density of the surrounding liquid, $R_{0}$ is the radius at half maximum (different from $R_{\ast}$), and $w$ the interfacial width.

For $n=80$, we find a solid density $\rho_{s}\simeq 1.11$ which would correspond (assuming bulk behavior) to a shear modulus $G\simeq 66$. On the other hand, the bulk equation of state yields $\rho_{s}^{\infty}\simeq 1.154$ and thus a significantly larger shear modulus $G_{\infty}\simeq 99$. This result marks our first finding to explain the discrepancy between $a_{W}$ and $a_{W}^{\infty}$. A smaller solid density results on one side in a smaller pressure difference between the two phases, thus decreasing the nucleation driving force $\Delta P$. On the other hand it results in a smaller shear modulus, increasing the elastic contribution to the nucleation work.

In Fig. 7(b), we show shear stress profiles from the same set of data for several droplet sizes $n$. We find that stress values at the center of droplets are about $3$ to $7$ times larger than the ambient shear stress $\sigma_{l}=\eta\dot{\gamma}$, to which all stress profiles decay. As observed in our density profiles, we find a gradual increase of the solid stress $\sigma_{s}=\sigma(0)$ as a function of the nucleus size $n$. This change is consistent with the increase of the density at the interface that locally increases the viscosity and thus the shear stress acting upon the nucleus. In fact, stress profiles can be well fitted by the same mean-field expression Eq. (17) as for the density profiles. This result marks our second finding, namely that the elastic work to form a droplet is going to be much larger than the prediction made using $\sigma_{l}$ in Eq. (5) (already about $10$ times larger for $n=80$). It also explains why one cannot extract the response coefficient by plotting nucleation works against $\sigma_{l}^{2}$.

IV.4.2 Local shear strain

Having evaluated the shear stress inside small droplets, we would like to characterize which elastic deformation “solid-like” particles have undergone. It is important to point out that there is not a unique and well-defined way to estimate local strains Falk and Langer (1998); Goldhirsch and Goldenberg (2002); Graner et al. (2008); Tsamados et al. (2009). Here, we adopt a method developed in the context of the rheology of glasses Falk and Langer (1998), which is based on an adjustable strain tensor that fits best the actual particle displacements over some time interval. In our methodology, we harvest configurations with strained droplets, see snapshot in Fig. 8(a), and switch off the imposed shear flow. Running short trajectories from these configurations, we can follow the elastic relaxation of the shear stress inside the solid phase and fit the local strain that particles undergo to release such a stress. In practice, we select a minimum of $400$ configurations with the same criteria as used for the density profiles. We run from them one trajectory of length $5\tau_{B}$. In Fig. 8(b), we show for $n=1000$ the average solid stress taken for $r<2d$ and ambient liquid stress taken for $r>3R_{0}$. We find that the stress, both in the solid and liquid, relaxes fast with relaxation times of less than $3\tau_{B}$ and $1\tau_{B}$, respectively.

During each run, we evaluate the local displacement of a particle $i$ between a reference configuration at the starting time $t_{0}$ and a time $t$ through the deviatoric strain Falk and Langer (1998) defined as

[TABLE]

Here, the sum runs over the $n$ closest neighbors of particle $i$, which are determined through a Voronoi tessellation at $t_{0}$. The deformation tensor $\mathbf{D}$ is determined by minimizing $D^{2}_{\text{min}}$. Because droplets can undergo body rotations, we evaluate the symmetric Green strain tensor $\mathbf{E}=\frac{1}{2}(\mathbf{D}^{T}\mathbf{D}-\mathbf{1})$. Finally, we can inspect the off-diagonal $E_{xy}$ to quantify shear deformations. In Fig. 8(c), we plot the distribution $P(2E_{xy})$ at $t=5\tau_{B}$ for particles within a sphere of a radius $R=2d$ from the center of mass of the droplet. We observe a Gaussian distribution centered at negative values indicating that particles have undergone a reverse shear transformation with average strain $\gamma=-2\langle E_{xy}\rangle$. In Fig. 8(d), we show the relaxation of the average strain in both the solid and liquid using the same criteria as for the stress shown in Fig. 8(b). The strain in the solid phase follows the stress decay and reaches a plateau for times $t>3\tau_{B}$, yielding an estimate for the solid strain $\gamma_{s}$. As expected, the strain experienced by liquid particles is negligible.

IV.5 Shear modulus

We collect the shear stress $\sigma_{s}$ [Fig. 9(a)] and strain $\gamma_{s}$ [Fig. 9(b)] at packing fraction $\phi\simeq 0.542$ and three values $\dot{\gamma}\simeq 0.036,0.06,0.084$ of the strain rate. The radius ranges from $2.5$ to $6$ particle diameters. As expected for a fixed droplet size, increasing $\dot{\gamma}$ results in an increase of both $\sigma_{s}$ and $\gamma_{s}$. Moreover, we find that both quantities scale approximately linearly with the droplet radius, which is a direct consequence of the increasing droplet density.

We can now extract estimates for the shear modulus of small droplets via $G=\sigma_{s}/\gamma_{s}$. In Fig. 9(c), we compare $G$ with both the bulk prediction $G_{\infty}(\rho_{s}^{\infty})$ and the bulk prediction at the actual droplet density, $G_{\infty}(\rho_{s})$. We find that both values increase slightly with $R_{0}$ and follow closely the trend of $G_{\infty}(\rho_{s})$ although being $30-40\%$ smaller. For small droplets, shear moduli are found to be close to $G\simeq 40$ and thus significantly smaller than the bulk value $G_{\infty}(\rho_{s}^{\infty})\simeq 100$.

Having extracted the stress inside droplets, for one density (using droplets composed of $80$ particles at $\dot{\gamma}\simeq 0.084$) we invert $W(\sigma_{s})=\Delta F_{0}(1+a_{W}\sigma_{s}^{2})$ to get a new estimate for $a_{W}$. In addition, we also correct our previous estimate of $a_{W}^{\text{CNT}}$ using $G$ and $\Delta P$ instead of $G_{\infty}$ and $\Delta P_{\infty}$, respectively. Here, the corrected pressure difference $\Delta P\simeq 0.8\Delta P_{\infty}$ is taken from Ref. 34. In Fig. 9(d), we compare our previous estimates with the corrected response coefficients. We now find a much better agreement between the extracted $a_{W}$ and the CNT prediction $a_{W}^{\text{CNT}}$.

IV.6 Elastic work

Since for small droplets we have access to the full effective free energy $F(n;\dot{\gamma})$, we can extract the excess

[TABLE]

due to the shear flow for all droplet sizes (not only critical droplets). We split this excess into the elastic work

[TABLE]

eliminating $G$ from Eq. (2) and further contributions $\Delta W$ that are not captured by the theory presented in Sec. II. In Fig. 10, we plot $\Delta F(n)$ together with the direct evaluation of $W_{e}$ for large droplets obtained from the seeded simulations, where we employ $V_{s}=\frac{4\pi}{3}R_{0}^{3}$ approximating the radius of tension $R_{s}$ by $R_{0}$. While covering different sizes, we find that both estimates follow the same trend although $\Delta F$ clearly overestimates the true elastic work $W_{e}$. This strongly indicates that the elastic contribution does not capture the entire change in the nucleation work for sheared melts and a positive term $\Delta W>0$ is still missing.

To gain further insight, we consider the scaling of $W_{e}$ with the number $n$ of solid particles. Employing bulk quantities, the stress $\sigma_{s}$, strain $\gamma_{s}$, and the density $\rho_{s}$ would be independent of $n$, which would result in a scaling $W_{e}\sim n$ since $V_{s}\sim n/\rho_{s}$. In contrast, we observe $W_{e}\sim n^{1+1/3}$, which agrees with the linear increase of $\sigma_{s}$, $\gamma_{s}$, and $\rho_{s}$ with respect to $R_{0}$ as seen in Fig. 9 and Fig. 7(a). Interestingly, we observe a different scaling for $\Delta F\sim n^{2/3}$, suggesting a work that is dominated by the droplet area. We confirm this result by extrapolating $W_{e}$ to small sizes and plotting the difference $\Delta W=\Delta F-W_{e}\sim n^{2/3}$ [inset of Fig. 10], which indeed behaves as a surface term. Moreover, we find $\Delta W\simeq k_{B}T$ for $n\sim 40$, which can already result in a significant change in the nucleation kinetics. For large droplets, we expect a crossover to a regime in which the elastic volume term $W_{e}$ will dominate the excess $\Delta F$.

V Discussions

Performing direct and seeded simulations of sheared hard spheres, we have demonstrated that an elastic contribution plays an important role in the change of the nucleation work for sheared liquids. Including this contribution, CNT makes near-quantitative predictions as long as one takes into account that (i) the density inside droplets does not reach its bulk value, $\rho_{s}<\rho_{s}^{\infty}$, and (ii) the solid shear stress is significantly larger than the ambient liquid shear stress.

Extracting independently the work difference $\Delta F=F-F_{0}$ from free energy calculations and the elastic work $W_{e}$ from computing the local stress and strain of small droplets, we find the two estimates to be consistent with each other, although they do not agree perfectly. We trace back this discrepancy to an additional work $\Delta W$ that the system has to spend in order to form a droplet. This suggests that the overall nucleation work to form a solid nucleus of size $n$ takes the form

[TABLE]

As discussed in Ref. 13, one expects that shear induces a global deformation of the nucleus into a more ellipsoidal shape, which we have neglected here. To lowest order, one expects an increase $\delta A$ of the droplet surface that scales as $\delta A\sim\dot{\gamma}^{2}$. Such an increase could explained why $\Delta W\sim n^{2/3}>0$, although one should bear in mind that it could also stem from an increase of the interfacial tension $\Gamma$. In fact, experimentally it has been observed for a colloidal gas-liquid interface that shear suppresses capillary waves and tends to increase $\Gamma$ Derks et al. (2006). Here, the microscopic picture is a local erosion at the interface, where particles rattle due to the stress and eventually dissolve into the liquid.

Here we have focussed on the pressure difference $\Delta P$ as the natural driving force of crystallization. For an incompressible solid, one can rewrite the driving force as the (absolute) difference of chemical potential $|\Delta\mu|=\mu_{l}-\mu_{s}$ between the solid and the liquid phase at the same ambient pressure $P_{l}$. For a quiescent liquid, increasing $P_{l}$ beyond the coexistence pressure $P_{\text{coex}}$, the two branches $\mu_{l}(P)$ and $\mu_{s}(P)$ move apart from each other increasing $|\Delta\mu|$ and thus leading to a faster nucleation process. Hence, to explain the suppression of nucleation under shear using the same framework one would require a non-equilibrium liquid chemical potential that becomes smaller. As discussed by Butler and Harrowell Butler and Harrowell (2002, 2003), there is no consistent definition of a non-equilibrium chemical potential that would predict such a reduction.

Indeed, we provide new numerical evidence in Appendix B that confirms the opposite scenario with an increase of $\mu_{l}$ as a function of $\dot{\gamma}$. To do so, we have employed the fast growth method Hendrix and Jarzynski (2001) to extract the insertion work $w_{\text{ex}}$ needed to place a particle in a (dense) melt. For a quiescent liquid this work reduces to the excess chemical potential $\mu_{l}^{ex}=\mu_{l}-\mu_{l}^{id}$, with $\mu_{l}^{id}$ the chemical potential of an ideal gas. Consistent with the flow symmetry, we find in the linear regime that $w_{\text{ex}}$ increases with $\dot{\gamma}^{2}$. Identifying this work with a non-equilibrium liquid chemical potential, it moves away from the solid branch and increases the driving force $|\Delta\mu|$, which is inconsistent with the observed suppression of nucleation.

VI Conclusions

In this paper, we have tested a possible extension of classical nucleation theory to model the crystallization of sheared liquids. We have demonstrated that one of the key ingredient to model the change of nucleation work as a function of the imposed strain rate $\dot{\gamma}$ is to take into account the elastic deformation of the newborn nucleus. Such an extension can already predict, on a qualitative level, many observations made in previous simulations, namely the drop of the nucleation rate and the increase of both the nucleation work and the critical nucleus size as a function of $\dot{\gamma}$.

To obtain quantitative predictions, one needs to evaluate the elastic work $W_{e}$, which contains the shear stress $\sigma_{s}$ inside the droplet and its shear modulus $G$. These two quantities are unknown, hence one has to employ approximations. As a first step, one might try to: (i) equating the solid shear stress with the ambient liquid stress $\sigma_{l}=\eta\dot{\gamma}$ and (ii) to treat small droplets as bulk phases and employ the bulk shear modulus $G_{\infty}$ of a macroscopic solid. In Sec. IV.2, we have demonstrated that these approximations do not lead to satisfactory results.

In the second part of this paper, we have employed a seeding method that generates configurations with larger droplets under shear, from which we have extracted the local solid density $\rho_{s}$, the shear stress $\sigma_{s}$, and the strain $\gamma_{s}$. Our first finding is a direct consequence of what is seen in quiescent crystallization, namely that the density at the center of droplets does not reach its bulk value, leading to an effective smaller shear modulus compared with the bulk value. Our second observation is the sharp increase of shear stress inside the nucleus compared with the ambient hydrodynamic stress $\sigma_{l}$. Employing the actual values $\sigma_{s}$ and $G$, we indeed find a much better agreement between numerics and CNT.

Finally, we have highlighted the presence of another term entering the nucleation work, $\Delta W$, and have argued that it arises from a positive contribution to the interfacial work. Due to scaling reasons, this term dominates the work for small nuclei but should become negligible compared to the bulk elastic work for large droplets. How the interfacial tension $\Gamma$ changes when the liquid layer at contact with the droplet is sheared remains to be investigated.

In the present manuscript, we have not addressed the regime of deeply supercooled liquids, in which the change of the nucleation work competes with the enhanced dynamics that increases the kinetic prefactor present in the nucleation rate. A natural extension of this work includes quantifying how the diffusion in droplet size space varies as a function of the liquid density and strain rates.

Acknowledgements.

We gratefully acknowledge ZDV Mainz for computing time on the MOGON supercomputer. We acknowledge financial support by the DFG through the collaborative research center TRR 146.

Appendix A Brownian dynamics vs. molecular dynamics

We perform molecular dynamics (MD) simulations in the NVT ensemble using the Lowe-Anderson thermostat Koopman and Lowe (2006) has in our previous work, Ref. 31. In these simulations, time is measured in units of $\sqrt{m\alpha^{2}/k_{B}T}$, with $m$ being the mass of particles. A velocity Verlet integrator is used with a time step $\Delta t=0.004$. The coupling between the system and the thermostat is regulated by the bath collision frequency $\Gamma=10$. Rates are again extracted via the plateau of the mean first passage time. We then rescale the MD nucleation time in units of the Brownian time $\tau_{B}^{MD}=d^{2}/D_{0}$, where

[TABLE]

is the bare self-diffusion coefficient taken from the Chapman-Enskog kinetic theory of gases Chapman et al. (1990). It allows a rescaling of the nucleation rate $k$ as well as the strain rate $\dot{\gamma}$. In Fig. 11, we compare nucleation rates extracted from BD and MD at various packing fractions. We find a good agreement between the two dynamics indicating the consistency between this work and our previous study in Ref. 31.

Appendix B Non-equilibrium chemical potential

The fast growth method Hendrix and Jarzynski (2001) is based on the Jarzynski equality

[TABLE]

relating the free energy difference $\Delta\mathcal{F}$ between two states $A$ and $B$ to the distribution of work $w_{\tau}$ performed on the system to move it from $A$ to $B$ within time $\tau$. Here, we evaluate the free energy difference between a liquid composed of $N$ and $N+1$ particles at a fixed volume $V$, which is nothing than the microscopic definition of the chemical potential $\mu_{l}$. The work $w_{\tau}$ is computed from a discrete protocol, where we progressively switch on the interaction between a tagged particle and the surrounding fluid through changing the parameter $\lambda$ from [math] to $1$. The work reads

[TABLE]

where $H_{\lambda_{t}}(\omega_{t})$ is the Hamiltonian at time $t$ with microstate $\omega_{t}$ (position and velocities of all particles). In practice, we use the protocol $\lambda(t)=(t/\tau)^{6}$. In the limit of instantaneous switching, $\tau\to 0$, one recovers the Widom insertion method Allen and Tildesley (2017). In Fig. 12(a), we show the (excess) liquid chemical potential $\mu_{l}-\mu_{l}^{\text{id}}$ as a function of the density $\rho_{l}$ for $\tau=10^{5}\Delta t$ and we compare these new results with the thermodynamic integration used in Ref. 40. We find a perfect agreement with our previous parametrization even beyond the freezing density $\rho_{f}$. In contrast, the Widom insertion method ($\tau\to 0$) is only able to compute chemical potential differences that are below $\sim 10k_{b}T$. Applying the fast growth method in a sheared liquid, we find a continuous increase of the insertion work $w_{\text{ex}}=-\ln(\langle\exp(-\beta w_{\tau})\rangle)$ as a function of the strain rate $\dot{\gamma}$, cf. Fig. 12(b).