How can dissipation constrain fluctuations beyond equilibrium?
Diffusion, structure and biased energy flows in a driven liquid
Based on a minimal model of periodically driven liquid, we explore how the dynamics and structure can be modified by tuning energy dissipation. We first put forward simple connections between tracer diffusion, density correlations and the rate of work applied by external nonequilibrium forces. Then, we show that atypical structures can be designed by biasing trajectories which preferentially dissipate or absorb energy on targeted particle bonds. Specifically, our analytical results and numerical simulations reveal that such energy flows in the liquid lead to a renormalization of pair interactions combined with spontaneous clustering. This interplay between dissipation and spatial organization can potentially be relevant for complex self-assembly.
Nonequilibrium forces can drive novel and specific pathways to modulate phase transitions and self-assembly in materials. The close connection between energy dissipation, powered by these forces, internal transport and spatial organization is especially apparent in living systems Toyabe et al. (2010); Fodor et al. (2016a); Battle et al. (2016); Gnesotto et al. (2018). As an example, the flagella motors of E. Coli exhibit a unique phenomenology combining ultra-sensitive response, adaptation, and motor restructuring as a function of the applied load Lele et al. (2013); Lan et al. (2012); Wang et al. (2017). Moreover, in vivo studies of the cellular cytoskeleton, as well as in vitro experiments on reconstituted systems, have also shown that motor-induced forces control a large variety of functionality in the cell Soares e Silva et al. (2011); Sanchez et al. (2012); Blanchoin et al. (2014); Murrell et al. (2015); DeCamp et al. (2015).
To elucidate the role of nonequilibrium forces in materials, it is crucial to examine how dissipation affects the emerging dynamics and structure. While equilibrium features are well established, progress in controlling systems with sustained dissipation has been hampered by a lack of general principles Cates and Tailleur (2015); Solon et al. (2015a); Nguyen and Vaikuntanathan (2016); Fodor et al. (2016b); Murugan and Vaikuntanathan (2017); Nardini et al. (2017); Nguyen and Vaikuntanathan (2018). Minimal models of active and driven systems provide analytical and numerically tractable test beds to investigate the interplay between dissipation and material properties far from equilibrium Marchetti et al. (2013); Han et al. (2017); Bechinger et al. (2016); del Junco et al. (2018); Fodor and Marchetti (2018). Such models have demonstrated how nonequilibrium driving can induce phase transitions and excite novel collective responses in soft media Vicsek et al. (1995); Tailleur and Cates (2008); Han et al. (2017); Nguyen and Vaikuntanathan (2016); van Zuiden et al. (2016). To rationalize such a phenomenology, recent theoretical work has proposed extending equilibrium concepts to active media, such as the definition of pressure Takatori and Brady (2015); Solon et al. (2015a, b). Besides, others have shown that changing dissipation strongly affects the internal transport and the density fluctuations of nonequilibrium liquids Cagnetta et al. (2017); del Junco et al. (2018); Nemoto et al. (2018), suggesting that controlling dissipation is a fruitful route to constraining material properties. Yet, despite these advances, any generic principle between dissipation and spatial organization still remains elusive.
In this paper, we explore how dissipation affects the dynamics and structure of a periodically driven liquid. Building on recent work Han et al. (2017); del Junco et al. (2018), we consider an assembly of Brownian particles, with only a subset being driven by an external time-dependent force :
where if and otherwise. The driven particles belonging to the set are referred to as tracers. We mainly focus on instances where the fraction of driven particles is lesser than the fraction of undriven particles, so that undriven particles are referred to as bath particles. The fluctuating term is a zero-mean white Gaussian noise with correlations , where and respectively denote the damping coefficient and the bath temperature, with the Boltzmann constant set to unity . The driving force is taken in two dimensions as
where and are respectively the amplitude and the frequency of the drive, such as the drive persistence reads . The relative strength of the drive is given by the Péclet number , where is the typical particle size Han et al. (2017); del Junco et al. (2018). In the absence of interactions (), the average position of driven tracers follows a periodic orbit, describing a circle in two dimensions.
In Sec. II, we first connect dissipation to liquid properties. Specifically, we consider the diffusion coefficient of a tagged driven tracer particle, and the density correlations between tracer particles and undriven bath particles in the liquid. At variance with del Junco et al. (2018), our relation between dissipation and diffusion is based on an explicit coarse-graining in terms of microscopic details Dean (1996); Démery and Dean (2011); Démery et al. (2014). Besides, we show that deviations from equilibrium tracer-bath correlations hold a quantitative signature of dissipation. These results clarify how dissipation modifies the transport and structure of the liquid, yet they do not provide any concrete intuition on how to stabilize particular configurations in terms of dissipation.
In Sec. III, we explore the interplay between energy flows and resultant organization in a manner independent of the details of the driving protocol. To this end, we sample trajectories which preferentially pump energy into, or extract energy from, interactions between specific classes of particles. The ensemble generated by such a sampling serves as a proxy for cases where energy is injected into particular modes by external forces. Using large deviation techniques Garrahan et al. (2007); Hedges et al. (2009); Jack and Sollich (2010); Pitard et al. (2011); Speck et al. (2012); Bodineau et al. (2012); Chetrite and Touchette (2013); Limmer and Chandler (2014); Nemoto et al. (2017), we predict perturbatively how interactions between the particles are modified in this ensemble. Direct sampling of biased trajectories, based on cloning algorithm Giardinà et al. (2006); Tailleur and Kurchan (2007); Hurtado and Garrido (2009); Nemoto et al. (2016); Ray et al. (2018); Klymko et al. (2018); Brewer et al. (2018), confirms that tuning energy flow into desired interactions changes the liquid structure in a controlled manner.
Together, our results illustrate how dissipation, induced by either external forces or energy bias, allows one to tailor specific structures in soft materials. Such a control can potentially result in design principles, for the construction of complex addressable self-assembly landscapes, at the cost of energy dissipation.
Ii Dissipation and liquid properties
To connect tracer dissipation and liquid properties, we first describe the dynamics of undriven particles in terms of a coarse-grained variable. Using standard techniques, the dynamics of the density field can be written as a non-linear Langevin equation Dean (1996). In the regime of weak interactions, the density fluctuations around the average density are Gaussian in terms of the following Hamiltonian Chandler (1993); Démery et al. (2014); Krüger and Dean (2017):
where . The conserved density dynamics reads
where and are respectively the field diffusion coefficient and the field damping coefficient. The term is a zero-mean Gaussian white noise with correlations . Owing to the linearity of the density dynamics (4), it can be readily written in Fourier space as
so that the field dynamics can be directly solved as
Considering the limit of dilute driven tracers, where interactions among them are negligible, their dynamics reads
with and referring to spatial dimension. As a result, (6) and (7) provide a closed dynamics for tracers only. It should only be valid for weak interactions a priori, yet previous works have shown that it remains qualitatively relevant even beyond this regime in practice Démery (2015); Martin et al. (2018). Indeed, Gaussian field theories for density fluctuations provide a very good description of simple liquids Chandler (1993).
To characterize the transport properties of the liquid in the presence of the driving forces, our first goal is to obtain an explicit expression, in terms of microscopic details, for the tracer diffusion coefficient:
We aim to explore connections between and dissipation, which is defined from stochastic thermodynamics as the power of the forces exerted by all tracers on solvent: . Here denotes here a Stratonovich product Sekimoto (1998); Seifert (2012). Substituting the dynamics (1), the dissipation can be separated into free-motion and interaction contributions as , where the rate of work reads
Note that we have used to simplify for a symmetric . Given that is the only non-trivial contribution to dissipation, connecting diffusion and dissipation simply amounts to expressing in terms of .
To set up a proper perturbation scheme, we scale the pair potential with a dimensionless parameter , which controls the coupling between tracer and bath equations of motion. In Appendix A, we obtain some explicit expressions for and to quadratic order in together with various limits of external driving parameters, such as amplitude or frequency. We first consider the limit of high frequencies , where the relaxation time scale is set here by density diffusion over tracer size , and of small Péclet number . In this regime, the rate of work and the increase of the diffusion coefficient , where is the equilibrium diffusion, scale like
These results are valid to quadratic order in and in . In the opposite limit, namely at low frequencies and small Péclet number , we get
Both and are now independent of the driving frequency . Again, this is valid to order and . Note that the scaled work coincides with the reduced equilibrium diffusion to this order Démery and Dean (2011); Démery et al. (2014), as expected from linear response. Our perturbation theory hence shows that scalings of and are identical, both in terms of the drive amplitude and of its frequency .
When the size of the bath particles is significantly smaller than the tracer size , which amounts to setting different pair potential for bath-bath and for bath-tracer interactions, one can safely neglect the variation of in (10-11), such as . Then, in both regimes and , the renormalization of the diffusion coefficient can be simply written in terms of the rate of work for as
Thus, the excess rate at which tracers move over their own size compared to equilibrium, set by the lhs of (12), is controlled by the rate at which work is applied on tracers, set by the rhs of (12). The proportionality factor is determined by the details of interactions and of the density fluctuations. This result corroborates some numerical observations obtained previously for a similar system, where composition-dependent diffusion constants can lead to phase transitions del Junco et al. (2018).
We now explore how dissipation can be connected to density correlations. To this end, we treat undriven bath particles without any approximation in what follows, instead of relying on the Gaussian density field theory for as we did above. Using Itô calculus, the average rate of change of the bath-tracer interaction potential can be written as
where denotes here an Itô product. Substituting the dynamics (1) and using , we get
where the -sum runs over all particles. From the steady-state condition , we then deduce that the rate of work , defined in (9), can be written in terms of density correlations as
and denotes a sum without the overlap of indices: and . The power balance (15) reflects how density correlations adapt to the presence of nonequilibrium forces. For a vanishing rate of work (), one recovers the first order of the equilibrium Yvon-Born-Green hierarchy, in its integral form, for two-component fluids Hansen and McDonald (2013). At finite rate of work (), the relation between the two-body correlation and the three-body terms is now implicitly constrained by dissipation.
Though being an exact result, (15) is not straightforward to test, either numerically or experimentally, due to the three-body correlations. In equilibrium, where tracer and bath particles are indistinguishable, we get . Assuming that this remains approximately valid in the driven case, for a dilute fraction of tracers, the rate of work can simply be written in terms of the force extered by bath particles on a tracer as
where we have used in the dilute limit. We simulate the dynamics (1) with the WCA potential , where denotes the Heaviside step function Weeks et al. (1971). Our measurements in Fig. 1 show that (17) is indeed a good approximation at small and small , namely when the drive only weakly perturbs the liquid. Besides, up to and , the relative deviation is at most of about . At variance with previous approaches Harada and Sasa (2005); Lander et al. (2012); Battle et al. (2016), which rely on probing the whole system, our results demonstrate that dissipation can actually be evaluated, with only a small error, by considering solely bath-tracer forces: the contribution of interaction forces among other particles is negligible for a dilute fraction of driven tracers.
To further evaluate the change in liquid structure induced by dissipation, we measure the deviation from equilibrium pair-correlation due to the driving forces. At a given , scaling by reveals that all curves almost collapse into a master curve for our numerical range , as shown in the middle column of Fig. 2. Given that the rate of work also exhibits such a scaling, this suggests the existence of an underlying relation between and : it appears that their ratio indeed gives a factor close to , as depicted in the right column of Fig. 2. Hence, this empirical relation shows that, considering the periodic drive (2) in the dilute limit, the dissipation can be directly estimated by comparing driven and equilibrium pair-correlations: this is clearly an asset with respect to invasive methods based on comparing fluctuations and response Harada and Sasa (2005); Mizuno et al. (2007); Fodor et al. (2015); Turlier et al. (2016); Ahmed et al. (2018).
Overall, the results of this Section illustrate how dissipation affects the transport properties of a driven liquid, measured in terms of diffusion coefficient, as well as its structural properties, characterized by density fluctuations. This immediately suggests that internal transport and emerging structure can be tuned externally by monitoring the nonequilibrium forces which power the dissipation. A related question then is whether specific structures and correlations can be stabilized by biasing particle trajectories which pump energy into or out of specific interactions. As mentioned previously, we regard such a bias as a way to mimic the effect of generic external drive. We now investigate this question.
Iii Interactions in biased ensembles
To explore how energy flows affect the microscopic interactions, we consider biasing the dynamics in terms of the rate of energy stored in the pair potential between specific subsets of particles. To sample such a biased ensemble, we rely on the framework of large deviation theory Chetrite and Touchette (2013); Jack and Sollich (2010). Formally, the biased ensemble is generated by introducing a weighting factor in the dynamical path probability. The trajectories and configurations in biased ensemble are then associated with a given statistics of , which can be tuned by monitoring . Such trajectory biases have already been used in a large variety of contexts, for instance, to explore dynamical heterogeneities in glassy systems Garrahan et al. (2007); Hedges et al. (2009); Pitard et al. (2011); Speck et al. (2012); Bodineau et al. (2012); Limmer and Chandler (2014); Nemoto et al. (2017), soliton solutions in high-dimensional chaotic chains Tailleur and Kurchan (2007); Laffargue et al. (2013) and, more recently, the clustering of self-propelled particles Cagnetta et al. (2017); Whitelam et al. (2018); Nemoto et al. (2018).
In what follows, we consider the effect of biasing the dynamics with where
and . The control parameters are specific for any given pair of particles . They can be chosen so that the rate of change of only specific inter-particle pair potentials are accounted for and used for biasing. Indeed, with the definition (18), we get , where denotes here an average in the biased ensemble. In the unbiased case (), there are no energy flows for any pair of particles, so that . At variance, by tuning positive (negative), one can now generate trajectories that preferentially inject (extract) energy into (from) some specific interactions by effectively realizing ().
Following the procedure in Jack and Sollich (2010); Chetrite and Touchette (2013), an auxiliary physical dynamics, with the same statistical properties as in the biased ensemble, can be constructed by solving the following eigenvalue equation
where the eigenvalue , parametrized by , is the scaled cumulant generating function appropriate to the biasing observable . The operator is the adjoint of the Fokker-Planck operator associated with the dynamics (1):
The force field of the auxiliary dynamics is then defined in terms of the eigenfunction as
In practice, computing is a highly non-trivial procedure for many-body systems, the explicit solutions obtained so far are only either perturbative or restricted to non-interacting systems Chetrite and Touchette (2013); Tsobgni Nyawo and Touchette (2016).
In our case, a simple expression can be obtained for the auxiliary force by solving (19) perturbatively in terms of the bias parameter . Specifically, we expand
where is the average in the unbiased case, and is the uniform eigenvector associated with the zero eigenvalue. Given that in steady state, the leading non-trivial order of (19) reads
As a result, the trajectories that extract (inject) energy from (into) some pair interactions can be generated with a physical dynamics where interactions are simply renormalized, at leading order, with a specific strength for each pair. For systems with purely repulsive pairwise interactions, modifying the interaction strength between targeted particle pairs is qualitatively consistent with the effect of external driving. For instance, phase separation in mixtures of driven and undriven particles, reported both experimentally and numerically del Junco et al. (2018); Han et al. (2017), can be rationalized in terms of an effective decrease of specific interactions between these particles.
Moreover, the effect of higher order bias can be anticipated by considering another biasing observable:
As shown in Appendix C, biasing the original dynamics (1) with is equivalent to biasing the auxiliary dynamics, specified by (24), with . In short, the original bias then results in two combined effects: (i) modifying some specific pair potentials, and (ii) optimizing the total squared force in the integrand of , which effectively tends to cluster particles.
To probe the range of validity of our perturbation, we consider an assembly of passive Brownian particles () where all the pairs between a subset and other particles are biased with the same strength , such as . To make connection with the setup used in the previous Sections, could for instance be the set of tracer particles. With this choice, we get where . In other words, we are biasing the flow of energy into bonds between particles in the subset and all other particles.
We compare measurements of obtained from simulations of the auxiliary dynamics (24) and from a direct sampling of the biased ensemble. The latter is implemented with a cloning algorithm where rare realizations are regularly selected and multiplied for efficient sampling Giardinà et al. (2006); Tailleur and Kurchan (2007); Hurtado and Garrido (2009); Nemoto et al. (2016); Ray et al. (2018); Klymko et al. (2018); Brewer et al. (2018). For convenience, interactions are now given by the soft-core potential . We observe a very good agreement between the two measurements for a range of . In particular, the agreement improves as is decreased, as reported in Figs. 3(a-b). This supports the validity of our perturbation in this regime, up to interaction change of about when . Note that the range of validity is asymmetric in .
To explore further the features of the biased ensemble, we now compare the density correlations of biased pairs obtained from both direct sampling and auxiliary dynamics. At , the agreement is good for the whole curve when , whereas a clear deviation appears beyond when , as shown in Figs. 3(c-d). In both cases, the region of particle overlap is well reproduced, as expected, since it is entirely controlled by the interaction change captured in (24). Yet, the tendency for particles to cluster, manifest in the peak at , comes as a higher order effect beyond this perturbation. Note that the peak value is comparable for , in agreement with (25) being symmetric in . Besides, at fixed , the structural modification becomes more dramatic as increases. Altogether, these results demonstrate that selecting trajectories which preferentially inject (extract) energy into (from) some specific pair interactions, by tuning , indeed modulates the liquid structure in a controlled manner for small bias and weak interactions.
Finally, to determine whether our bias can stabilize configurations which are significantly distinct from the unbiased system, we probe numerically the effect of large bias from direct sampling. The particles spontaneously tend to cluster for both positive and negative , as shown in Fig. 4 and movies in mov (). This confirms the propensity of trajectories to maximize interaction forces at high bias, in agreement with (25). The shape of clusters differs depending on the sign of : a micelle-like structure featuring the particles in at the core (blue) surrounded by others (red) appears for , whereas clusters have a random composition for . Again, this agrees with (24), where the interactions strength is either increased or decreased, respectively for and . Overall, these results establish a reliable proof of principle for the design of self-assembled structures based on biasing energy flows into or out of particular pair interactions.
Developing techniques to characterize and control the behavior of systems operating far from equilibrium remains a central and outstanding problem. In this paper, we have demonstrated that specifying the amount of energy dissipated by nonequilibrium forces allows one to constrain the dynamics and structure of such systems. Despite the apparently complex interplay between internal dissipation and emerging properties, their exist simple connections between tracer diffusion, density correlations and dissipation in a driven liquid. It remains to investigate whether similar results can be obtained in more complex systems, for instance including anisotropic building blocks, such as active liquid crystals or driven chiral objects Joshi et al. (2017); van Zuiden et al. (2016); Nguyen et al. (2014).
Motivated by these findings, we have also put forward a novel method to control liquid structure: it relies on biasing the dynamics with the pair potential between individual components. For a minimal case study, we have shown that the emerging biased configurations, sampled with state-of-the-art numerics Giardinà et al. (2006); Tailleur and Kurchan (2007); Hurtado and Garrido (2009); Nemoto et al. (2016); Ray et al. (2018); Klymko et al. (2018); Brewer et al. (2018), could be rationalized using recent methods of large deviation theory Chetrite and Touchette (2013); Jack and Sollich (2010). Importantly, since our analytic results encompass the case of a specific bias for each pair, it might be possible to promote the spontaneous self-assembly of more complex structures in such biased ensembles. For instance, inspired by recent works Murugan et al. (2015); Stern et al. (2017), one could potentially design appropriate energetic landscapes, in terms of the pair-specific bias parameters, to selectively stabilize some target molecules. Finally, our results extend straightforwardly to active systems Marchetti et al. (2013); Cates and Tailleur (2015); Bechinger et al. (2016); Fodor and Marchetti (2018), where the driving force has now an independent dynamics for each particle, opening the door to control self-organization in biological context Sperandio et al. (2002); Brown and Johnstone (2001); Wedlich-Söldner and Betz (2018).
Acknowledgements.This work was granted access to the HPC resources of CINES/TGCC under the allocation 2018-A0042A10457 made by GENCI and of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the program Investissements d’Avenir supervised by the Agence Nationale pour la Recherche. SV and LT were supported by the University of Chicago Materials Research Science and Engineering Center, which is funded by National Science Foundation under award number DMR-1420709. SV acknowledges support from the Sloan Foundation and startup funds from the University of Chicago. ÉF benefits from an Oppenheimer Research Fellowship from the University of Cambridge, and a Junior Research Fellowship from St Catherine’s College.
Appendix A Dissipation and diffusion
This appendix is devoted to the derivation of the dissipation rate and the diffusion coefficient of a driven tracer, as defined in Sec. II. We employ a perturbative treatment at weak interactions, originally introduced for a particle driven at constant force in Démery and Dean (2011); Démery et al. (2014). To this aim, the tracer-bath interaction potential is scaled by a small dimensionless parameter in what follows. Besides, we focus on the regime of dilute tracers, so that interactions among them, either direct or mediated by the bath, can be safely neglected .
The dynamic action associated with the tracer dynamics (6-7) follows from standard path integral methods Martin et al. (1973); De Dominicis (1975). It can be separated into contributions from the free tracer motion and from interactions, respectively denoted by and :
where is the tracer diffusion coefficient in the absence of interactions (), and is the process conjugated with the tracer position . For weak interactions , any average value can be then expanded in terms of as , where is the average taken with respect to only. As a result, determining the first correction from interactions in any observable amounts to computing the corresponding average .
Considering the dissipation rate , the leading order is , and the first correction reads . Given the explicit form of in (26), the correlations of interest are
where we have used and . Expanding at small , we deduce
Substituting the explicit expression of the drive (2), and integrating over and , we obtain
We now turn to deriving the diffusion coefficient . It is defined in terms of the mean-squared displacement (MSD) as . At leading order, the MSD reads . To obtain the first order, we need to compute the following correlations
where we have used again that is Gaussian in terms of , yielding
Expanding at small , we get
where refers to the diffusion coefficient in the absence of driving force (). From the explicit time integrations, we then obtain
Appendix B Numerical simulations
In Sec. II, numerical simulations of the dynamics (1) are performed using the LAMMPS simulation package in a two-dimensional box , where is the particle diameter, with periodic boundary conditions at average density . The driving force only acts on ten percent of the particles with dynamics (2). The time step is . The system is first relaxed for time steps, and later equilibrated for . We evaluate average values over trajectories with duration of at least . Parameter values: , , , .
In Sec. III, a custom code of molecular dynamics, based on finite time difference, is used to perform the simulations in a two-dimensional box with periodic boundary conditions. The pair potential of particle pairs is biased in a system of particles in total. To sample the biased ensemble, we use the cloning algorithm described in Appendix A of Nemoto et al. (2016). The time interval for cloning is and the number of clones is . The time step is , the initial relaxation time is , and the total simulation time is . Parameter values: , , (Fig. 4), .
Appendix C Equivalence of biased ensembles
In this Appendix, we demonstrate the equivalence between two dynamical biased ensembles: (i) the original dynamics (1) biased with , where is defined in (18), and (ii) the auxiliary dynamics with force field (24) biased with defined in (25). This amounts to showing that the path probabilities associated with (i) and (ii) are similar. The corresponding dynamic actions, respectively denoted by and , are given here in terms of the Lagragians and defined as