Confinement-induced demixing and crystallization

We simulate a strongly size-disperse hard-sphere fluid confined between two parallel, hard walls. We find that confinement induces crystallization into n-layered hexagonal lattices and a novel honeycomb-shaped structure, facilitated by fractionation. The onset of freezing prevents the formation of a stable glass phase and occurs at much smaller packing fraction than in bulk. Varying the wall separation triggers solid-to-solid transitions and a systematic change of the size-distribution of crystalline particles, which we rationalize using a semi-quantitative theory. We show that the crystallization can be exploited in a wedge geometry to demix particles of different sizes.

An important aspect has, however, not yet received much attention, which is the effect of particle sizedispersity on the above phenomenology. In laboratory experiments on colloids a size-dispersity occurs naturally from synthesis. Additionally, for the study of structural relaxation in supercooled liquids and glasses a sizedispersity must be introduced to prevent crystallization even at very small packing fractions [6,12,13]. Despite their popularity as model glass formers recent observations in bulk have revealed that Gaussian-distributed hard spheres [36][37][38][39] as well as the often used Kob-Anderson model [40,41] form crystals already in the supercooled regime. This crystallization is induced by a process called fractionation which describes the separation of a homogeneous fluid into different liquid or crystalline fractions with very different particle-size distribu- * Corresponding author: gerhard.jung@uibk.ac.at tions. Since walls or other boundaries lead to heterogeneous crystallization which is known to strongly effect and often enhance nucleation processes [42,43] it should be expected that confinement could also have an impact on crystallization and fractionation. Understanding the equilibrium structural properties of confined size-disperse particles thus remains a critical outstanding problem and is essential for future studies of simple glass formers and colloids in confinement.
In this work we study the crystallization of sizedisperse hard spheres in a slit geometry. We use two complementary approaches: First, we perform exact eventdriven molecular dynamics (EDMD) simulations with enhanced sampling techniques to reveal a complex phase behaviour, displaying amongst others a new phase which is shaped like a honeycomb lattice. Significantly, we find formation of stable crystals at much smaller packing fractions than reported for bulk systems, and show that the crystallization is enabled by a confinement-controlled fractionation. These surprising results show additionally that the multiple-reentrant glass transition reported in Ref. [6] is only metastable. Second, we rationalize our findings with a semi-quantitative theory allowing us to provide a deep insight into the mechanisms that drive the confinement-induced crystallization. Most importantly, our study reveals a very general technique for the demixing of size-disperse particles with excluded volume interactions.

II. METHODS
The investigated system consists of two parallel, hard walls separated by a distance H in the z-direction and periodic boundary conditions in the other two dimensions. The slab is filled with size-disperse hard spheres at packing fraction φ. In the spirit of recent simulations in bulk [39] we employ a Gaussian distribution with average particle diameter σ and dispersity δ = 0.15. Here, σ sets the length scale of the system and the dispersity δ is defined as the ratio of the standard deviation of the particle diameter to the average particle diameter.
A. Event-driven molecular dynamics simulations The system is simulated using a hybrid scheme of event-driven molecular dynamics simulations (EDMD) in the NVT ensemble [44][45][46] and the recently proposed swap Monte Carlo algorithm (SWAP) [47][48][49]. SWAP introduces "swap events" that are triggered every ∆t = 0.5 σv −1 th time steps. Each swap event consists of N s = 50, 000 Monte Carlo moves in which two particles are randomly selected and their diameters are exchanged. The move is accepted if the particles exhibit no overlap with other particles after the swap. This corresponds to the typical Monte Carlo acceptance criterium for hard spheres and is consistent with detailed balance. The hybrid scheme is implemented in the EDMD code DynamO [46]. Its usage is essential to accelerate the dynamics and systematically ensure that the equilibrium state of the system is attained. We observe typical acceptance ratios of Γ = 0.05, showing that a significant proportion of the particles is swapped in each swap event. The slab consists of N = 13068 particles for any wall separation 2.0 ≤ H/σ ≤ 3.0 resulting in a minimum longitudinal box size of L x = L y = 67.15 σ, which is large enough to eliminate any finite size effects. The system is equilibrated for N e = 3 · 10 10 events, which roughly corresponds to a total equilibration time of t eq = 7.0 · 10 4 σv −1 th . To determine the size-distribution of particles in the solid phase we define the rotationally invariant bond-order parameter d 6 (m, n) = φ 6,m φ 6,n . Here, we introduced the local bond-orientational parameter, φ 6,m = N −1 m Nm m =1 e 6iθ mm , as the relative orientation of particle m with its N m neighbors [50][51][52]. Neighbors are all particles m within a cut-off distance r < 1.4 σ and r ⊥ < 0.1 σ, with ⊥ denoting the direction perpendicular to the walls. Using these invariant measures we define a particle m as crystalline if it has at least four neighbours n for which the rotationally invariant bond-order parameter d 6 (m, n) > 0.7. Above this threshold, the local environment of particle m is strongly correlated with the one of neighbor n which indicates a solid-like structure. To identify the big particles in the honeycomb phase 3 we set r < 1.8 σ and require only three crystalline neighbours. We report the packing fraction φ c at which the onset of crystallization is observed, as well as the average particle diameterā and dispersity δ c of particles in the solid phase [50][51][52].

Fundamental measure theory (FMT)
FMT is based on density functional theory and minimizes the functional for the grand free energy leading to the equation [53,54], Here, n i (z) is the local number density, µ i the chemical potential [55] and V i (z) the external wall potential of component i. Each component has a different diameter a i = a min + (i + 0.5)∆a, which allows us to emulate a size-dispersity by creating a mixture of m different components, similar to Ref. [6]. We use m = 25 with ∆a = 0.024 σ and a min = 0.7 σ. We also ensured that the density profiles and the resulting free energies only marginally change (∆F < 0.02 k B T ) upon an increase of m. For the excess free energy functional of a hard-sphere mixture F ex we use the White Bear Version II [54]. To determine the density profile n i (z), Eq. (1) is solved selfconsistently in an iterative procedure [54]. Different from Ref. [6] we ensure in our algorithm that the chemical potential µ i is adapted in each iteration such that the resulting slit-averaged concentrations of every component n i (a i ) = dz n i (z) precisely match a Gaussian profile with mean σ and variance δ. The free energy per particle of component i can then be evaluated from its ideal and excess contributions, The constant λ i is set to λ 3 i = σ 4 ∆a −1 , such that ρ i (z) = λ 3 i n i (z) = σ 4 ∆a −1 n i (z) is the density distribution of the disperse mixture in the spirit of Ref. [36]. We carefully checked this expression and ensured that it is independent of a min and ∆a and indeed gives the correct ideal gas contributions (e.g. the expected values for the mixing free energy).
The results for the free energy per particle of a hardsphere mixture in the liquid state are shown in Fig. 1. It can be observed that the maximum of the distribution is shifted to around 1.07 σ. This is because the accessible width of the slab for hard spheres is L = H − a, which means that the density of big particles (with large a and thus small L) in the accessible region is higher than the one of the of small particles. This similarly holds for the ideal free energy in the liquid phase F fl (a). Furthermore there is a pronounced non-monotonic dependence of the free energy on the wall separation H. In Ref. [4,6] a similar non-monotonic behavior in the structure factor was explained with commensurate and incommensurate packing. Commensurate packing describes the favorable packing of spheres in k layers for wall separations H/σ ≈ k, k ∈ N, while incommensurate packing relates to wall separations close to half-integer values. The latter is less favorable since single particles have to be included between the k layers which also leads to a slowing-down of diffusion.

Cell theory
To determine the free energy of particles in the crystalline phase we employ cell theory [8,[56][57][58]. As input we chose five different lattice types (2 , 3 , 3 as well as the square lattices 2 and 3 ) which can be described by two or three lattice parameters c i . One of these parameters is fixed to set the desired packing fraction φ. For each lattice type C we create a minimal structure of the crystals in a slab geometry with wall separation H and packing fraction φ and place particles on the lattice sites with diameterā. Only one of the lattice sites is occupied with a test particle of diameter a. For this particle we determine the free volume cell using a Monte Carlo procedure and thus determine its free energy based on the assumptions used in cell theory. For each H, φ,ā and a we minimize this free energy with respect to the free lattice parameters c i . In crystals with three layers, there is a difference in the free volume cell of particles in the centre of the lattice or at the wall. In this case we determine the contributions separately and minimize the averaged free energy per particle.

Determination of the stability diagram
Having determined the free energies in the fluid and the crystalline phase, we assume a top-hat distribution f (a |ā, δ c ) of particles in the crystalline phase with size-dispersity δ c and create a list of free energy differences, Here we have to account explicitly for the mixing free energy ∆F mix (δ c ) = − ln( √ 12δ c ) of the particles in the crystal [59,60] and include a shift parameter ∆F shift since cell theory only provides an upper bound for the free energy. The latter was set to ∆F shift = −1.80 in Ref. [8] to match the value for monodisperse discs in the 2D limit. We adopted this value in our work. For every wall separation H we report the minimal packing fraction φ c for which the liquid and crystalline free energies are equal, ∆F = 0. Additionally, we determine the corresponding crystalline type C, the mean particle diameter a and dispersity δ c . In a monodisperse system a coexistence region of liquid and crystalline phases can be determined from the free energies using the common tangent construction. The underlying assumption is that inside the coexistence region both the liquid and the crystalline phase have a constant packing fraction (and obviously a constant particlesize distribution). The fractionation in size-disperse systems, however, implies that the particle-size distribution in the liquid phase changes with the onset of crystallization. In Refs. [36,61] is was shown that these systems can be handled using the moment free energy method, leading to highly coupled non-linear equations. In the case of confinement the application of this method is unfeasible since the free energies can only be determined numerically and it is not clear how they could be rewritten such that they only depend on four moments of the (position-dependent) density distribution. Here, we thus rely on a semi-quantitative scheme and identify the "critical" packing fraction φ c in a system with wall separation H as the packing fraction at which the free energies in the crystalline and the liquid phase are equal.

III. STABILITY DIAGRAM IN SLAB GEOMETRIES
We first study the crystallization of hard spheres confined in slabs with different wall separations 2.0 < H/σ < 3.0.. The system displays an onset of crystallization at packing fractions around φ ≈ 0.5, much smaller than the bulk value φ b ≈ 0.6 [39]. Four different crystalline phases can be observed, namely an amorphous liquid phase for small packing fraction φ, a two-layered hexagonal structure (2 ) for small H, a three-layered hexagonal structure (3 ) for large H and an intermediate, honeycomb-shaped phase (3 ). A visualization of the different phases is shown in Fig. 2, the stability diagram is displayed in Fig. 3. The sequence of solid-to-solid transitions 2 → 3 → 3 upon a change of wall separation is reminiscent of similar transitions in the case of monodisperse hard spheres or charged colloids [8,62,63]. Differ- Color code indicates particle diameter with a < 0.9 σ (blue), a > 1.1 σ (red) and interpolation between. The sphere diameter is scaled by a factor of 0.5 and the lower layers are partially transparent for the sake of visualization. ent from these monodisperse systems, the size-dispersity leads to a significantly enlarged domain where the hexagonal structures are stable. In the small intermediate domain between 2 and 3 we observe the honeycomb phase, which consist of three layers of small particles forming honeycomb cells and big particles in the upper and lower layer which fill these cells, visualized in Fig. 2 (3 ). The formation of this honeycomb structure leads to the suppression of a square-lattice phase which was observed in the case of small or vanishing dispersity [8,12]. In fact, for δ = 0.15 the square-lattice is only stable in simulations at very high packing fraction φ > 0.53 in coexistence with a liquid and a honeycomb phase.
Comparing the critical packing fractions presented in Fig. 3 to the glass transition line determined in Ref. [6] shows that the glass is metastable and can only coexist with a crystalline phase. Interestingly, both the glass transition line and the critical packing fraction φ c exhibit Dashed line corresponds to the glass transition line of the same system as determined in Ref. [6]. Colors indicate transitions to different crystalline phases: Two hexagonal layers (2 , red), honeycomb phase (3 , green) and three hexagonal layers (3 , pink).
a pronounced non-monotonicity leading to the possibility of observing reentrant phase transitions by changing the wall separation [6,64].
The size distribution of particles in the emergent crystals at packing fraction φ φ c is strongly correlated with the wall separation H, as shown in Fig. 4A. The average particle diameterā for crystalline particles exhibits a linear dependence on H for one crystalline phase and a distinct jump upon a solid-to-solid transition. It can also be observed that the honeycomb phase is indeed interpolating between the two hexagonal structures since the diameters of the big and small particles in the honeycomb phase perfectly follow the linear dependence of the 2 and 3 hexagonal phases respectively. The size-dispersity of the crystalline particles is around δ c ≈ 0.05 − 0.06, much smaller than the dispersity of the host fluid δ = 0.15, which indicates a significant fractionation (see Fig. 4B). The dispersity of crystalline particles in confinement is very close to the one for crystals in bulk δ c 0.07, and similar for all of the crystalline phases we observed.
A major impact of size-dispersity on the stability diagram is the enlargement of domains with stable hexagonal structures compared to the monodiperse case. This is facilitated by a linear increase of the average diameterā with wall separation H, seen in Fig. 4A, leading to a constant ratio H/ā. Consequently, in the extreme case of a flat size-distribution of the host fluid, H would reduce to a mere length scale setting the average particle sizeā, leading to the existence of a single crystalline phase and a constant critical packing fraction φ c . In reality, however, this picture is oversimplified since the particle-size distribution of the host fluid already defines the length scale σ. Therefore, the values ofā are bound to σ(1 − δ) ā σ(1 + δ). Upon a further increase of H (and thusā) a discontinuous solid-to-solid transition is triggered by increasing the number of crystalline layers and reducing the average diameter of crystalline particles.
We additionally expect that the maximum in the particle-size distribution of the host fluid would lead to the existence of a minimum in φ c which implies a reentrent melting transition. The minimum should be approximately located at a wall separation H for which the average particle diameterā = σ. Comparing Figs. 3 and 4B, however, shows that this minimum is shifted and placed aroundā = 1.1 σ. To rationalize this observation we consider the free energy distribution F fl (a) calculated with fundamental-measure theory, which indeed exhibits a shifted maximum around a ≈ 1.07 σ as discussed in Sec. II B. Additionally, the differences between commensurate and incommensurate packing in the slab leads to an increase of the liquid free energy for H ≈ 2.5 σ (see Fig. 1).
Using the semi-quantitative theory we can therefore explain the emergence and location of the reentrant melting transition. The depth of the minimum in the critical packing fraction φ c is, however, exaggerated in the theory, which also manifests itself in deviations from the expected dispersity in Fig. 4 B). That being said, the theory shows a near-quantitative agreement with the simulation results for all investigated quantities.
This analysis shows that similar to the observations in bulk [36,39], the crystallization of disperse particles in confinement relies on fractionation. In both cases we would expect that the free-energy barriers that have to be overcome for the particles to diffuse and order would be significantly larger compared to the case of monodisperse crystallization, which could lead to long nucleation times. For high dispersity this is likely to reach time scales much longer than typical laboratory experiments. However, in confinement, the onset of ordering has already been observed in standard EDMD simulations (see the long wavelength peak in the structure factor for H = 2.30σ shown in Ref. [4]). This is likely due to the strong enhancement of heterogeneous nucleation rates compared to homogeneous nucleation [42,43]. Additionally, as has been discussed for bulk systems [39], one possible strategy for an experimental realization is the usage of softer potentials which strongly enhance diffusion and thus speed up the crystallization [65]. We thus believe that the described confinement-induced crystallization is also observable in laboratory experiments. The observed phenomenon is particularly important also for other applications involving external shear stresses. It has been widely observed that shear induces fractionation in channels even at moderate packing fractions, which is for example an important phenomenon for food processing [66]. We expect that similar fractionation effects arise in dense packings, through a combination of shear-induced crystallization [67] and the confinement-induced fractionation detailed here.

IV. CONFINEMENT-CONTROLLED FRACTIONATION IN A WEDGE GEOMETRY
We also use event-driven molecular dynamics to simulate size-disperse particles in a wedge geometry. To create the wedge one of the walls is tilted by a small angle α creating a linear height profile H(x) = H 0 + tan(α)x. The wedge is closed at H(x) = 4.0 σ by a vertical flat wall and only the y-direction is periodic with L y = 60 σ. The wedge consists of 25,000 particles (for α = 1 • ) leading to an average packing fractionφ = 0.51. We use N s = 100, 000 Monte Carlo moves per swap event and employ a similar equilibration time as used in the slab geometry. For the purpose of evaluation the wedge is separated into slices of size ∆H = 0.05 σ and local averages are calculated to determine the profiles shown in Fig. 5. We perform simulations for four independent wedges per tilt angle α to estimate the uncertainty of the presented profiles.
We have seen that confinement has a strong influence on the particle-size distribution. Different from the fractionation observed in bulk systems [36], there is a systematic dependence of the average particle diameter on the wall separation. We can highlight this by tilting one of the walls by a small angle α = 1 • , thus creating a wedge geometry, seen in Fig. 5A [6, 63]. The non-parallel walls lead to a variation in the local wall separation H(x) which now depends on the x-coordinate. The average packing fraction is set toφ = 0.51. After equilibration, separate liquid and crystal regions are observed in the wedge, displaying the same phase behaviour as discussed for the . C)-F) Show local averages of mean particle diameterā, local dispersity δ l , number of crystalline neighbors per particle Nc and packing fraction φ. The averages and errors were determined from four independent runs of the same system. slab geometry (compare Figs. 3 and 5B). The results are independent of the tilt angle α and the inclusion of a liquid reservoir (see App. B). We find that the average particle diameterā exhibits a pronounced sawtooth-like dependence on the local wall separation H(x) while the dispersity attains values between δ l ≈ 0.05 − 0.07 in the crystalline regions (number of crystalline neighbors N c > 3.0) and δ l ≈ 0.1 − 0.15 in the liquid regions (N c < 1.0) (see Fig. 5C-E).
We also observe large local variations in the packing fraction φ(H) with wall separation both in the liquid and crystalline regions (see Fig. 5F). Two very different factors lead to these local density variations. First, even at low average packing fraction in a pure liquid, oscillations in the local packing emerge due to the layering of the particles [6], which explains the behavior of φ(H) in the liquid regions. Second, the local variation of the average particle diameterā(H) in the crystalline regions leads to mechanical stress in the crystals. Assuming a uniform lattice constant we could expect a quadratic increase in the packing fraction φ(H), sinceā scales with H. However, a close inspection of the snapshot in Fig. 5B) reveals that the system responds to the mechanical stress by bending the lattice "planes" (see App. C), which significantly reduces the variation in packing fraction. As a consequence, for smallā the packing fraction increases approximately linearly in the crystalline regions, and it is nearly constant for largeā.

V. CONCLUSION
We have studied the crystallization of size-disperse hard spheres in a confined geometry. The simulation results uncover solid-to-solid transitions between different crystalline phases and a non-monotonic dependence of the critical packing fraction for the liquid-to-solid transition. The latter implies the existence of a reentrant melting scenario. Moreover, we have revealed a systematic dependence of the particle-size distribution on the wall separation in the crystalline phase. To rationalize theses findings we have presented a semi-quantitative theory with which we can explain all features of the stability diagram and the emergent particle-size distributions.
Due to the universality of the underlying mechanisms which lead to the described crystallization in confinement we expect that our results generally describe the phase behaviour of particles for which excluded-volume interactions dominate. We anticipate that this holds in particular for the confinement-controlled demixing and fractionation. We have therefore devised a generally applicable technique for the purification of size-disperse particles by insertion into a wedge geometry. In the future it would be interesting to extend our studies to laboratory experiments on charged or hard-sphere-like colloids [68,69] or soft mesoscopic particles [65]. Furthermore it would be exciting to analyze the influence of external shear stresses or different boundary conditions on the reported crystallization behavior. For example, it is expected that rough walls would significantly increase the critical packing fraction of the liquid-to-solid transition and should thus enable the existence of a stable glass phase. evaluated with shifted free energy differences, ∆F shift = (−1.80 ± 0.25) k B T to estimate uncertainties originating from the assumptions behind cell theory. It should be noted that the value of 0.25 k B T is a lower estimate of uncertainties in the cell theory and that systematic errors can obviously not be discussed quantitatively. The most important conclusion from this analysis is that while a small shift of the free energy leads to quantitative differences, especially for the critical packing fraction, the qualitative features stay the same. It also indicates that the exaggeration of the minimum in the critical packing fraction (and the emergent maximum in the dispersity) is a systematic error of the fundamental measure theory. This error can also be revealed when comparing the density profile n(z) predicted by FMT to event-driven computer simulations as has already been discussed in Ref. [6]. Here we find that these deviations remain even without any differences in the setting up of the dispersity. We conclude that strong inhomogeneous density profiles due to broken translational symmetry are not perfectly described by FMT and lead to an overestimation of the free energy for incommensurate packing and high packing fraction. That being said, the theory shows very good agreement with simulations, despite the discussed uncon- trolled approximations, and thus delivers several important explanations for the observed phenomena.
Appendix B: Dependence on tilt angle and particle reservoir in a wedge geometry The resulting profiles for the wedge with α = 1 • were presented and discussed in Sec. IV, here we compare this data to a wedge with α = 2.6 • (see Fig. 7). The only major difference between the profiles is the suppression of the honeycomb phase for α = 2.6 • . The reason is that the region in which this phase would self-assemble becomes too small to create stable crystals. Additionally, a significant layering is observed in the packing fraction φ at large H due to the vertical wall at H = 4.0 σ (see Fig. 7D). Apart from these differences, the profiles are practically identical.
A typical laboratory experiment would correspond to a grand-canonical ensemble, because the wedge would most probably be coupled to a "bulk" reservoir. Due to the high packing fraction it is, however, impossible to perform a grand-canonical simulation. Here, we study the effect of the ensemble by including a reservoir in our simulations, adjacent to the wedge. This reservoir has a slab geometry with wall separation H = 3.62 σ, and a volume twice that of the wedge. The packing fraction was chosen such that the wedge itself has the same average packing fractionφ = 0.51 as the ones without reservoir. The reservoir is in a liquid phase withā = 0.985 σ and δ = 0.145 and thus sufficiently close to the parti- cle distribution in the bulk. The results are presented in Fig. 8. Only marginal differences between the wedges with and without a reservoir can be observed. The reservoir increases the available number of big particles which are, as discussed before, more prone to crystallize and thus the separation effect is actually enhanced. We thus expect that the effects presented in this manuscript are independent of the ensemble and also hold in the thermodynamic limit, corresponding to a wedge with L x → ∞.
Appendix C: Curvature of the lattice planes in a wedge geometry We also discussed in Sec. IV the role of local variations in particle diameterā(H) inside the crystalline phases which requires the lattice "planes" to exhibit a certain curvature. For a one-layered hexagonal structure under the mild assumption thatā(x) = H(x) we can derive from straightforward geometrical calculations a radius of curvature, R, for lattice "planes" in the crystalline structures. The sketch in Fig. 9A highlights all necessary quantities and a circular sector with the desired radius of curvature R. We find, . (C1) From this we can estimate the radii of curvature R α for the lattice planes marked in Fig. 9B,C under the assumption that H is approximately constant and find R α=1 • ≈ 70 σ and R α=2.6 • ≈ 30 σ in very good agreement with the measured values R m α=1 • ≈ 75.4 σ and R m α=2.6 • ≈ 29.2 σ, respectively.