Phase separation dynamics in a symmetric binary mixture of ultrasoft particles

Phase separation plays an role in determining the self-assembly of biological and soft-matter systems. In biological systems, liquid-liquid phase separation inside a cell leads to the formation of various macromolecular aggregates. The interaction among these aggregates is soft, i.e., these can significantly overlap at a small energy cost. From the computer simulation point of view, these complex macromolecular aggregates are generally modeled by the so-called soft particles. The effective interaction between two particles is defined via the generalized exponential potential (GEM-n) with n = 4. Here, using molecular dynamics simulations, we study the phase separation dynamics of a size-symmetric binary mixture of ultrasoft particles. We find that when the mixture is quenched to a lower temperature below the critical temperature, the two components spontaneously start to separate. Domains of the two components form, and the equal-time order parameter reveals that the domains grow in a power-law manner with exponent 1/3, which is consistent with the Lifshitz-Slyozov law for conserved systems. Further, the static structure factor shows a power-law decay with exponent 4 consistent with the Porod law.


I. INTRODUCTION
In many biological and polymeric systems liquid-liquid phase separation in binary (or more component) systems plays an eminent role in triggering self-assembly processes and the morphology of systems (see, e.g., [1,2]).In biological cells, such phase separation scenarios are induced by the differences in the molecular affinity or chemical potential of the species, thus macromolecules condense into high-and low-density phases (see, for example, [3,4]).These biological systems are in general characterized by a high density and the effective interaction among the constituent molecular assemblies can be considered as soft, allowing thereby significant spatial overlap between the entities [5,6].
In computer simulations, such mutually penetrable, complex macromolecular aggregates can conveniently be modeled via so-called ultrasoft particles, i.e., potentials that are characterized by a finite energy penalty at zero separation.These effective interactions emerge as the degrees-of-freedom of the constituting microscopic particles (such as atomic or molecular entities) are traced out [5,6].Among these ultrasoft interactions the so-called generalized exponential model (with index n) [7]-GEMn (with typically n ≥ 2) -ranges among the most familiar ones, due to the simplicity of its functional form and its wide-spread use in various investigations (see, e.g., [7][8][9][10][11][12][13][14][15].For n > 2 GEM-n particles form -despite their mutual repulsion -stable clusters of overlapping particles: these clusters are polydisperse in the disordered phase and condense upon increasing the densities in stable cluster crystals, i.e., BCC or FCC lattices, where each lattice site is populated by a well-formed cluster of overlapping particles.The phase diagram of such systems (notably for n = 4) is well-explored both at high [7][8][9] as well as at low * gerhard.kahl@tuwien.ac.at temperatures [13].The mechanism behind this unusual phase behaviour and of various other interesting properties [11] is well-understood and documented [12].The particle dynamics in the crystalline state is even more intriguing as particles keep hopping from one cluster to another, preserving the overall crystalline structure [16].Thus cluster crystals represent an interesting model for defect-rich clusters with interesting structural and mechanical properties (see, e.g., [14,15,17].
Ever since the publication of Ref. [7] a few attempts have been made to extend investigations of such a system to more component mixtures of ultrasoft particles.Molecular dynamics simulations have revealed that a polydisperse mixture of ultrasoft particles shows a cluster glass phase: the structure of the clusters (represented via the positions of their centers of mass of clusters) remains disordered and their dynamics becomes drastically slowed down, both features that represent signatures of the emergence of a glassy phase [18].On the other hand, a size-asymmetric binary mixture of ultrasoft particles shows micro-segregation of the two components in the clusters [18].Finally, a few investigations on the ordered phases of binary mixture can be reported (see, e.g., , despite clear evidence of the emergence of cluster crystals for such systems [19,20]. However, phase separation in these mixtures has not yet been studied, a phenomenon which is of utmost importance in understanding the self-assembly occuring in many biological and polymeric systems induced by phase separation.Also, the kinetics of phase separation, i.e., the evolution of the dynamics of the system when rendered thermodynamically unstable by a sudden change of temperature, has not been studied, yet.Several theoretical, simulation, and experimental studies have been dedicated to understand phase separation kinetics in polymer blends and soft-colloid mixtures (see, e.g., [21][22][23]).However none of these have dealt with the phase separation kinetics of a mixtures of ultrasoft particles for which interesting phenomena can be expected due to their cluster-forming capacity.
In this contribution we study the phase separation dynamics of a size-symmetric, equimolar binary mixture of ultrasoft particles (with species A and B).Using molecular dynamics simulations [24] in the canonical ensemble and considering a fairly large number of particles, we determine the tentative phase diagram of this mixture in the temperature-density plane and identify the coexistence curve and an estimate for the critical temperature T c .At high temperatures, the two components of the mixture remain in a homogeneously mixed state, while below T c the species spatially separate into A-rich and Brich phases.The two phases are identified by recording the spatial distribution of the difference in the concentrations of A and B particles [25].When quenched from a high-temperature mixed state to a subcritical temperature, the mixture spontaneously phase separates, and domains of A and B-rich regions are formed.In an effort to quantify this dynamic process we correlate at the same instant the local concentrations at two locations in the system which are separated by a distance r via a correlation function C(r, t).This function provides evidence of a self-similar growth (when scaled with a time-dependent length l(t) which estimates the size of the emerging clusters) which shows asymptotically for large distances r a power-law growth with an exponent equal to 1/3 [26].This feature is consistent with the predictions of the Lifshitz-Slyozov law for conserved systems [27] and has been extensively studied for Lennard-Jones fluids [28][29][30][31].Furthermore, the spatial Fourier transform of this correlation function, a time-dependent structure factor shows at high wave vectors a power-law decay with an exponent equal to 4, a feature which is consistent with Porod's law which signifies the presence of sharp interfaces in the system [32][33][34].
The remainder of the manuscript is organized as follows: in Section II, we introduce the model and provide details about the simulation method and the related protocols.Results are presented and discussed in Section III, while the final section, Section IV, contains the summary of results, concluding remarks, and an outlook to related future investigations.

II. MODEL AND SIMULATION DETAILS
For the interactions between the ultrasoft particles we have used the generalized exponential (GEM-n) potential, whose functional form is given by (see also [7]) In accordance with previous publications we choose n = 4.In Eq. (1) r ij represents the distance between particles i and j (with i, j = 1, ..., N ), N being the total number of particles in the system.The indices α and β can assume the values A and B, representing the two species of the particles.ǫ αβ are the energy parameters and σ αβ are the range parameters of the respective interactions.As we deal with a (size-)symmetric mixture we impose that σ AA = σ BB = σ AB ≡ σ.Further, ǫ AA = ǫ AB ≡ ǫ and ξǫ = ǫ AB with ξ = 1.5.Further N = N A + N B , N A and N B (= N A ) being the number of particles of species A and B, respectively.In our simulation-based investigations the potential is truncated at a distance r c = 2.2σ with Φ(r c ) ≃ ǫ 6.7 × 10 −11 .All results presented in this contribution are based on simulations of ensembles with N = 524 288 particles; we have deliberately chosen such a huge ensemble in an effort to reduce size effects as much as possible.The cubic volume V of the system (with box length L) is chosen in such a way that the (reduced, dimensionless) density, defined by ρ * = ρσ 3 = σ 3 N/V = 2. Time t, temperature T , and density ρ are given in units of σ m/ǫ, k B T /ǫ and ρσ 3 , respectively, where m is the mass of the particles and k B is the Boltzmann constant.For simplicity we set henceforward ǫ, σ, m, and k B to unity.
To simulate the system, (non-equilibrium) molecular dynamics (MD) simulations have been performed in an NVT ensemble, taking advantage of the LAMMPS package [24].The temperature of the ensemble is maintained via a dissipative particle dynamics (DPD) thermostat, which is equivalent to the non-conservative portion of a DPD force field [35].Within DPD, the equations of motion of the particles are given by, In the above relations the r i and p i are the positions and the momenta of a particle with index i.The forces acting on atom i due to other atoms (with index j) are given by the sum over the conservative forces F C ij , the dissipative forces F D ij , and the random forces F R ij : the F C ij can be calculated from the interparticle potential defined in Eq. (1).Further, the F D ij and the F R ij are given by, is the distance between particles with indices i and j, rij = r ij /|r ij | is the unit vector connecting these two particles, and is the relative velocity between particles i and j.Further, γ is the friction coefficient, which is set to unity.ω D (r) and ω R (r) are distance-dependent weight functions vanishing for r > R c [36,37].The usual choice of the weight functions for the continuous stochastic differential equation within the DPD algorithm is given by, For simplicity we choose the cutoff radii R c for these functions as R c = r c .Further, the Θ ij = Θ ij (t) represent uniformly distributed random numbers with a Gaussian statistics, i.e., Θ ij (t) = 0 and Θ ij (t)Θ kl (t ′ ) = (δ ik δ jl + δ il δ jk )δ(t − t ′ ), with the brackets denoting ensemble averages.Initially all the A and B particles are randomly distributed within our simulation box; due to high density of our system it is very likely that the particles overlap.From this initial configuration we start our simulation and equilibrate it during 10 6 MD steps ∆t at a temperature T = 3.0; for the time increment we have used ∆t = 0.005.The equilibrated system is then simulated over 10 6 ∆t, storing configurations in intervals of 10 5 ∆t.All these configurations serve in the following as independent initial configurations in subsequent MD runs.The system is equilibrated over 10 6 ∆t at the respective temperatures T .Eventually the equilibrated systems are simulated over 10 6 ∆t in order to obtain the required information of the system.

III. RESULTS
In this Section, we present our results for the system (as specified in the preceding Section II).We will start by characterizing the morphologies of the system that we have obtained across different temperature ranges; we will then proceed to a discussion on the impact of a quenching process (from high to a low, sub-critical temperature) on the properties of the system.

A. Equilibrium morphologies
We first examine the equilibrium morphologies of the system as they occur at different temperatures.In Fig. 1, we show snapshots of the system in different perspectives views (as specified in the caption) at three different temperatures, T = 1.8, 1.4, and 1.0.At T = 1.8 and 1.4, the two components of the system form a homogeneously mixed phase, while at T = 1.0, a clear spatial separation into an A-and a B-rich phase becomes visible.
In an order to characterize the structure of the system, we have calculated the radial distribution functions g αβ (r) (α, β = A or B) [38] of the two species, which are defined as, where the brackets denote an ensemble average and the prime indicates that the contribution is not considered in case α = β and i = j.The g αβ (r) are depicted in Fig. 2 as functions of distance r for two temperatures (as labeled) in the homogeneous, mixed phase (cf.panels (a) and (b) in Fig. 1  the corresponding snapshots).The g αβ (r) show the typical behaviour as encountered in a liquid phase: the functions oscillate at small r-values and then tend towards unity at large distances.At small distances (i.e., for r 1), the g αβ (r) assume finite values (close to unity), indicating the presence of clusters of overlapping particles, a typical feature occurring in systems of ultrasoft particles.
As systems of ultrasoft particles are prone to form clusters at any density [7-9, 11-14, 16, 39], we have to properly specify what we mean with a cluster in the disordered phase.To this end we have defined a distance r cl , being the position of the first (local) minimum in the g αβ (r) -see Fig. 2 -as a rough estimate for the spatial extent of a cluster; this distance amounts to r cl ≃ 0.6.For all particles pertaining to the same cluster we then define the center of mass, r CM , of this clusters of overlapping particles via, In contrast to the partial radial distribution functions g αβ (r), the radial distribution function of the clusters, g c (r), plotted in the insets of the panels of Fig. 2, vanishes at small distances, akin to the radial distribution function of a system with steeply repulsive interactions [38].This function indicates that the centers of mass of the clusters do not overlap, i.e., that the clusters behave essentially as strongly repulsive, effective particles.

B. Phase separation and critical behavior
In an effort to characterize the morphologies of the system (as they have already been shown in Fig. 1 for different temperatures) we subdivide the simulation box into cubic cells of length 2σ and count the number of A-particles (N i A ) in the cell with index i and relate this number to the total number of particles contained in this cell, i.e., (N i A + N i B ); the fraction of these two numbers is the local concentration Ψ i (which might also be termed Obviously, Ψ i ∈ [0, 1].A typical color-coded map of Ψ i is shown in Fig. 3 for selected snapshots, taken at three different temperatures (as labeled): the upper panels show again perspective views of the simulation box, while the lower panels display cross sections of the respective snapshots, taken at z = L/2.
In an effort to learn more about the phase separation of the system into A-and B-rich phases, we have calculated the probability distribution of the Ψ i in the system at different temperatures, which we term P (x A ).This is done by subdividing the range [0, 1] into 500 bins and by recording the number of occurrences of Ψ i in each of these pin; in this manner we obtain -after normalizing -P (x A ).The function P (x A ) is shown in the top panel of Fig. 4, as obtained for different temperatures (as labeled).Obviously, in the homogeneous phase (occurring at high temperatures), P (x A ) shows a pronounced, symmetric peak centered at x A = 0.5.In contrast, in the heterogeneous regime (i.e, at low temperatures), two peaks appear (which are located at x A -values symmetric to x A = 0.5): here, the peak at the smaller x A -value corresponds to the A-poor (or B-rich) phase, while the peak at the larger x A -value corresponds to the A-rich (or B-poor) phase.The fact that the distribution function shows two pronounced peaks (and vanishes in between) for the lowest temperature (i.e., at T = 1.0) provides evidence that the system is clearly phase separated into regions, hosting either a B-rich (peak at x A ≃ 0.06) and an A-rich phase (peak at x A ≃ 0.94).
If we then plot the x A -values of the peaks in P (x A ) for the different temperatures we obtain the phase diagram, shown in the bottom panel of Fig. 4. Below the critical temperature, T c , we see the characteristic coexistence line, separating the two above mentioned phases.The critical temperature T c can be extracted from these data by fitting the difference in the concentrations of the two coexisting phases (indices '(1)' for the B-rich and '(2)' for the A-rich phase, respectively), i.e., (x A ) (as obtained from the positions of the maxima of the distribution function P (x A ) with the usual functional form [25]: here, B and T c are fitting parameters, anticipating that the phase separation scenario at hand pertains to the 3D-Ising universality class, i.e., β ≃ 0.325.The inset of the bottom panel of Fig. 4 provides evidence that (x A − x

C. Phase separation dynamics
We proceed by exploring the non-equilibrium dynamics of the mixture when quenched from a high temperature to a low, subcritical temperature.We start off from the equilibrated mixture at T = 3.0 which is in a homogeneous phase.We instantaneously quench the system down to T = 0.68 (≃ 0.5T c ), i.e., to a state located inside the coexistence region of the phase diagram.The two components of the system demix, forming A-rich and B-rich domains.In an effort to visualize the emerging morphologies, we assign a value +1 to A-rich (sub-)cells (introduced above) and a value −1 to B-rich (sub-)cells.Snapshots of the mixture and their coarse-grained versions (with binary values) at T = 0.5 T c are shown in Figs. 5 and 6, respectively.In each of the figures panels (a), (b), and (c) show the snapshot of the mixture in after different time-spans after the quench (as labeled).The respective top panels show the the snapshots in perspective views, while the bottom panels are cross sections of the respective upper panels, taken at z = L/2.We track the spatial evolution of these domains with time t by calculating the correlation function C(r, t) of the parameter Ψ(r, t) (previously introduced as Ψ i )r being the position vector of the sub-cell within the simulation box -, evaluated at the time t; where the time-scale starts once the quench of the system has been realized; the above function is defined as: This function correlates the quantity Ψ(r, t), evaluated for two different sub-cells within the simulation box at the same time t: R and (R + r) are the positions of two sub-cells and r is the vector between these two cells points.
The above correlation function C(r, t) = C(r, t), calculated for different values of t (namely t = 100, 300 and 500) is shown in the top panel of Fig. 7 at a temperature T = 0.5 T c .C(r, t) is displayed as a function of r (with the range in r being limited by the simulation box).The correlation function initially decays linearly for all t-values considered, then fluctuates around zero, and eventually tends towards zero.
The correlation functions collapse -when their rdependence is appropriately scaled via a time-dependent scaling factor l(t) -on a single master curve, as shown in the bottom panel of Fig. 7. Thus, introducing thereby the function g(r).The spatial scaling factor for the distance is the (time-dependent) average size of the domains, l(t), which can be estimated as follows: a measure for l(t) is obtained by identifying the first zero of the function C(r, t) (see, e.g., [28]).This length is shown in Fig. 8 as a function of t.For t 80 the length l(t) shows a power-law behaviour with an exponent 1/3 [41,42], a behaviour which is consistent with the theoretical predictions for conserved systems, found, e.g., by Lifshitz and Slyozov [27].
We further define the spatial Fourier transform of the correlation function C(r, t] via S(k, t) = S(k, t) can be considered as a time-dependent structure factor.d(= 3) represents the dimension of the system Further, k is the wave vector and g(p) is the Fourier transform of g(x), i.e., The panels of Fig. 9 show the function S(k, t), either as a function of k -top panel -or of the scaled wave vector, i.e. kl(t) -bottom panel; S(k, t) is shown for three different values of time t (as labeled), evaluated at the temperature T = 0.5T c ; note the semi-logarithmic representations.The unscaled functions show pronounced peaks at k-values, which indicate (via their respective inverse values) the size of the above mentioned domains.Replotting S(k, t) as a function of the scaled wave vector, i.e., of kl(t), we find that the functions, calculated for different t-values collapse on a single master curvesee bottom panel of Fig. 9, showing thus essentially the function g(kl(t)) as defined above in Equation (13).Further we observe that the S(kl(t), t) show for kl(t) 1.6 a power-law type decay with an exponent 4(= d + 1).This observation is consistent with Porord's law [32], which predicts for larger k-values a decay of the scaled function S(kl(t), t) with an exponent (d + 1) signifying the presence of sharp interfaces in the mixture.

IV. CONCLUSIONS AND OUTLOOK
With the help of extensive molecular dynamics simula-(both with respect to ensemble size and simulation length), carried out in an NVT ensemble we have studied the liquid-liquid phase separation scenario in an equimolar, size-symmetric binary mixture of ultrasoft particles (with species labeled A and B).Consistent with previous publications we use the generalized exponential model of index n = 4 for the interparticle interaction.In an effort to keep the temperature fixed we use a DPD thermostat, inbuilt in the LAMMPS program package.While the interactions among like particles are identical, the cross interactions is by a factor of 1.5 stronger, i.e., more repulsive, inducing thereby the desired phase separation scenario.
In a first step we have investigated the equilibrium morphologies of the system via static correlation functions and trace out the phase diagram of the system via calculating and analysing the distribution functions of the local concentration of particle species A. An analysis of these distribution functions in combination with an extrapolation of the coexistence densities (and anticipating that the phase separation scenario pertains to the Ising 3D universality class) yields a critical point at a temperature T c = 1.35.Snapshots of the system confirm the expected behaviour at super-critical, at close-to-critical and at sub-critical temperatures.The structure of the system is analysed for a wide range of temperatures (both above and below T c ) in terms of the (partial) radial distributions and the radial distribution function of the centers of mass of the clusters of overlapping particles formed in the system.At supercritical temperatures the radial distribution function of the clusters shows the expected behaviour: a finite value at short distances, indicating the formation of clusters of mutually overlapping particles as well as oscillations at long distances.The cluster correlation function shows a behaviour akin to a system of strongly repulsive particles, indicating that the clusters behave as effective, mutually repulsive particles.
In a further step we have explored the dynamics of the phase separation process within the system.To this end we have rapidly (i.e., instantaneously) quenched the ensemble from a high temperature to the subcritical temperature T = 0.5T c .In a time-dependent process, the initially homogeneously mixed system then phase separates into A-and B-rich phases, a process which is tracked over time via the time-and space-dependent correlation function C(r, t) of the local concentration: to this end we have subdivided the simulation box in small cubic subcells and have recorded the local concentration of either species.This function correlates at the same instant the concentrations in two sub-cells, separated by a distance r.As the properties of the system is tracked along the time scale, a time-dependent (average) size of clusters of particles, l(t), which grows (from intermediate times onwards) in time via a power law (consistent with the predictions of Lifshitz and Slyozov).Scaling the distance between two sub-cells via l(t) makes the correlation functions calculated for different time instants collapse on a single, timeindependent master curve.Likewise, the time-Fourier transform of C(r, t), the time-dependent structure factors collapse -when the wave-vector is scaled by the cluster size -on a single master curve which shows at large wave-vectors a power-law decay, as predicted by Porod's law.
With the knowledge achieved from this particular, symmetric case we are ready to push forward our in-vestigations on binary mixtures of ultrasoft particles in two directions: on one side we plan to introduce asymmetry into the system, namely both with respect to the size of the particles as well as to the interaction strength.On the other side we will expose the system to other non-equilibrium conditions, such as the exposure to shear forces.Investigations in both directions are currently on their way.

FIG. 1 .FIG. 2 .
FIG. 1. Snapshots of the system at hand at different temperatures.Upper panels: perspective views of snapshots of the system at (a) T = 1.8,(b) T = 1.4,and (c) T = 1.0.Lower panels: cross sections of the respective upper panels, taken at z = L/2.A-and B-particles are shown in yellow and blue, respectively.

FIG. 3 .
FIG. 3. Color-coded representation of Ψi (as defined in the text) for selected snapshots.Upper panels: perspective views of snapshots of the system at (a) T = 1.8,(b) T = 1.4,and (c) T = 1.0.Lower panels: cross sections of the respective upper panels, taken at z = L/2.The actual value of Ψi ∈ [0, 1] of a cell with index i can be extracted from the respective colorcodes, depicted on the right hand side of the snapshots.

FIG. 4 .
FIG. 4. Top panel (a): probability distribution, P (xA) (as defined in the text) as a function of xA for selected, different temperatures (as labeled).Bottom panel (b): symbols -values of the coexistence densities x (1) A and x (2) A (as identified as maxima of P (xA)) for the different (sub-critical) temperatures investigated in this study.The red line displays the function specified in Eq. (10), using the fitting parameters given in the text.The inset shows x (2) A − x of temperature T ; the dotted line displays again the function specified in Eq. (10).

FIG. 5 .FIG. 6 .
FIG.5.Snapshots of the system at different times t after the system has been quenched down to T = 0.5 Tc (as detailed in the manuscript).Upper panels: perspective views of snapshots of the system at (a) t = 100, (b) t = 300, and (c) t = 500.Lower panels: cross sections of the respective upper panels, taken at z = L/2.A-and B-particles are shown in yellow and blue, respectively.

FIG. 7 .
FIG. 7. Top panel: correlation function C(r, t) as defined in Eq. (11) as a function of r, evaluated for different values of time t (as labeled) after the quench; the considered temperature isT = 0.5Tc.Bottom panel: correlation function C(r, t), shown in the top panel, now as a function of the scaled distance r/l(t), as defined in the text.

FIG. 9 .
FIG. 9. Time-dependent structure factor S(k, t) as defined in Eq. (13) as a function of the wave vector k -top panel -and of the scaled wave vector kl(t) -bottom panel; S(k, t) is evaluated for different values of t (as labeled).The dotted line in the bottom panel indicates a power-law decay in k with ∼ k −4 .Note the semi-logarithmic presentation of the data.