Superadiabatic demixing in nonequilibrium colloids

Abstract

Dispersed colloidal particles that are set into systematic motion by a controlled external field constitute excellent model systems for studying structure formation far from equilibrium. Here we identify a unique demixing force that arises from repulsive interparticle interactions in driven binary colloids. The corresponding demixing force density is resolved in space and in time and it counteracts diffusive currents which arise due to gradients of the local mixing entropy. We construct a power functional approximation for overdamped Brownian dynamics that describes superadiabatic demixing as an antagonist to adiabatic mixing as originates from the free energy. We apply the theory to colloidal lane formation. The theoretical results are in excellent agreement with our Brownian dynamics computer simulation results for adiabatic, structural, drag and viscous forces. Superadiabatic demixing allows to rationalize the emergence of mixed, laned and jammed states in the system.

Introduction

The concept of mixing entropy in multi-component systems is profound for our understanding and quantitative description of complex systems. In equilibrium, mixing entropy can overcome strong energetic repulsion that arises due to internal interactions between the constituent particles. It is a driving mechanism, if not an antagonist, for a wide range of structuring and self-assembly phenomena^{1,2,3,4,5,6,7,8,9,10,11,12}. The systems and phenomena where mixing entropy plays a crucial role cover a broad range, ranging from the fundaments of liquid state theory to specific applications.

Mixing entropy was shown to be relevant for the determination of the equation of state of binary mixtures of hard spheres¹. It determines phase behavior, such as entropy-driven phase transitions of colloid-polymer mixtures², liquid-liquid immiscibility in polymer solutions³, and the demixing transition in athermal mixtures of colloids and flexible self-excluding polymers⁴. For the liquid crystal physics displayed by anisotropic particles, mixing entropy plays a role for the isotropic-nematic-nematic phase separation for bidisperse rod-like particles⁵, and for rod-plate mixtures, where biaxiality competes with demixing⁶. Mixing entropy is relevant for percolation in binary and ternary mixtures of patchy colloids⁷, and for selectivity in spatially inhomogeneous binary fluid mixtures⁸.

Important technological applications where mixing entropy is relevant include the capacitive mixing for harvesting the free energy of solutions⁹, and in ‘blue energy’ from ion adsorption and electrode charging in sea and river water¹⁰. Although in nonequilibrium the situation is much less clear, gradients of position-resolved mixing entropy fields were shown to be relevant e.g., for the dynamics of liquid films with soluble surfactant¹¹. Beyond the soft matter realm, mixing entropy is relevant for the stability of metallic alloys¹².

Here, we identify and describe a competing effect that occurs in genuine nonequilibrium and that can counteract diffusive forces generated by the mixing entropy in a similar way that explicit interparticle repulsion does in equilibrium. We show that this effect is “superadiabatic” in character, i.e., it acts above all effects that can be understood on the basis of an equilibrium (“adiabatic”) reference state and its free energy. As we show, superadiabatic demixing is a genuine nonequilibrium effect and a corresponding unique superadiabatic force density distribution can be identified that acts spatially and temporally resolved in nonequilibrium systems. As both a relevant application and a demonstration of the concept we revisit the well-studied phenomenon of colloidal lane formation in oppositely driven binary mixtures, where for the first time we are able to rationalize quantitatively, on the basis of a physical model of the underlying superadiabatic effect, the emergence of nonequilibrium structure formation in this system.

Results

Lane formation

We show in Fig. 1 characteristic snapshots of a colloidal mixture of two species with repulsive interparticle interactions driven in opposite directions. At low driving, Fig. 1a, the diffusive forces generated by the entropy of mixing dominate and the system remains in a homogeneous state with both species flowing through each other. At high driving however the species segregate into two lanes, Fig. 1b. Lane formation constitutes a genuine nonequilibrium self-organization process that has attracted much interest in the literature since it occurs in strikingly different areas, ranging from colloidal systems^{13,14,15,16,17,18}, plasmas^19,20, and lattice models²¹, to different kinds of living systems, such as bacteria in channels²², ant trails²³, and groups of pedestrians²⁴. Laning has also been studied when the external driving directions of the two species are non-parallel to each other²⁵, have high shear rates²⁶, or other characteristics such as rotating magnetic field in channels, and periodic driving with different friction coefficients^27,28,29. Studies were devoted to the influence of noise and of hydrodynamic interactions^30,31, as well as the characteristics of the transition toward laning in two-dimensional systems³². Although the primary focus is on purely repulsive model systems, such as model suspensions of charged colloids, laning has also been shown to appear in systems governed by attractive³³ and by dipolar³⁴ interparticle interactions.

Fig. 1: Characteristic snapshots of a two-dimensional colloidal mixture in a homogeneous (gravitational) external field.

Particles of species 1 (orange) possess a positive buoyant mass, i.e., the external force acting on them is directed downwards. Particles of species 2 (blue) have a negative buoyant mass and hence upwards directed external force acting on them. The arrows indicate the direction of the different forces in the system. The color of the arrow shows the species on which the forces act, and its size indicates the relative magnitude of the force as compared to the other forces. At low external driving (a) the system is homogeneous. Forces only act in the flow $\hat{{\bf{y}}}$-direction. At high external driving (b) the system segregates into lanes. The forces acting in the gradient $\hat{{\bf{x}}}$-direction (diffusive, adiabatic, and structural) balance each other. The small white arrows indicate the direction of the velocity inside the lanes.

Physical mechanisms for laning were attributed to rectification of diffusion on the particle scale³⁵. The pair correlation functions were considered to be the key observables to quantify the laning phenomenon³⁶. Theoretical treatments have been scarce. Based on phenomenological versions of dynamical density functional theory Chakrabarti and coauthors^14,15 obtained a dynamical instability as a steady-state bifurcation of the density field around the onset of structural inhomogeneity in the simulations. Poncet et al.³⁶ used linearized stochastic density functional theory for the pair correlation functions. Kohl et al.³⁷ used a microscopic approach, based on Kirkwood’s superposition approximation as a closure relation, to also calculate nonequilibrium pair correlations in strongly interacting driven binary mixtures. Wächtler et al.³³ performed a stability analysis based on dynamical density functional theory.

Despite this progress, a systematic theoretical framework that would describe laning has not yet been formulated. Here, we present a comprehensive theory for driven binary mixtures that operates on the level of one-body correlation functions, such as (partial) density and current distributions. The general framework constitutes a closed theory, based on the power functional concept for mixtures^38,39. We present an approximate form of the superadiabatic free power functional, which is a functional generator of the nonequilibrium force distributions. Crucially, we demonstrate that the superadiabatic demixing force fields are kinematic functionals, i.e., depend not only on the density, but also on the flow fields. The theory allows us to quantitatively describe the stability of lane formation.

Dynamical one-body description

Let the positions of N particles in d spatial dimensions be indicated by r₁, …r_N ≡ r^N. We consider systems that are composed of different species, labeled by α. The set ${{\mathcal{N}}}_{\alpha }$ contains all particle indices i that belong to species α. The dynamics of the system are governed by the coupled Langevin equations for overdamped Brownian motion,

$$\gamma {\dot{{\bf{r}}}}_{i}=-{\nabla }_{i}u({{\bf{r}}}^{N})+{{\bf{f}}}_{{\rm{ext}},i}({{\bf{r}}}_{i},t)+{{\boldsymbol{\xi }}}_{i}(t),$$

(1)

where γ is the friction constant against the implicit solvent (we consider γ to be the same for all particles), the overdot denotes a time derivative, ∇_i indicates the derivative with respect to r_i, u(r^N) is the interparticle interaction potential, f_ext,i(r, t) is an external force field acting on particle i at position r at time t; here r is a generic position coordinate. Furthermore ξ_i(t) is a stochastic white noise force with vanishing mean, 〈ξ_i(t)〉 = 0, and auto-correlation given by

$$\langle {{\boldsymbol{\xi }}}_{i}(t){{\boldsymbol{\xi }}}_{j}(t^{\prime} )\rangle =2\gamma {k}_{{\rm{B}}}T{\bf{1}}{\delta }_{ij}\delta (t-t^{\prime} ),$$

(2)

where the angles denote an average over the realizations of the noise, k_B indicates the Boltzmann constant, T is absolute temperature, 1 is the d × d-unit matrix, and δ(⋅) is the Dirac delta distribution. In practice, we discretize the equations of motion, Eq. (1), with a finite time step Δt, and integrate in time using the Euler algorithm.

We define the one-body density and current distribution functions, respectively, for each species α as

$${\rho }_{\alpha }({\bf{r}},t)=\left\langle \sum _{i\in {{\mathcal{N}}}_{\alpha }}\delta ({\bf{r}}-{{\bf{r}}}_{i})\right\rangle ,$$

(3)

$${{\bf{J}}}_{\alpha }({\bf{r}},t)=\left\langle \sum _{i\in {{\mathcal{N}}}_{\alpha }}\delta ({\bf{r}}-{{\bf{r}}}_{i}){{\bf{v}}}_{i}\right\rangle ,$$

(4)

where the angles denote an average both over the noise as before, but also over the set of initial states; the velocity of particle i is given by a symmetric time derivative, v_i(t) = (r_i(t + Δt) − r_i(t − Δt)) ∕ (2Δt) (see e.g., appendix A of ref. ⁴⁰ for the derivation). We obtain the microscopically resolved velocity profile of species α as the ratio

$${{\bf{v}}}_{\alpha }({\bf{r}},t)=\frac{{{\bf{J}}}_{\alpha }({\bf{r}},t)}{{\rho }_{\alpha }({\bf{r}},t)}.$$

(5)

The continuity equation holds individually for each species,

$$\frac{\partial {\rho }_{\alpha }}{\partial t}=-\nabla \cdot {{\bf{J}}}_{\alpha },$$

(6)

where ∇ indicates the derivative with respect to r.

The one-body equation of motion translates the sum of all forces that act into local motion. The dynamics on the one-body level are specified by the (exact) force density balance relationship

$$\gamma {{\bf{J}}}_{\alpha }=-{k}_{{\rm{B}}}T\nabla {\rho }_{\alpha }+{{\bf{F}}}_{{\rm{int}}}^{\alpha }+{\rho }_{\alpha }{{\bf{f}}}_{{\rm{ext}}}^{\alpha },$$

(7)

where the three terms on the right hand side constitute the force densities due to thermal diffusion, internal interactions, and external influence. Due to the overdamped character of the dynamics, the sum of these forces directly induces a current, cf. the left hand side of Eq. (7). Here, the internal force density field arises from the interparticle interaction potential u(r^N) and is defined via the average

$${{\bf{F}}}_{{\rm{int}}}^{\alpha }({\bf{r}},t)=-\left\langle \sum _{i\in {{\mathcal{N}}}_{\alpha }}\delta ({\bf{r}}-{{\bf{r}}}_{i}){\nabla }_{i}u({{\bf{r}}}^{N})\right\rangle .$$

(8)

Furthermore the external force field ${{\bf{f}}}_{{\rm{ext}}}^{\alpha }({\bf{r}},t)$ in Eq. (7) acts on particles of species α, hence the external force acting on particle $i\in {{\mathcal{N}}}_{\alpha }$, as appearing in the Langevin equation, Eq. (1), is given by ${{\bf{f}}}_{{\rm{ext}},i}({{\bf{r}}}_{i},t)={{\bf{f}}}_{{\rm{ext}}}^{\alpha }({{\bf{r}}}_{i},t)$.

The internal one-body force density distributions can be split into adiabatic and superadiabatic one-body contributions^38,39,41, according to

$${{\bf{F}}}_{{\rm{int}}}^{\alpha }={{\bf{F}}}_{{\rm{ad}}}^{\alpha }+{{\bf{F}}}_{\sup }^{\alpha }.$$

(9)

Here, the adiabatic force density distribution ${{\bf{F}}}_{{\rm{ad}}}^{\alpha }({\bf{r}},t)$ is defined as the internal force density that occurs in a corresponding “adiabatic” equilibrium system that is defined by having the same partial density distributions as the nonequilibrium system, ${\rho }_{\alpha }^{{\rm{ad}}}({\bf{r}})={\rho }_{\alpha }({\bf{r}},t)$, ∀α. The adiabatic force density is also defined by Eq. (8), but with the crucial alteration that the average is now taken over an equilibrium probability distribution (i.e., that of the adiabatic system). Here, the external “adiabatic” one-body potentials ${V}_{{\rm{ad}}}^{\alpha }({\bf{r}})$ act on species α in order to stabilize the given partial density distributions, via (conservative) force fields $-\nabla {V}_{{\rm{ad}}}^{\alpha }({\bf{r}})$.

The superadiabatic force density distribution ${{\bf{F}}}_{\sup }^{\alpha }({\bf{r}},t)$ in Eq. (9) contains therefore all nonequilibrium effects which arise due to the presence of the flow in the system.

Dividing Eq. (7) by the partial density ρ_α and using Eqs. (5) and (9) we obtain the species-resolved force balance:

$$\gamma {{\bf{v}}}_{\alpha }=-{k}_{{\rm{B}}}T\nabla \mathrm{ln}{\rho }_{\alpha }+{{\bf{f}}}_{{\rm{ad}}}^{\alpha }+{{\bf{f}}}_{\sup }^{\alpha }+{{\bf{f}}}_{{\rm{ext}}}^{\alpha }.$$

(10)

The total density profile ρ(r, t) and the total current distribution J(r, t) can be obtained, respectively, by summation of Eqs. (3) and (4) over all species, i.e.,

$$\rho ({\bf{r}},t)=\sum _{\alpha }{\rho }_{\alpha }({\bf{r}},t),$$

(11)

$${\bf{J}}({\bf{r}},t)=\sum _{\alpha }{{\bf{J}}}_{\alpha }({\bf{r}},t).$$

(12)

In analogy to the partial velocity, Eq. (5), the mean total velocity is then obtained as

$${\bf{v}}({\bf{r}},t)=\frac{{\bf{J}}({\bf{r}},t)}{\rho ({\bf{r}},t)}.$$

(13)

Summing the species-resolved continuity equation, Eq. (6), over all species and using Eqs. (11) and (12) yields the total continuity relation as

$$\frac{\partial \rho }{\partial t}=-\nabla \cdot {\bf{J}}.$$

(14)

Analogously, summing the species-resolved force density balance, Eq. (7), over all species yields the total force density balance,

$$\gamma {\bf{J}}=-{k}_{{\rm{B}}}T\nabla \rho +{{\bf{F}}}_{{\rm{int}}}+{{\bf{F}}}_{{\rm{ext}}},$$

(15)

where the total internal force density field is ${{\bf{F}}}_{{\rm{int}}}=\sum _{\alpha }{{\bf{F}}}_{{\rm{int}}}^{\alpha }$, with the force density ${{\bf{F}}}_{{\rm{int}}}^{\alpha }$ acting on species α defined via Eq. (8) and the total external force density given by ${{\bf{F}}}_{{\rm{ext}}}=\sum _{\alpha }{\rho }_{\alpha }{{\bf{f}}}_{{\rm{ext}}}^{\alpha }$.

We next divide Eq. (15) by the total density profile ρ(r, t) and use the definition of the total velocity profile, Eq. (13), in order to obtain

$$\gamma {\bf{v}}=-{k}_{{\rm{B}}}T\nabla \mathrm{ln}\rho +{{\bf{f}}}_{{\rm{int}}}+{{\bf{f}}}_{{\rm{ext}}},$$

(16)

where we have defined the total internal and external force fields, f_int(r, t) and f_ext(r, t), respectively, as

$${{\bf{f}}}_{{\rm{int}}}={{\bf{F}}}_{{\rm{int}}}/\rho ,$$

(17)

$${{\bf{f}}}_{{\rm{ext}}}=\sum _{\alpha }\frac{{\rho }_{\alpha }}{\rho }{{\bf{f}}}_{{\rm{ext}}}^{\alpha }.$$

(18)

As we demonstrate below, f_int(r, t) is one crucial one-body field that enables us to rationalize the nonequilibrium behavior of driven mixtures.

We restrict ourselves in the following to binary mixtures, such that the species are labeled by α = 1, 2. We view the total force field f_int to act on both species in the same way, and also introduce a differential force density G_int(r, t) which drives the two species through and against each other. Hence using f_int and G_int we express the underlying species-resolved internal force density fields ${{\bf{F}}}_{{\rm{int}}}^{\alpha }$ for α = 1, 2 as

$${{\bf{F}}}_{{\rm{int}}}^{(1)}={\rho }_{1}{{\bf{f}}}_{{\rm{int}}}+{{\bf{G}}}_{{\rm{int}}},$$

(19)

$${{\bf{F}}}_{{\rm{int}}}^{(2)}={\rho }_{2}{{\bf{f}}}_{{\rm{int}}}-{{\bf{G}}}_{{\rm{int}}}.$$

(20)

We can invert this relationship and use Eq. (17) to obtain

$${{\bf{G}}}_{{\rm{int}}}=\frac{{\rho }_{2}}{\rho }{{\bf{F}}}_{{\rm{int}}}^{(1)}-\frac{{\rho }_{1}}{\rho }{{\bf{F}}}_{{\rm{int}}}^{(2)}.$$

(21)

Hence Eqs. (17) and (21) express a variable transformation from the species-resolved force densities ${{\bf{F}}}_{{\rm{int}}}^{(1)}({\bf{r}},t)$ and ${{\bf{F}}}_{{\rm{int}}}^{(2)}({\bf{r}},t)$ to the total force field f_int(r, t) and the differential force density G_int(r, t). As we demonstrate below, each of the new one-body fields describes a unique physical effect. Briefly, f_int contains effects that are also present in a one-component system, while G_int contains the genuine mixture contributions, such as drag and superdemixing.

The dynamics of the total density field are still governed by Eq. (16). To describe the mixture contributions, we introduce the density difference ρ_Δ and the current difference J_Δ, respectively, via

$${\rho }_{\Delta }={\rho }_{2}-{\rho }_{1},$$

(22)

$${{\bf{J}}}_{\Delta }={{\bf{J}}}_{2}-{{\bf{J}}}_{1},$$

(23)

where the continuity equation ∂ρ_Δ ∕ ∂t = −∇ ⋅ J_Δ is readily obtained from subtraction of the species-resolved continuity equations, Eq. (6). We define the differential external force density as

$${{\bf{F}}}_{{\rm{ext}}}^{\Delta }={{\bf{F}}}_{{\rm{ext}}}^{(2)}-{{\bf{F}}}_{{\rm{ext}}}^{(1)}\equiv {\rho }_{2}{{\bf{f}}}_{{\rm{ext}}}^{(2)}-{\rho }_{1}{{\bf{f}}}_{{\rm{ext}}}^{(1)}.$$

(24)

The corresponding equation of motion is obtained by building the difference of the species-resolved one-body force density balance, Eq. (7), for α = 1 and 2, which yields

$$\gamma {{\bf{J}}}_{\Delta }=-{k}_{{\rm{B}}}T\nabla {\rho }_{\Delta }+{{\bf{f}}}_{{\rm{int}}}{\rho }_{\Delta }+2{{\bf{G}}}_{{\rm{int}}}+{{\bf{F}}}_{{\rm{ext}}}^{\Delta },$$

(25)

where, as before, f_int(r, t) is given by Eq. (17) and G_int(r, t) is given by Eq. (21). The structure of Eq. (25) is crucial for the dynamics of driven mixtures. The differential current (left hand side) emerges from four different types of force density (right hand side): the first term is the ideal diffusive force density field due to gradients in the density difference. The second term generates a transport effect on ρ_Δ which is induced by the presence of f_int. The third term is a genuine mixture contribution that acts directly on the density difference; recall the definition of G_int, Eq. (21). The fourth term is due to the external influence, cf. Eq. (24).

We next split the fields into adiabatic contributions (i.e., those that can be understood on the basis of a corresponding equilibrium system with the same partial density profiles) and superadiabatic contributions of genuine nonequilibrium character, according to

$${{\bf{f}}}_{{\rm{int}}}={{\bf{f}}}_{{\rm{ad}}}+{{\bf{f}}}_{\sup },$$

(26)

$${{\bf{G}}}_{{\rm{int}}}={{\bf{G}}}_{{\rm{ad}}}+{{\bf{G}}}_{\sup }.$$

(27)

Due to linearity, the same variable transformation as before, Eqs. (17) and (21), relates the terms on the right hand side with the species-resolved force densities, Eq. (9). Hence for the adiabatic contributions:

$${{\bf{f}}}_{{\rm{ad}}}={{\bf{F}}}_{{\rm{ad}}}/\rho ,$$

(28)

$${{\bf{G}}}_{{\rm{ad}}}=\frac{{\rho }_{2}}{\rho }{{\bf{F}}}_{{\rm{ad}}}^{(1)}-\frac{{\rho }_{1}}{\rho }{{\bf{F}}}_{\sup }^{(2)},$$

(29)

where the total adiabatic force density is defined as ${{\bf{F}}}_{{\rm{ad}}}=\sum _{\alpha }{{\bf{F}}}_{{\rm{ad}}}^{\alpha }$. For the superadiabatic contributions:

$${{\bf{f}}}_{\sup }={{\bf{F}}}_{\sup }/\rho ,$$

(30)

$${{\bf{G}}}_{\sup }=\frac{{\rho }_{2}}{\rho }{{\bf{F}}}_{\sup }^{(1)}-\frac{{\rho }_{1}}{\rho }{{\bf{F}}}_{\sup }^{(2)},$$

(31)

where the total superadiabatic force density is defined as ${{\bf{F}}}_{\sup }=\sum _{\alpha }{{\bf{F}}}_{\sup }^{\alpha }$.

We summarize by inserting the adiabatic-superadiabatic splitting, Eqs. (26) and (27), into the velocity equation of motion, Eq. (16), and into the differential current, Eq. (25), which yields, respectively,

$$\gamma {\bf{v}}=-{k}_{{\rm{B}}}T\nabla \mathrm{ln}\rho +{{\bf{f}}}_{{\rm{ad}}}+{{\bf{f}}}_{\sup }+{{\bf{f}}}_{{\rm{ext}}},$$

(32)

$$\begin{array}{ll}\gamma {{\bf{J}}}_{\Delta }\,=\,-{k}_{{\rm{B}}}T\nabla {\rho }_{\Delta }+{\rho }_{\Delta }{{\bf{f}}}_{{\rm{ad}}}+{\rho }_{\Delta }{{\bf{f}}}_{\sup }\\ \,\,\,\,\,\,\,\,+\,2{{\bf{G}}}_{{\rm{ad}}}+2{{\bf{G}}}_{\sup }+{{\bf{F}}}_{{\rm{ext}}}^{\Delta }.\end{array}$$

(33)

Equations (32) and (33) form the basis of our subsequent classification of the different types of occurring physical effects.

Ideal mixture

In order to highlight the fundamental nonequilibrium effects, we particularize to an “ideal” mixture, in which the internal interactions do not depend on the type of particle. Formally, this implies that the value of the internal interaction potential u(r^N) is unchanged under permutation of the positions. For pair potentials (as we consider below) the inter-species pair potentials ϕ_αα(r) are identical to each other and are identical to the cross interaction potential ${\phi }_{\alpha \alpha ^{\prime} }(r)$ between particles of different species α and $\alpha ^{\prime}$. Hence, ${\phi }_{\alpha \alpha ^{\prime} }(r)=\phi (r)$, where ϕ(r) is a universal function. Given that all the internal interactions are the same, the adiabatic force field only has a nontrivial component, ${{\bf{f}}}_{{\rm{ad}}}^{\alpha }={{\bf{f}}}_{{\rm{ad}}}$, where f_ad is universal (independent of species). No differential component occurs in the adiabatic system, hence G_ad = 0, which simplifies the dynamics, cf. Eq. (33). Importantly, this symmetry does not apply to the external force field. Hence, in general ${{\bf{f}}}_{{\rm{ext}}}^{\alpha }({\bf{r}},t)\,\ne\,{{\bf{f}}}_{{\rm{ext}}}^{\alpha ^{\prime} }({\bf{r}},t)$ for $\alpha\,\ne\,\alpha^{\prime}$, which imprints differences into the kinematic one-body fields ρ_α and J_α during the time evolution. Such driving is far from trivial. One might picture differently colored particles being driven through each other.

Next, we consider a special simplifying situation consisting of steady states characterized by a single direction of the flow along which all vector fields point, i.e., v, v₁, v₂, J, J₁, J₂ are all collinear. In the simulation results shown below this is the $\hat{{\bf{y}}}$-direction (the flow is in the “vertical” direction). All gradients in the system also share a common direction, which is orthogonal to the flow. In the simulations, this is the $\hat{{\bf{x}}}$-direction (the system is “horizontally” inhomogeneous). Such “perpendicular flow” geometry forms a class of common nonequilibrium situation encompassing e.g., simple shear flow, steady flow between parallel plates etc.

We split the superadiabatic forces into viscous (subscript “visc”), drag, and structural (subscript “struc”) contributions according to

$${{\bf{f}}}_{\sup }={{\bf{f}}}_{{\rm{visc}}}+{{\bf{f}}}_{{\rm{struc}}},$$

(34)

$${{\bf{G}}}_{\sup }={{\bf{G}}}_{{\rm{drag}}}+{{\bf{G}}}_{{\rm{struc}}},$$

(35)

where, as we show, f_visc is a viscous force field that arises from total shear motion, f_struc is a structural force field that sustains gradients of the total density distribution, G_drag is a superadiabatic force density distribution that describes the friction effect that occurs due to counterflow (i.e., when v_Δ ≠ 0), and G_struc is the superadiabatic demixing force density. We demonstrate the validity of this interpretation below. On a formal level, the splitting in Eqs. (34) and (35) is uniquely defined by f_visc and G_drag acting collinear with the flow, and f_struc and G_struc acting perpendicular to it.

We can now separate the equations of motion according to the direction of the forces. In the flow direction, using Eqs. (34) and (35), we obtain from Eqs. (32) and (33), respectively,

$$\gamma {\bf{v}}={{\bf{f}}}_{{\rm{visc}}}+{{\bf{f}}}_{{\rm{ext}}},$$

(36)

$$\gamma {{\bf{J}}}_{\Delta }={\rho }_{\Delta }{{\bf{f}}}_{{\rm{visc}}}+2{{\bf{G}}}_{{\rm{drag}}}+{{\bf{F}}}_{{\rm{ext}}}^{\Delta },$$

(37)

which indicates that the total (scaled) flow (γv) arises from a competition of the external driving (f_ext) and the viscous forces (f_visc). The differential current (γJ_Δ) is counteracted by a transport effect that the viscous forces exert on the density difference (ρ_Δf_vis), and the drag that any counterflow of the two species induces (2G_drag), and it is driven by the external differential force density (${{\bf{F}}}_{{\rm{ext}}}^{\Delta }$).

In the gradient direction, we also use Eqs. (32) and (33), respectively, to obtain

$$0=-{k}_{{\rm{B}}}T\nabla \mathrm{ln}\,\rho +{{\bf{f}}}_{{\rm{ad}}}+{{\bf{f}}}_{{\rm{struc}}}$$

(38)

$$0=-{k}_{{\rm{B}}}T\nabla {\rho }_{\Delta }+{\rho }_{\Delta }({{\bf{f}}}_{{\rm{ad}}}+{{\bf{f}}}_{{\rm{struc}}})+2{{\bf{G}}}_{{\rm{struc}}}.$$

(39)

Here, remarkably, the external forces do not appear explicitly. Equation (38) constitutes a force balance relationship, where the sum of the diffusive force, the adiabatic force, and the structural force cancel, such that no motion occurs. Equation (39) is crucial for understanding the physics of driven mixtures: The sum of the diffusive force density difference,  −k_BT ∇ ρ_Δ, and the effect of total adiabatic and structural force fields on the density difference, ρ_Δ(f_ad + f_struc), is balanced by the structural differential force density, 2G_struc. This implies that G_struc ≠ 0 can induce a density contrast ρ_Δ ≠ 0 in the system.

We can now revert to the species-labeled description and start from the general relationship

$${{\bf{v}}}_{\alpha }\,=\,\frac{\rho {\bf{v}}\,\pm\,{{\bf{J}}}_{\Delta }}{\rho\,\pm\, {\rho }_{\Delta }},$$

(40)

where the upper sign (+) applies to species 1, and the lower sign (−) applies to species 2. Application to the perpendicular geometry results in the following balance relations in the flow direction and in the gradient direction

$$\gamma {{\bf{v}}}_{\alpha }={{\bf{f}}}_{{\rm{visc}}}\pm \frac{{{\bf{G}}}_{{\rm{drag}}}}{{\rho }_{\alpha }}+{{\bf{f}}}_{{\rm{ext}}}^{\alpha }\quad ({\rm{flow}}),$$

(41)

$$0={{\bf{f}}}_{{\rm{ad}}}+{{\bf{f}}}_{{\rm{struc}}}\pm \frac{{{\bf{G}}}_{{\rm{struc}}}}{{\rho }_{\alpha }}-{k}_{{\rm{B}}}T\nabla \mathrm{ln}{\rho }_{\alpha }\quad ({\rm{grad}}),$$

(42)

where again  ±  corresponds to α = 1, 2. As we have four equations for the four unknown superadiabatic fields, we can solve for these, with the results:

$${{\bf{f}}}_{{\rm{visc}}}=\gamma {\bf{v}}-{{\bf{f}}}_{{\rm{ext}}},$$

(43)

$${{\bf{G}}}_{{\rm{drag}}}=\frac{\gamma }{2}({{\bf{J}}}_{\Delta }-{\rho }_{\Delta }{\bf{v}}+{\rho }_{\Delta }{{\bf{f}}}_{{\rm{ext}}}-{{\bf{F}}}_{{\rm{ext}}}^{\Delta })$$

(44)

$${{\bf{f}}}_{{\rm{struc}}}={k}_{{\rm{B}}}T\nabla \mathrm{ln}\rho -{{\bf{f}}}_{{\rm{ad}}},$$

(45)

$${{\bf{G}}}_{{\rm{struc}}}=\frac{{k}_{{\rm{B}}}T}{\rho }\left({\rho }_{1}\nabla {\rho }_{2}-{\rho }_{2}\nabla {\rho }_{1}\right),$$

(46)

where f_ext is the total external force field given by the linear combination, Eq. (18), and ${{\bf{F}}}_{{\rm{ext}}}^{\Delta }$ is the differential external force density defined via Eq. (24); both fields depend on the partial density profiles. As all quantities on the right hand side of Eqs. (43–46) are accessible in simulations, we have formulated a means to obtain a complete splitting of the nonequilibrium force contributions in the system.

The total species-resolved superadiabatic force field is then given by

$${{\bf{f}}}_{\sup }^{\alpha }={{\bf{f}}}_{{\rm{visc}}}\pm \frac{{{\bf{G}}}_{{\rm{drag}}}}{{\rho }_{\alpha }}+{{\bf{f}}}_{{\rm{struc}}}\pm \frac{{{\bf{G}}}_{{\rm{struc}}}}{{\rho }_{\alpha }},$$

(47)

where the first and the second term on the right hand side act in flow direction, and the third and the fourth term act in the gradient direction; again  ±  corresponds to α = 1, 2.

At the end of the results section we also analyze briefly a second case in “planar” geometry, in which there is only a single direction in the system ($\hat{{\bf{y}}}$) along which the flow occurs and also all gradients point. The system is transitionally invariant in $\hat{{\bf{x}}}$. As the flow occurs in the direction of the gradients, in general, such situations will be nonstationary, i.e., time-dependent.

Power functional approximation

Power functional theory³⁸ (PFT) is an exact variational approach to nonequilibrium dynamics in overdamped Brownian systems. PFT reduces to density functional theory in equilibrium. Within PFT, the superadiabatic force is obtained as a functional derivative with respect to the current,

$${{\bf{f}}}_{\sup }^{\alpha }=-\frac{\delta {P}_{t}^{{\rm{exc}}}[\{{\rho }_{\alpha ^{\prime} },{{\bf{J}}}_{\alpha ^{\prime} }\}]}{\delta {{\bf{J}}}_{\alpha }}.$$

(48)

Here ${P}_{t}^{{\rm{exc}}}$ is the superadiabatic excess power functional³⁸, which is a functional of the density and current profiles of both species {ρ_α, J_α}. The superadiabatic force is therefore also a functional of both types of one-body fields: density profiles and current profiles.

Approximations for ${P}_{t}^{{\rm{exc}}}$ in monocomponent systems based on a series expansion of gradients of the velocity field have been proposed in refs. ^42,43. The resulting superadiabatic forces describe viscous⁴² and structural effects⁴³. In a binary mixture, besides the velocity gradients, the velocity difference between both species v_Δ = v₂ − v₁ is a further key ingredient to construct approximated power functionals⁴⁴.

The following approximation for ${P}_{t}^{{\rm{exc}}}$ reproduces semiquantitatively all types of superadiabatic forces in the system considered here:

$$ {P}_{t}^{{\rm{exc}}}=\frac{{C}_{{\rm{visc}}}}{2}\int {\mathrm{d}}{\bf{r}}{\rho }^{2}{\left(\nabla \times {\bf{v}}\right)}^{2}$$

(49)

$$ + \frac{{C}_{{\rm{drag}}}}{2}\int {\mathrm{d}}{\bf{r}}{\rho }_{1}{\rho }_{2}{{\bf{v}}}_{\Delta }^{2}$$

(50)

$$-{C}_{{\rm{struc}}}\sum _{\alpha }\int {\mathrm{d}}{\bf{r}}{\rho }^{2}{{\bf{v}}}_{\Delta }^{2}\nabla \cdot {{\bf{J}}}_{\alpha }$$

(51)

$$\ \ \ \ \ \ \ \ \ \ \ \ +\ {D}_{{\rm{struc}}}\sum _{\alpha }\ \int {\mathrm{d}}{\bf{r}}{\rho }^{2}\left[(\nabla \ \times \ {{\bf{v}}}_{\alpha })\ \times \ {{\bf{v}}}_{\Delta }\right]^{\prime} \ \cdot \ {{\bf{v}}}_{\Delta },$$

(52)

where ρ² indicates the square of the total density profile, as defined Eq. (11). The functional consists of four different terms, each one accounting for one type of superadiabatic force. The spatial dependence is local in density and semi-local (via the gradient) in velocity. The functional form in Eqs. (49–52) is based on expressions that contain low powers of the density and velocity profiles. The expressions in Eqs. (49) and (51) are of velocity gradient form⁴², with the former describing viscous⁴² and the latter describing structural effects⁴³. Furthermore, the differential velocity v_Δ, which captures the degree of interflow of the two components, appears quadratically in Eq. (50) as is appropriate to describe drag effects⁴⁴. It couples in Eq. (51) to the divergence of the current, which creates a structural force. Finally, Eq. (52) generates the superdemixing force, as we will demonstrate below.

The prefactors of each term in Eqs. (49–52) are positive constants that represent the amplitude of the respective superadiabatic force contribution. We use these coefficients as fitting parameters (Numerical values are given below in the caption of the respective figure). In reality, these prefactors depend on microscopic features of the model, as the interparticle potential and are, in general, functionals of the density. The prime in Eq. (52) indicates the following time integral

$$\left[(\nabla \times {{\bf{v}}}_{\alpha })\times {{\bf{v}}}_{\Delta }\right]^{\prime} \ =\ {\int }_{0}^{t}{\mathrm{d}}t^{\prime} K(t-t^{\prime} )(\nabla \times {\bf{v}}^{\prime} )\times {\bf{v}}^{\prime} ,$$

(53)

with K(t) being a temporal convolution kernel with normalization ∫dtK(t) = 1 and the primes on the right hand side indicate dependence on $t^{\prime}$. We do not investigate the specific form of K(t) in the present work further, as we are primarily concerned with steady states. The effective role that the kernel plays in steady state is that the functional derivative in Eq. (48) does not act on the primed terms.

The viscosity term. Equation (49), is the lowest order term in a power series expansion in the velocity gradient⁴² and accounts for the Stokes-like viscous force. A similar term containing the square of the divergence of the average velocity, instead of the square of the curl, does not produce any force in the current setup, since the velocity profile is free of divergence. The superadiabatic viscous force that results from differentiating Eq. (49) is species-independent and given by

$${{\bf{f}}}_{{\rm{visc}}}=\frac{{C}_{{\rm{visc}}}}{\rho }\left(\nabla {\rho }^{2}\cdot \nabla {\bf{v}}-\nabla \cdot {\rho }^{2}{(\nabla {\bf{v}})}^{\top }\right)$$

(54)

$$=\frac{{C}_{{\rm{visc}}}}{\rho }{\nabla }_{x}({\rho }^{2}{\nabla }_{x}{v}_{y}(x))\hat{{\bf{y}}},$$

(55)

where the leftmost derivatives act on each entire expression, and the second equation, Eq. (55), has been particularized for our current setup ${\bf{v}}({\bf{r}})={v}_{y}(x)\hat{{\bf{y}}}$ and constitutes a force in flow direction.

The second term, Eq. (50), describes the drag of particles of one species due to the flow of particles of the other species⁴⁴. The resulting superadiabatic force is species-dependent. The corresponding force density is:

$${{\bf{G}}}_{{\rm{drag}}}={C}_{{\rm{drag}}}{\rho }_{1}{\rho }_{2}{{\bf{v}}}_{\Delta }.$$

(56)

The structural superadiabatic force generated by the third term, Eq. (51), is species-independent,

$${{\bf{f}}}_{{\rm{struc}}}=-{C}_{{\rm{struc}}}\nabla {\left(\rho {{\bf{v}}}_{\Delta }\right)}^{2}.$$

(57)

As we will see, this force balances the sum of the adiabatic force and the diffusive force of the total density gradient.

The last term, Eq. (52), is also structural and the only one generating the superdemixing. The force density obtained via functional differentiation of Eq. (52) is

$${{\bf{G}}}_{{\rm{struc}}}=-{D}_{{\rm{struc}}}{\rho }^{2}\sum _{\alpha }\left[(\nabla \times {{\bf{v}}}_{\alpha })\times {{\bf{v}}}_{\Delta }\right]^{\prime} ,$$

(58)

where only the functional derivative with respect to the unprimed velocity difference in Eq. (52) produces a force due to the structure of Eq. (53). Before presenting results, we first lay out the specifics of the microscopic model that we consider.

Brownian dynamics simulations

In the following we explicitly consider an equimolar two-dimensional binary mixture with a total of N = 1066 particles (${{\mathcal{N}}}_{1}={{\mathcal{N}}}_{2}=532$) interacting via the same quasi-hard core pair potential ϕ(r) = ϵσ/r³⁶. Here, σ denotes the characteristic particle size, which sets our length scale, and ϵ sets the energy scale. The particles are subject to a homogeneous gravitational field in $\hat{{\bf{y}}}$-direction,

$${{\bf{f}}}_{{\rm{ext}}}^{\alpha }=-{m}_{\alpha }g\hat{{\bf{y}}},$$

(59)

where g is the acceleration due to gravity, $\hat{{\bf{y}}}$ is a unit vector along the vertical $\hat{{\bf{y}}}$-direction, and m_α is the buoyant mass of species α. The only difference between the two species is that they have opposite buoyant masses, i.e., m₁ = −m₂. Although we make use of a gravitational field to illustrate the process of lane formation, other external fields such as magnetic and electric fields can also be used.

We simulate the system using Brownian dynamics simulations performed in a square simulation box of length L/σ = 35 with periodic boundary conditions in both directions. The total bulk density is hence ρ_bσ² = Nσ²∕L² = 0.87. The temperature is fixed to k_BT/ϵ = 0.5. The particle trajectories are integrated in time via the Euler algorithm with a discrete time step Δt/τ = 3.0 × 10⁻⁵ in units of the reduced time τ = σ²γ/ϵ. The simulations of the laned state are initialized with the particles randomly located in a homogeneous state. Once the steady state is reached the system is sampled for a total time of  ~10⁶ τ.

Laned state

The magnitude of the external force that drives the two species against each other, m_αg, highly influences the behavior of the system. Without gravity, particles of both species form an effectively homogeneous one-component system in equilibrium, in which the only difference between the species is an arbitrary label. At low values of gravity, the steady state remains homogeneously mixed, and the particles are driven slowly through each other, see the characteristic snapshot in Fig. 1a for a value of the external field m₁g = 0.05 ϵ/σ. No forces act in the $\hat{{\bf{x}}}$-direction. The internal interactions lead to weak superadiabatic drag forces in the (flow) $\hat{{\bf{y}}}$-direction that partially counteract the external force.

For sufficiently high external driving, such as e.g., m₁g = 1 ϵ/σ, the two species segregate into two lanes, see Fig. 1b. Each lane is characterized by the concentration of one species being significantly higher than that of the other (minority) species. In each lane there is a current of particles parallel to the external force corresponding to the dominant species. As shown above, there occur superadiabatic drag and viscous forces opposing the external forces, see Eqs. (36) and (37). In addition, in the (gradient) $\hat{{\bf{x}}}$-direction there is a balance of forces, cf. Eqs. (38) and (39). The sum of the ideal diffusion and the adiabatic force tries to mix the two species. These forces are, however, canceled by structural superadiabatic forces that keep the system in the demixed laned steady state. A brief animation of the particle flow in the homogeneous and in the laned stated is shown in Supplementary Movie 1.

The density profiles of the laned state are shown in Fig. 2a. There, the right lane (x > 0) is dominated by species 1 that possesses a positive buoyant mass, i.e., the external force is directed downwards. In the left lane (x < 0) the majority of particles belongs to species 2 with negative buoyant mass and therefore an upwards directed external force. Although here species 1 predominantly occupies the right lane, the state is degenerated in the sense that a swap of the two lanes can be compensated by a shift by L/2 in the $\hat{{\bf{x}}}$-direction (due to periodic boundary conditions). The total density profile, Fig. 2a, increases toward the centers of the lanes and decreases at the interfaces, although it is quite homogeneous as compared to the single species density profiles.

Fig. 2: Lane formation.

Kinematic fields: (a) total density profile ρ (black dotted line) and partial density profiles ρ₁ (orange solid line) and ρ₂ (blue dashed line), (b) flow component of the total velocity profile v^y (black dotted line), see Eq. (13), and of the partial velocity profiles ${v}_{1}^{y}$ (orange solid line) and ${v}_{2}^{y}$ (blue dashed line). Species-resolved forces (c) and force densities (d) acting in flow direction, Eq. (10): external ${{\bf{f}}}_{{\rm{ext}}}^{\alpha }$ (dashed lines), and superadiabatic ${{\bf{f}}}_{\sup }^{\alpha }$ (solid lines), which is the sum of viscous and drag, cf. Eq. (47). Species-resolved forces (e) and force densities (f) acting in gradient direction: sum of adiabatic and diffusive forces (dotted lines) and superadiabatic (solid lines), cf. Eq. (42). In all plots the line color indicates the species: orange (blue) for species 1 (2). The arrows indicate the direction of the force or force density. Black circles indicate points where the force or force density vanishes. All plots are for driving m₁g = 1.0ϵ∕σ and all values are presented as functions of the x-coordinate.

In Fig. 2b we show the velocity profiles. Only the flow component of the velocity profiles ${v}_{\alpha }^{y}$ does not vanish in steady state. The total velocity profile clearly shows the motion of particles in opposite directions in the two lanes, irrespective of the species. The total velocity is directed downwards in the right lane and upwards in the left lane, while being zero at the interfaces, as is consistent with the symmetry of the system.

Forces acting in flow direction

We focus first on the species resolved force balance equation, Eq. (10). The flow components of the forces and of the force densities are presented in Fig. 2c, d, respectively. In what follows we focus the discussion on species 1. Due to the symmetry of the system the behavior of species 2 follows directly. Only two types of forces act in flow direction, the constant external force field ${{\bf{f}}}_{{\rm{ext}}}^{\alpha }$, Eq. (59), and the superadiabatic force field due to interparticle interactions ${{\bf{f}}}_{\sup }^{\alpha }$. Adiabatic (${{\bf{f}}}_{{\rm{ad}}}^{\alpha }$) and diffusive ($-{k}_{{\rm{B}}}T\nabla \mathrm{ln}{\rho }_{\alpha }$) forces appear only when gradients of the density profile are present. As the system is homogeneous in the flow direction, such contributions do not occur in this direction. Therefore, the internal forces in flow direction are purely superadiabatic.

For those particles located in the lane for which its species is in the majority (e.g., species 1 in the right lane in Fig. 2b), the velocity and the external force are parallel, while the internal force opposes both. For those particles in the lane for which its species is in the minority (e.g., species 2 in the right lane in Fig. 2b) the velocity and the external force are antiparallel. This is caused by a strong internal force, mostly created by particles of the opposite species which dominates the lane, that drags the particles of the minority species. The force densities in flow direction are shown in Fig. 2d. Force and force densities in gradient direction are depicted in Fig. 2e, f, respectively. The force densities are the actual contributions to the current of particles (up to a multiplicative constant given by the inverse friction coefficient) and are therefore a macroscopically accessible quantity.

Next, we split the species-resolved superadiabatic forces, Eq. (47), in flow direction in its species-dependent and species-independent contributions, Eq. (41), which are the viscous force f_visc, Eq. (43), and the drag  ±G_drag, Eq. (44), respectively. Both are of the superadiabatic kind.

The forces are shown in Fig. 3a, b, and the force densities in Fig. 3c. The species-dependent term  ±G_drag is the dominant part and describes the drag of one species due to motion of the other species. For example, for species 1 this force density is high in the lane of particles of species 2 (left lane) and points upwards. The species-independent force f_visc points in the direction opposite to the total velocity v, as corresponds to a viscous effect. The amplitude of the viscous force is relatively small as compared to that of the drag force. However, the maxima (absolute value) of f_visc occur at the peaks of the partial densities. As a result, the viscous force density ρ_αf_visc introduces a relevant correction to the dominant species-dependent force density, see Fig. 3c.

Fig. 3: Superadiabatic forces in flow directions as a function of the x − coordinate.

a Species-independent superadiabatic force f_visc (dark solid green), species-dependent superadiabatic force acting on species 1 (solid orange), G_drag∕ρ₁, and 2 (solid blue),  −G_drag∕ρ₂. b Enlarged view of the viscous force. c Species-dependent force density G_drag acting on species 1 (solid orange), and species-independent force density acting on species 1 (orange-dark green), f_viscρ₁. The arrows (colored according to the corresponding data set) indicate the direction of the force. The symbols show the theoretical predictions for the viscous force (a), (b) force and the drag force density (c). We have set C_visc = 2ϵτσ² and C_drag = 8.5ϵτ in Eqs. (55) and (56), respectively.

A comparison between the Stokes-like viscous force predicted by our power functional approximation, Eq. (55), and the viscous force measured in simulations is shown in Fig. 3b. The predicted force has the correct sign everywhere (opposite of the total flow). However, due to the spatially local approximation taken in Eq. (55) and the strong driving conditions, the agreement is not perfect. The drag force predicted by PFT, Eq. (56), also shows a good agreement with the simulation data, see Fig. 3c.

Forces acting in gradient direction

In the gradient direction $\hat{{\bf{x}}}$, which is perpendicular to the flow, no external force is applied and the system is in steady state with no flow along this direction, cf. Eq. (42). The density modulation (formation of lanes), Fig. 2a, is solely created by the flow along $\hat{{\bf{y}}}$, Fig. 2b.

The species-resolved net force balance, Eq. (10), vanishes in the $\hat{{\bf{x}}}$-direction, and both adiabatic ${{\bf{f}}}_{{\rm{ad}}}^{\alpha }$ and superadiabatic ${{\bf{f}}}_{\sup }^{\alpha }$ forces contribute to the internal force. For species 1 the sum of adiabatic and diffusive, $-{k}_{{\rm{B}}}T\nabla \mathrm{ln}{\rho }_{\alpha }$, forces attempts to move particles in the center of the simulation box to the left and particles at the borders of the box to the right, see Fig. 2e, f. Hence, as could be expected, the sum of the diffusive and the adiabatic forces pushes particles outside of the majority lane and therefore tries to mix the two species by eliminating density gradients. The superadiabatic force in gradient direction balances the diffusion and the adiabatic forces by pushing particles into their majority lane and out of their minority lane, Fig. 2e, f. The superadiabatic force is therefore the only force that creates and sustains the demixing and segregation into lanes in steady state. Superadiabatic forces are therefore crucial for the description of lane formation, as expected given the nonequilibrium nature of the phenomenon.

The adiabatic forces are obtained by sampling the internal forces in the adiabatic equilibrium system. The adiabatic equilibrium system is also a symmetric mixture with the same quasi-hard interparticle interactions regardless of the species. The only difference between the real nonequilibrium and the adiabatic equilibrium systems is that the external nonconservative driving force is switched off in the adiabatic equilibrium system and replaced by conservative external forces generated by two potentials ${V}_{{\rm{ad}}}^{\alpha }$ that reproduce the same density profiles as in nonequilibrium. The adiabatic potentials are obtained following the procedure described in “Methods” section, where also results are shown for the adiabatic potentials and the density profiles in the adiabatic system.

In the adiabatic system, the external potential counteracts both the ideal diffusive force $-{k}_{{\rm{B}}}T\nabla \mathrm{ln}{\rho }_{\alpha }$ and the adiabatic (internal) force ${{\bf{f}}}_{{\rm{ad}}}^{\alpha }$. The diffusive force and the adiabatic force have two very distinct effects, which are schematically illustrated in Fig. 4a. Quantitative plots of both force and force density profiles for species 1 are shown in Fig. 4b, c, respectively. Force and force density profiles for species 2 are depicted in Fig. 4d, e. The diffusive force is the only force that attempts to mix both species. In contrast, the adiabatic force field is the same for all particles, independently of the species, due to the symmetry of the internal interactions. That is ${{\bf{f}}}_{{\rm{ad}}}^{\alpha }={{\bf{f}}}_{{\rm{ad}}}$, as discussed before. The adiabatic force tries to eliminate the density gradients of the total density profile ρ = ρ₁ + ρ₂ by pushing particles of both species toward the interfaces between the two lanes.

Fig. 4: Adiabatic forces.

a Partial snapshots of the laned state at gravity m₁g = 1.0ϵ∕σ. The arrows indicate the direction of the forces acting on species 1 and 2. The diffusive force field $-{k}_{{\rm{B}}}T\nabla \mathrm{ln}{\rho }_{\alpha }$ (violet arrows) tries to mix both species. The adiabatic force field (green) is the same for both species. Force (b) and force density (c) profiles acting on species 1 in $\hat{{\bf{x}}}$-direction as a function of x: total force: adiabatic plus diffusion (orange dotted), diffusive force (violet dashed), and adiabatic force (green solid). Panels d and e are the same as panels b and c, respectively, but for species 2: total force (blue dotted), diffusion (violet dashed), and adiabatic (green solid).

We next split the species-resolved superadiabatic forces, Eq. (47), into the species-independent f_struc, Eq. (45), and the species-dependent  ±G_struc, Eq. (46), contributions in gradient direction. Both contributions are of the structural type in the sense that they sustain the density gradients and are non-dissipative. Fig. 5a illustrates schematically the different effects of f_struc and G_struc. The actual profiles are displayed in Fig. 5b, c.

Fig. 5: Superadiabatic forces in gradient direction.

a Partial snapshots of the laned state (m₁g = 1.0ϵ∕σ). The arrows indicate the direction of the species-independent f_struc and dependent  ± G_struc superadiabatic forces in gradient direction. The brightness of the particles in the upper snapshot indicates the total density modulation (brighter regions are less dense). Species-independent superadiabatic force f_struc (b) and species-dependent force density G_struc acting on species 1 (c). Lines are simulation data and symbols denote the power functional approximation given in Eqs. (45) and (58) with parameters C_struc = 10.5ϵτ²σ², and D_struc = 0.73ϵτ². d Negative of the species-independent superadiabatic force−f_struc (dark green dashed) and the adiabatic excess force f_ad (bright green solid). e Negative of the thermal diffusion force associated to the total density (blue) and sum of the species-independent superadiabatic force and the adiabatic force (black dotted).

The species-independent superadiabatic force in gradient direction, f_struc, pushes particles from the interfaces toward the centers of the lanes, creating the total densty modulation, see Fig. 5b. f_struc is almost the opposite of the adiabatic force f_ad, as shown in Fig. 5d where we plot the adiabatic force, and  −f_struc. These two terms, however, do not exactly cancel each other as the sum of f_struc and the adiabatic force shows a clear signal (see Fig. 5e) that cancels exactly the diffusive force orginated by the gradient of the total density profile, cf. Eq. (45). This diffusive term also describes a species-independent effect.

The species-dependent force density +G_struc transports particles of species 1 away from the left lane into the right lane, see Fig. 5c. The term  −G_struc acts on particles of species 2 in the opposite way. This superadiabatic contribution is precisely the opposite of the contribution due to the local demixing, compare Fig. 5c and the diffusive force density shown in Fig. 4c. Therefore, the species-dependent term G_struc sustains the demixing of the two species in two lanes in steady state.

The sum of both superadiabatic structural terms in gradient direction cancels the diffusive and the adiabatic contributions, cf. Eq. (42). Our power functional approximation, Eqs. (57) and (58), is in very good agreement with the simulation data for both f_struc and G_struc as shown in Fig. 5b, c, respectively.

Jammed state

We conclude the results section with a brief discussion of the jammed state, see Fig. 6a. The jammed state can be found at intermediate values of gravity, and in it flow and gradient directions coincide. This state occurs if the two species block each other during the diffusive process induced by the external field. The state is homogeneous in $\hat{{\bf{x}}}$-direction, and it is characterized by a significant change of the total density in $\hat{{\bf{y}}}$-direction. The opposite external driving of the two species leads to the formation of a compressed region with similar concentrations of both species, high total density, and incipient crystalline order. This region is surrounded by two regions rich in particles of only one species that press against the compressed area. In the jammed state, strong superadiabatic forces act against the external fields. In addition, the adiabatic and diffusive forces try to smooth the total density profile. In contrast to the laned state, the demixing in the jammed state occurs along the direction of the external driving force. In our simulations the jammed state is metastable, but it persists for very long times, which allows for systematic study of the forces.

Fig. 6: Jammed state formed at intermediate values of the external driving.

a Characteristic snapshot. The orange (blue) arrows indicate the forces acting on species 1 (2). b Total density profile ρ (black dotted line), and partial density profiles ρ₁ (solid orange) and ρ₂ (blue dashed). c Total velocity profile v^y (black dotted), and partial velocity profiles ${v}_{1}^{y}$ (orange solid) and ${v}_{2}^{y}$ (blue dashed). d Current profiles: total current profile j^y (black dotted), and partial current profiles ${j}_{1}^{y}$ (orange solid) and ${j}_{2}^{y}$ (blue dashed). e Forces acting in $\hat{{\bf{y}}}$-direction: external (dashed), diffusive plus adiabatic (dotted), and superadiabatic (solid). The color indicates the species: orange (blue) for species 1 (2). f Adiabatic and diffusive forces acting in $\hat{{\bf{y}}}$-direction acting on species 1: Sum of both forces (orange dotted), adiabatic (green solid), and ideal diffusion (purple dashed). Data obtained for a value of gravity m₁g  = 0.7ϵ∕σ. All quantities are presented as a function of the y-coordinate.

To characterize the jammed state we set the gravity to m₁g = 0.7ϵ ∕ σ and split the time evolution of the system in time windows of 20τ. Only those time windows in which the difference between the maximum and the minimum of the total density is higher than 0.2σ² are considered to be in the jammed state. We have found by visual inspection of microstates that this is a sufficiently accurate criterion for identifying states that one would classify as being transiently jammed. For flowing states and for less well-characterized states, the density distribution is typically more homogeneous, such that the chosen criterion removes these states from the sampling. Using several simulations, we average over these relevant windows for a total sampling time of 10⁵ τ. The density profiles, Fig. 6b, show the demixing along the $\hat{{\bf{y}}}$-axis. The particles of species 1 (2) accumulate at the top (bottom) of the simulation box. In the center of the box there is a region with similar concentrations of both species. The total density shows a small modulation with a global maximum located in the middle of the box (in each time window we translate the data such that the maximum density occurs in the the middle of the box). The velocity profiles are shown in Fig. 6c. Only the $\hat{{\bf{y}}}$-component of the velocity is non-zero. The total velocity is almost negligible across the entire systemf. The velocity of species 1 almost vanishes everywhere except at the bottom of the simulation box, i.e., in the region rich in particles of species 2. Here, the velocity of species 1 is negative. Hence, particles of species 1 move toward the top of the simulation box (recall that we use periodic boundary conditions) where the species 1 dominates. Exactly the opposite behavior is observed in species 2, as expected due to symmetry considerations.

The total current vanishes, see Fig. 6d, and the partial currents are small and indicate that the system is still undergoing compression. Given that the partial currents are not constant, it is clear the system is not in steady state. Nevertheless, the currents are very small and rather homogeneous, which justifies the assumption of a quasi-steady state. For other systems the jammed state has been observed to be steady. This is the case of a 2D hard disk system with very low opposite driving at zero temperature⁴⁵, and also on a square lattice model²¹.

The species-resolved forces in $\hat{{\bf{y}}}$-direction split into diffusion, adiabatic, superadiabatic, and external, cf. Eq. (10). Plots of each contribution are shown in Fig. 6e. We compute the adiabatic force via the adiabatic potential, as described in “Methods” section.

The (constant) external force for both species point in opposite directions. The sum of the thermal diffusion and the adiabatic force is for both species positive at the top and mostly negative at the bottom of the simulation box. Hence, this force field tries to smooth the density modulations by transporting particles away from the center toward the low density regions. The superadiabatic force for species 1 is everywhere directed upwards against the external force and has a maximum at the bottom of the sample, the region dominated by species 2.

The thermal diffusion and the adiabatic force acting on species 1 are shown in Fig. 6f. The adiabatic force is positive at the top of the box and negative at the bottom, moving particles away from the center in order to smooth the total density profile. Similarly, the ideal diffusion force tries to smooth the partial density profile of species 1.

Recall that in the laned state, the diffusive force is the only force counteracting the demixing, while the adiabatic force has a phenomenologically irrelevant effect. In the jammed state, however, the adiabatic force has a prominent effect, while the diffusion term is only a small contribution to the forces present in the adiabatic state. The flow in the jammed state is curl-free instead of divergence-free (which is the case in in the laned state). Hence, the leading terms in the expansion of the excess free power for the viscous and the species-independent structural forces will differ from those in Eqs. (49) and (52), respectively.

Discussion

By splitting the nonequilibrium superadiabatic forces of a binary colloidal mixture into species-dependent and independent contributions, we have identified a structural force contribution that induces demixing of two different species in nonequilibrium. This superadiabatic force is responsible for the formation of lanes in oppositely driven binary colloidal systems. The force counteracts the ideal thermal diffusive force that arises as a consequence of the mixing entropy generated by a gradient in composition between the two lanes.

The entropy-driven demixing that occurs in equilibrium mixtures of hard bodies, see e.g., ref. ⁴⁶, is a well-known phenomenon. There, the interparticle interactions generates an internal adiabatic force that counteracts the ideal entropy of mixing (diffusive force) and hence sustains the density gradients. The mechanism behind the demixing in lane formation is completely different. Here, the adiabatic force only leads to a minor modulation of the total density and does not attempt to demix the species. Although we have analyzed ideal symmetric mixtures, we expect this to be true also in asymmetric mixtures that could be considered in future studies.

Lane formation is a purely superadiabatic effect. Consequently a direct theoretical description of laning via dynamical density functional theory^47,48,49, in which only adiabatic forces are considered, is not possible and phenomenological ad-hoc corrections must be included¹⁴. For the related, but different case of colloidal suspensions in confinement, Aerov and Krüger argued⁵⁰ that dynamical density functional is insufficient to describe the effects of shear.

We have presented here a different approach based on power functional theory³⁸, an exact variational approach for nonequilibrium situations. Power functional theory incorporates superadiabatic effects via a functional of both density and velocity profiles. Our power functional approximation reproduces semiquantitatively all superadiabatic forces present in the system, and, to the best of our knowledge, it constitutes the first theoretical approach that describes lane formation from first principles.

Although we have focused our study on lane formation, our approach to analyze the force contributions and our power functional approximation are general and applicable to other out-of-equilibrium situations in multicomponent colloidal systems. In particular, our work paves the way to study superadiabatic forces inside channels where ratchet-like wall shapes can determine the direction of the motion within the lanes⁵¹, and the description using power functional theory of the complete dynamics of sedimentation of binary mixtures from an initial uniform state to the formation of the sedimentation-diffusion-equilibrium stacking sequence^52,53. The power functional approach presented here can also help to understand the differences and similarities between lane formation in oppositely driven mixtures and in one-component sheared systems^54,55. There, the formation of lanes in steady state is induced by a shear field and likely sustained by a superadiabatic force. It would be interesting to relate the laning effect to the physics of active Brownian particles. The power functional approach is well equipped for systematic investigation of the underlying coupled nature of the nonequilibrium many-body systems, as the different types of driving (whether active^44,56,57,58 or external as in the present case) enter explicitly; the superadiabatic free power functional, however, is a universal object, independent of the specific form of driving.

Methods

Adiabatic construction

The internal force field is a functional of the density and the velocity profiles and can be split into adiabatic and superadiabatic contributions. Mathematically, the adiabatic force is the term that depends only on the density profile. Physically, the adiabatic force represents the internal force of an equilibrium system that shares the same density profile(s) as the actual out-of-equilibrium system. In ref. ⁴¹ Fortini et al. presented the first simulation method for sampling the superadiabatic forces in a monocomponent BD nonequilibrium system. The method relies on numerically finding the adiabatic external potential that generates the desired density profile in equilibrium. A related and improved iterative (custom flow) method to numerically find the adiabatic potential is presented in ref. ⁴⁰. We follow here another approach to find the adiabatic potential, based on classical density functional theory (DFT)⁴⁷. In DFT the adiabatic force is a functional of the one-body density profiles, and given by the functional derivative

$${{\bf{f}}}_{{\rm{ad}}}^{\alpha }=-\nabla \frac{\delta {F}_{{\rm{exc}}}[\{{\rho }_{\alpha ^{\prime} }\}]}{\delta {\rho }_{\alpha }},$$

(60)

where ${F}_{{\rm{exc}}}[\{{\rho }_{\alpha ^{\prime} }\}]$ is the excess (over ideal) Helmholtz free energy functional.

The adiabatic system is obtained by finding the species-dependent adiabatic external potentials ${V}_{{\rm{ad}}}^{\alpha }(x,t)$, which are constructed such that the equilibrium one-body densities for both species in the adiabatic (equilibrium) system equal those in the nonequilibrium system. For the laned system, the system is homogeneous in the flow $\hat{{\bf{y}}}$-direction and therefore the adiabatic potentials can only depend on the x—coordinate. In the adiabatic (equilibrium) system, the force balance reads

$$-{k}_{{\rm{B}}}T\nabla \mathrm{ln}{\rho }_{\alpha }+{{\bf{f}}}_{{\rm{ad}}}^{\alpha }(x)-\nabla {V}_{{\rm{ad}}}^{\alpha }(x)=0.$$

(61)

In order to sample the adiabatic reference system of the laned state, we use Monte Carlo simulations with the same values of ${{\mathcal{N}}}_{\alpha }$, L, and k_BT as the nonequilibrium BD system. We initialize the adiabatic system with a configuration taken from the BD simulations and then equilibrate the system. The maximum displacement each particle is allowed to move in one Monte Carlo step is set such that 25–30% of the moves are accepted. We then sample data for  ~10¹¹ Monte Carlo sweeps (MCS). A MCS is defined as an attempt to individually move all particles in the system.

In refs. ^40,41 iterative methods for finding the adiabatic potential that enters in Eq. (61) are presented. An initial guess for the potential is iteratively adjusted until the density profile equals the corresponding counterpart in the nonequilibrium system. Here, we have to match two density profiles, one for each species. As the modulation of the total density is small, a calculation of the adiabatic potentials using a simple density functional based on the equation of state of scaled-particle theory (SPT) for hard-disks yields excellent results. Within SPT the compressibility factor of an homogeneous system of hard disks is Z = (1 − η)^−2 59, with η = ρ_bπσ²∕4 the total packing fraction of the system. Here, ρ_b is the bulk number density. An expression for the excess free energy per particle ψ_exc follows from Z via integration⁶⁰:

$${\psi }_{{\rm{exc}}}({\rho }_{b})=\int _{0}^{\eta }{\mathrm{d}}\eta ^{\prime} \frac{Z-1}{\eta ^{\prime} }=-\mathrm{ln}(1-\eta )+\frac{\eta }{1-\eta }.$$

(62)

Using the above expression we construct an approximation for the excess free energy functional:

$$\beta {F}_{{\rm{exc}}}[\rho ]=\int {\mathrm{d}}x\rho (x){\psi }_{{\rm{exc}}}(\bar{\rho }),$$

(63)

in which the excess free energy is evaluated in a weighted density given by

$$\bar{\rho }(x)=\frac{1}{\sigma } \int _{x-\sigma /2}^{x+\sigma /2}{\mathrm{d}}x^{\prime} \rho (x^{\prime} ).$$

(64)

Since the mixture is symmetric only the total density ρ = ρ₁ + ρ₂ enters the expression for the excess free energy. The resulting adiabatic force f_ad via functional differentiation is species-independent and directed along the gradient direction

$${{\bf{f}}}_{{\rm{ad}}}(x)={{\bf{f}}}_{{\rm{ad}}}^{\alpha }(x)=-\nabla \frac{\delta {F}_{{\rm{exc}}}[\rho ]}{\delta \rho (x)}.$$

(65)

Equations (61) and (65) yield an expression for the species-dependent adiabatic potential

$$\begin{array}{ll}\beta {V}_{{\rm{ad}}}^{\alpha }(x)\,=\,-{\mathrm{ln}}\left({\rho }_{\alpha }{\Lambda }_{\alpha }\right)-{\psi }_{{\rm{exc}}}(\bar{\rho })\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-\frac{1}{\sigma }\int dx^{\prime} \rho (x^{\prime} )\frac{\partial {\psi }_{{\rm{exc}}}(\bar{\rho })}{\partial \bar{\rho }},\end{array}$$

(66)

where Λ_α is the (irrelevant) thermal wavelength of species α and the spatial argument of $\bar{\rho }$ and η was omitted for clarity.

Figure 7a displays results for the adiabatic potentials. As expected, the potential for each species is high in the region of the respective minority lane and low in the majority lane, with a total difference between maximum and minimum of approximately 1.3ϵ (equal to 2.6k_BT). The external forces created by the two adiabatic potentials segregate the two species in the equilibrium system. In Fig. 7b, we present the total density sampled in the adiabatic system and that in the nonequilibrium system. Both systems have almost the same density profiles, which shows the quality of the adiabatic construction. The species-dependent density profiles in the adiabatic system are also in very good agreement with those in out-of-equilibrium (not shown).

Fig. 7: Adiabatic potentials.

a Adiabatic external potentials ${V}_{{\rm{ad}}}^{\alpha }$ in the adiabatic (equilibrium) system for species 1 (orange solid) and species 2 (blue dashed) as a function of x. b Total density profile ρ sampled via Brownian Dynamics in nonequilibrium steady state (black dotted) and via Monte Carlo in the corresponding adiabatic equilibrium system (violet solid).