# Flow induced crystallization of penetrable particles

###### Abstract

For a system of Brownian particles interacting via a soft exponential potential we investigate the interaction between equilibrium crystallization and spatially varying shear flow. For thermodynamic state points within the liquid part of the phase diagram, but close to the crystallization phase boundary, we observe that imposing a Poiseuille flow can induce nonequilibrium crystalline ordering in regions of low shear gradient. The physical mechanism responsible for this phenomenon is shear–induced particle migration, which causes particles to drift preferentially towards the center of the flow channel, thus increasing the local density in the channel center. The method employed is classical dynamical density functional theory.

## I Introduction

The creation of a colloidal crystal generally proceeds via the slow process of nucleation and growth within a dense, crowded environment oxtoby (). Technological applications requiring periodic order on long length scales (e.g. photonic crystals) are hindered by the defects and grain boundaries which inevitably develop during the growth process. It is thus desirable to identify ways by which external fields can be employed to control the nucleation dynamics and thus optimize the quality of the resulting crystal. One way to achieve this aim is to drive the system out of equilibrium using mechanical deformation of the sample.

Since the pioneering light scattering investigations of Ackerson and Pusey pusey () it has been known that the judicious application of an oscillatory shear strain can induce three-dimensional crystalline order in colloidal fluids. In a similar spirit, Besseling et al. have performed real-space confocal microscopy experiments on colloidal mixtures under oscillatory shear besseling (). Crucially, in both cases the bulk density of the experimental sample lies below that of the equilibrium freezing transition, such that these shear–induced crystals represent true out-of-equilibrium states; the microstructure relaxes back to equilibrium following the cessation of the flow. Although Besseling et al. besseling () supported their experimental findings with simulation data, there currently exists no first-principles theoretical approach to address this phenomenon. In a recent study, Peng et al. have used molecular dynamics to study the influence of steady simple shear on the crystal growth kinetics voigtmann (). However, in contrast to Refs. pusey (); besseling (), the thermodynamic statepoints considered in Ref. voigtmann () were all in the crystal region of the phase diagram.

The aforementioned studies of flow–induced crystallization pusey (); besseling () considered simple shear flow, where the (time-dependent) shear gradient is constant in space. However, there are many situations of interest for which the shear gradient is a function of position. For example, in the commonly encountered case of Poiseuille flow along a channel the shear gradient is a linear function of the distance from the channel center. It is well–known that when the shear–rate varies significantly on the scale of a particle diameter, then the particles will begin to exhibit a biased diffusion towards regions of low gradient: shear–induced migration LA (). The physical origin of this effect is that the collision frequency of a given particle with it’s neighbours is not isotropically distributed over the surface of the particle; surface regions subject to a higher shear–rate will experience, on average, a greater number of collisions than those regions subject to a lower shear-rate MSB (). In the case of Poiseuille flow this mechanism leads to an increase in density at the centre of the channel.

In this paper we will explore how the particle migration induced by Pouseuille flow can interact with the underlying equilibrium free energy to generate nonequilibrium ordered states in which a crystalline region appears at the center of the channel. In addition to presenting a novel kind of nonequilibrium state arising from the coupling of shear–rate gradients to free energy minima, we can imagine that these nonequilibrium steady states may be relevant for particle transport in certain microfluidic devices. The method we employ is dynamical density functional theory, modified along the lines of Refs. KB1 (); KB2 (); SKB () to incorporate the non-affine particle motion necessary to capture particle migration effects.

The paper will be organized as follows: in section II we explain the model, in section III we give an overview of the employed theory, in section IV we report our results and finally, in section V, we provide comments and an outlook for future work.

## Ii Model

In order to investigate the process of shear–induced crystallization we choose to employ a simple, minimal model of penetrable particles, namely the generalized exponential model (GEM). Within this model the purely repulsive pair interaction potential has an exponential form, with an exponent which can be adjusted to tune the strength of the repulsive force between particles. The continuous nature of the potential makes it particularly well suited for numerical studies. We consider the GEM model with an exponent of , described by the pair-potential

(1) |

where and set the length and energy scales, respectively.

Ordering phenomena in GEM models, especially concerning quasi-crystalline order andy_quasi_crys_1 (); andy_quasi_crys_2 (), solidification andy_solidification (), generation of defects and disorder from quenching andy_quenching (), crystallisation under confinement at interfaces and in wedges Archer_GEM () have been recently addressed. An important feature of penetrable particle models, such as the GEM-8, is that at high densities they exhibit a so–called ‘cluster crystal’ phase. Cluster crystals differ from more familiar hard–sphere type crystals in that each density peak can contain, on average, more than one particle. While the multiple occupancy of a given lattice site costs a potential energy of order per particle pair, the resulting gain in entropy dominates to minimize the grand potential. This mechanism makes the study of crystallization in penetrable particle models more subtle than the case of simple hard spheres, for which only entropy plays a role. However, penetrabel particle systems have the significant benefit that simple mean field approaches can yield accurate results.

We consider a two-dimensional GEM-8 fluid system confined between a pair of parallel repulsive walls. We choose to work in two dimensions as this is the minimal situation in which crystallization can be studied and reduces the numerical demands of integrating the DDFT equations - three dimensional calculations would be very expensive. Moreover, the majority of the literature concerning the equilibrium GEM-8 model addresses the two-dimensional system.

The suspension is confined to a channel by the following external potential field

(2) |

where is the distance between the walls. In a dynamical calculation the infinitely repulsive walls described by (2) are equivalent to imposing a no-flux boundary condition. We choose the shear deformation such that the velocity field takes a Poiseuille form, with stick boundary conditions. The velocity field is given by

(3) |

where is the maximal velocity (which occurs in the center of the channel).

## Iii Theory

The dynamical density functional theory (DDFT) provides a convenient approximate method to study the dynamical response of the particle number density, , to external force fields. The original formulation of DDFT, which can be derived from either the Langevin marconi_tarazona_1 (); marconi_tarazona_2 () or Smoluchowski (archer_evans, ) descriptions, considered time–dependent external potential fields. Generalization to treat external flow fields (e.g. shear) was first made by Rauscher et al. by incorporating an additional term representing the affine solvent flow field rauscher (). However, it subsequently became apparent that this simple approach ignores important non–affine particle motion, which renders the approach incapable of describing interaction–induced currents orthogonal to the affine flow. As a consequence, the DDFT of Ref. rauscher () cannot capture the physics of laning, particle migration or other phenomena arising from a coupling between flow and interparticle interactions.

These shortcomings were addressed by Krüger, Brader and Scacchi in Refs. KB1 (); KB2 (); SKB (), who reintroduced the missing non-affine motion into DDFT using a dynamical mean–field approximation. In it’s most recent form the DDFT equation is given by

(4) |

where is the solvent velocity field, is the mobility and is the Helmholtz free energy containing details of the interparticle interactions and the external potential field.

The Helmholtz free energy can be split into ideal, excess and external field contributions

(5) |

The ideal gas part is known exactly and is given by

(6) |

where we have set the thermal wavelength equal to unity. For the GEM-8 model presently under consideration the excess part can be well described by the mean–field functional

(7) |

The central approximation made here is that the two body density can be written in factorized form, i.e. . Despite the simplicity of this approximation, it has recently been pointed out that it performs better than one would expect chacko ().

The velocity field can be decomposed into an affine term and a fluctuation term, according to

(8) |

Krüger and Brader have shown KB1 (); KB2 () that the fluctuation term in the velocity field can be approximated by the mean–field form

(9) |

where is a time-independent Kernel depending on the interparticle interaction potential and is the position dependent shear flow. While the convolution form of the mean–field term is both physically plausible and computationally appealing, it nevertheless represents an empirical add-on to the microscopically derived original versions of DDFT marconi_tarazona_1 (); marconi_tarazona_2 (); archer_evans () - a true microscopic derivation would indeed be very useful. As a consequence the detailed form of the convolution kernel remains unspecified and lacks a rigorous microscopic prescription. One can, however, give a physical interpretation to the kernel in terms of binary collision events, which facilitates the development of approximations.

The product of the kernel with the shear–rate in (9) has the units of velocity. This non-affine velocity arises from the force of interaction between a pair of particles when they are driven into each other by the flow field (‘flow–interaction coupling’) and thus encodes a non–affine motion. Physically, one can envisage ‘rolling over’ type motion as two spherically symmetric particles attempt to follow the affine flow as closely as possible. Noting that for over damped systems the velocity and force are equivalent (up to a friction coefficient), a simple ansatz for the kernel is

(10) |

In general one would expect the coefficient to have a density dependence (not necessarily local). However, in the absence of detailed information about this function we choose the simplest possibility and set . We will see that this already leads to some interesting phenomenology and leaves open the possibility for fine tuning in order to fit data from numerical simulations or experiments.

The equilibrium state corresponds to the long-time limit of equation (4) in the absence of flow. Equivalently, the equilibrium density can be obtained by minimizing the grand potential

(11) |

where is the chemical potential. The minimization generates the Euler-Lagrange equation

(12) |

where and the one-body direct correlation function is given by a convolution with the interaction potential

(13) |

The bulk phase diagram shown in Fig. 1 has been calculated by solving Eq. (12) in the case of vanishing external potential.

For densities above the spinodal the uniform liquid becomes linearly unstable with respect to the inhomogeneous crystal–phase density distribution. As previously mentioned, because of the soft–core nature of the particles the crystal phase consists of clusters of particles occupying the sites of a two-dimensional hexagonal lattice. The phase diagram shown in Fig. 1 applies to an infinite system and we note that small discrepancies can occur in our numerical DDFT calculations arising from finite size effects. As the crystallising system is rather sensitive to the exact location of the statepoint we have been careful to keep to a minimum finite size effects arising from the confinement/periodic boundaries employed in our dynamical calculations. For more details about boundary effects, we refer the reader to Archer_GEM (), where a detailed study of crystallisation under confinement at interfaces and in wedges is performed.

## Iv Results

We prepare the system at a thermodynamically stable state point close to the crystal phase, as shown in Fig. 1. In Fig. 2 we show the equilibrium (unsheared) density distribution of a system at this statepoint; average density , and with size . The average density is calculated by integrating over the density profile and dividing by the system volume. Note that to improve visualization of the density profile we saturate the colorbar at a density value of . This initial density distribution was calculated by solving equation (12) at chemical potential . As we are working quite close to the freezing phase boundary there is clear evidence of pre-crystallisation at each wall. We choose to perform our calculations quite close to the phase boundary for two reasons: (i) it is interesting to see how switching on the Poiseuille flow will influence the pre-crystal regions at each wall and (ii) the flow–induced crystallisation phenomena we wish to investigate can be observed at relatively low shear rates; high shear rates generate stability issues which make difficult an accurate numerical solution of the DDFT equation.

The integration step used here is . To avoid numerical issues, we switch on the flow using a ramp function for the first Brownian time units. After this time, the flow reaches is maximum value with . At this moment, the system is only slightly distorted, with the majority of the distortion around the walls. The density peaks become smeared out and a stronger laning effect is visible over the entire system, see Fig. 3.

From onwards the system undergoes a steady Poiseuille flow. For about further time units the density does not undergo any significant change. At (see Fig. 4) the density distribution starts to change considerably as the migration mechanism pushes more particles to the center of the channel and the system undergoes a symmetry breaking along the flow direction in the central region of the channel. The ‘stripes’ which emerge have a parabolic form similar to that of the external flow, but can only be observed within a relatively short time-window. These structures serve as a precursor to the formation of a more crystalline region at the center of the channel. As time progresses each of the stripes develops some internal structure and distinct density peaks begin to emerge: Fig. 5 shows the situation at .

Within approximately 10 time units the density peaks stabilize to form a steady–state and no further qualitative change in the density is observed. In Fig. 6 we show a representative steady–state density distribution at . Because the system is under constant external drive the steady–state crystalline structure at the channel center is subject to continuous deformation; the density peaks try to follow the affine flow, but have to ‘squeeze past’ their neighbours in order to do so. Despite this constant rearrangement a general hexagonal packing structure can be identified at any given time.

To provide an alternative visualisation of this phenomenon we show in Fig. 7 one-dimensional density distributions obtained by integrating the two-dimensional density along the direction of flow. Data is shown for five different times to illustrate important stages of the time-evolution.

Following the onset of flow the density gradually increases in the center of the system (from brown to yellow lines) as a consequence of shear-induced particle migration. Once the local density around the channel center approaches the bulk crystallization phase boundary density peaks start to form in the two-dimensional density distribution, which lead to the strong oscillations in the integrated data shown in Fig. 7 for (c.f. Fig. 5). For later times we enter a steady–state situation and the oscillations get smoothed out as a result of averaging over the continuous distortion of the crystal structure when in steady–state.

Within the framework of DDFT the flow-induced structural changes discussed above can be related to changes in the free energy of the system. The fact that we can analyze dynamic phenomena using the equilibrium concept of a free energy is a consequence of the adiabatic approximation underlying the DDFT.

In Fig. 8 we show the time-evolution of the Helmholtz free energy. For the driven system under consideration we recall that there is no ‘H-theorem’ and that the free energy is not required to decrease on approach to the steady state; the applied shear flow is constantly adding energy to the system. The points in the figure indicate times for which the two-dimensional data is shown in the previous figures. The free energy increases steadily up to around and then remains constant up to around . Up until this point in time the density does not show significant change, other than the erosion of pre-crystal structures at each wall. At around the first indication of a symmetry breaking occurs as the ‘stripes’ start to appear in the center of the channel and this is reflected by a sudden increase in the free energy. At the stripes become unstable and density peaks develop, leading to another jump in free energy. For later times the local crystal is formed and the value of the free energy oscillates. These oscillations are due to the continual deformation of the crystal as the nonuniform flow pushes the central peaks at a larger velocity than those located to the left or right of the center. We note that the details of these oscillations are somewhat influenced by the finite system size employed in our calculations.

## V Discussion

The DDFT is a well established method to calculate approximately the relaxation to equilibrium of the one-body density in overdamped systems. Its direct application to driven systems is problematic due to the neglegt of nonaffine particle motion; a phenomenological correction to the theory, as embodied by the flow kernel, must be introduced to account for these in a mean-field fashion. We have presented a new approximation for the nonaffine velocity field, which is appropriate for treating systems with soft, penetrable interparticle interactions.

By applying our theory to the GEM-8 model under Poisseuile flow at statepoints close to the freezing transition we have investigated the interaction between crystallisation and shear-induced particle migration. The latter mechanism is not accounted for in standard DDFT and only enters as a result of our mean-field treatment of the nonaffine velocity. Following the onset of flow we observe that, after a period of transient dynamics, a local crystalline steady-state can form at the center of the channel. We anticipate that this phenomenon, which we have not seen reported elsewhere, will be generic for a wide class of soft particle models. Whether hard-spheres or similiar systems with strongly repulsive interactions would also exhibit this effect remains an open question.

The present study presents much opportunity for development and future investigation. Of considerable interest would be stochastic simulations to establish the limitations of our phenomenological DDFT results. It would also be very interesting to investigate the impact of microstructural ordering (laning and local crystallization) on the rheology of the system. For example, one could enquire whether the transition to a nonequilibrium crystalline state increases the throughput in channel flow at a prescribed pressure. These and other questions will be the subject of further research.

On a more fundamental level, we are currently employing methods of bifurcation/stability analysis to investigate in more analytic detail the onset of symmetry breaking in ‘non-affine DDFT’, i.e. DDFT with a convolutional flow kernel term, such as that considered in the present work. We have so far focussed on the laning instability under simple shear flow instability (), but we anticipate that similar techniques could be fruitfully employed to study in detail, and with reduced numerical effort, more general flow induced crystallization and ordering phenomena. Work in this direction is currently in progress.

## Vi Acknowledgement

We thank Andrew J. Archer for interesting discussions and the Swiss National Science Foundation for financial support under the grant number 200021-153657/2.

## References

- (1) D. W. Oxtoby, in fundamentals of inhomogeneous fluids, ed. D. Henderson, p. 407 (1992)
- (2) B. J. Ackerson and P. Pusey, Phys.Rev.Lett. 61 1033 (1988)
- (3) T. H. Besseling et al., Soft Matter 8 2931 (2012)
- (4) H. L. Peng, D. M. Herlach and Th. Voigtmann, arXiv:1707.0551v1 (2017)
- (5) D. Leighton and A. Acrivos. J.Fluid.Mech., 181 415 (1987)
- (6) P. Marmet, A. Scacchi, J.M. Brader, Mol.Phys. 115 1691 (2017)
- (7) A. Scacchi, M. Krüger, J.M. Brader, J. Phys.: Condens. Matter 28 244023 (2016)
- (8) J. M. Brader and M. Krüger, Mol.Phys. 109 1029 (2011)
- (9) M. Krüger and J. M. Brader, EPL 96 68006 (2011)
- (10) A. J. Archer, A.M. Rucklidge, E. Knobloch, PRL 111 165501 (2013)
- (11) A. J. Archer, A. M. Rucklidge, and E. Knobloch, PRE 92 012324 (2015)
- (12) A.J. Archer, M.C. Walters, U. Thiele and E. Knobloch, PRE 90 042404 (2014)
- (13) A.J. Archer, M.C. Walters, U. Thiele and E. Knobloch, Springer Proceedings in Mathematics and Statistics 166 1 (2016)
- (14) A. J. Archer and A. Malijevský, J. Phys.: Condens. Matter 28 244017 (2016)
- (15) U. M. B. Marconi and P. Tarazona, J. Chem. Phys. 110 8032 (1999).
- (16) U. M. B. Marconi and P. Tarazona, J. Phys.: Condens. Matter 12 A413 (2000).
- (17) A.J. Archer and R. Evans, J. Chem. Phys. 121 4246 (2004).
- (18) M. Rauscher, A. Dominguez, M. Krüger and F. Penna, J. Chem. Phys. 127 244906 (2007)
- (19) A.J. Archer, B. Chacko and R. Evans, J. Chem. Phys. 147, 034501 (2017).
- (20) A. Scacchi, A.J. Archer and J.M. Brader, arXiv:1711.00662 (2017).