# Competition between lanes and transient jammed clusters in driven binary
mixtures

Honghao Yu  Yusuf Hamied Department of Chemistry, University of Cambridge,
Lensfield Road, Cambridge CB2 1EW, United Kingdom  Robert L. Jack  Yusuf
Hamied Department of Chemistry, University of Cambridge, Lensfield Road,
Cambridge CB2 1EW, United Kingdom  Department of Applied Mathematics and
Theoretical Physics, University of Cambridge, Wilberforce Road, Cambridge CB3
0WA, United Kingdom

###### Abstract

We consider mixtures of oppositely driven particles, showing that their non-
equilibrium steady states form lanes parallel to the drive, which coexist with
transient jammed clusters where particles are temporarily immobilised. We
analyse the interplay between these two types of non-equilibrium pattern
formation, including their implications for macroscopic demixing perpendicular
to the drive. Finite-size scaling analysis indicates that there is no critical
driving force associated with demixing, which appears as a crossover in finite
systems. We attribute this effect to the disruption of long-ranged order by
the transient jammed clusters.

##  I Introduction

Non-equilibrium systems can form complex patterns in their steady states,
including flocking of birds [1, 2, 3], active nematics [4, 5], and motility-
induced phase separation [6, 7, 8]. These patterns are sustained by the
continuous injection of energy into the system, which is dissipated as heat.
The formation and evolution of large coherent structures can sometimes be
captured by deterministic equations [9, 10, 11], but some chaotic fluctuating
patterns necessitate a stochastic description, using methods of non-
equilibrium statistical physics. These systems challenge existing theories:
How should emergent patterns be measured and classified? Can one identify
distinct dynamical phases, analogous to equilibrium systems? Can the patterns
be controlled by external perturbations?

A common example of a fluctuating pattern is the emergence of _lanes_ , when
two species of particles are driven in opposite directions: particles of the
same species tend to follow each other. In physics, this occurs for colloidal
systems in electric fields [12, 13, 14], and in plasmas [15, 16, 17, 18], and
ionic liquids [19]. Similar structures are also familiar from pedestrian flow
[20, 21, 22, 23, 24], and ant foragers [25], and oscillatory driving can also
induce other patterns [26, 27, 28].

Laning behaviour can be captured in computer simulations of model systems
where particles interact by repulsive short-ranged potentials and move by
Brownian dynamics. Experiments on laning involve more complex interactions,
such as electrostatic or hydrodynamic interactions in colloids [13, 14, 27],
or non-trivial decision-making by pedestrians [20, 21, 22, 23]. Nevertheless,
simple models can capture generic laning behaviour, indicating that this
robust collective phenomenon [12, 29, 30, 31, 32, 33] is amenable to the
methods of statistical physics.

A generic mechanism for laning is enhanced lateral diffusion (ELD) [34, 35,
36, 32, 37, 38, 39, 40]: collisions between oppositely driven particles
promote rapid diffusive motion perpendicular to the field. This favours
accumulation in lanes, where such collisions are avoided. Despite considerable
work on these non-equilibrium steady states, open questions remain about their
generic properties, such as the possible existence of long-ranged order, and
de-mixing of the system into domains, perpendicular to the drive [30, 31, 32,
41, 39].

This paper presents computer simulations of a model system that exhibits
laning. The main phenomenology is outlined in Fig. 1, with details given
below. Red and blue particles are driven in opposite directions by a force
that is parameterised by the Peclet number (Pe). Increasing Pe from zero, one
observes increasingly inhomogeneous steady states (labelled II-IV), which may
be contrasted with the homogeneous equilibrium mixture at $\mathrm{Pe}=0$
(state I). We find three main results: (i) Using extensive simulations for
different system sizes and driving strength, we characterise for the interplay
between the laning state III and the demixed state IV. (ii) While state IV
occurs in finite systems, we show that sufficiently large systems are always
mixed: hence demixing is not a phase transition. (iii) We characterise the
large density fluctuations in states II and III in terms of transient jammed
clusters (TJCs), which tend to disrupt the lanes, suppressing long-ranged
order. Previous work [30] analysed states II and III by focussing on long-
ranged order parallel to the drive; here we analyse phase separation
perpendicular to the drive, showing that state IV is suppressed in large
systems. Hence, we interpret laning behaviour in these systems as a
combination of ELD (favouring lane formation), and transient clustering
(disrupting laning).

Figure 1: (a) Example of a transient jammed cluster. (b) Steady state
configurations of a system with $L_{\parallel}=150d_{0}$, varying
$\mathrm{Pe}$, as shown. Supplemental Material movies 1-3 show trajectories
from steady states II-IV [42].

##  II Model

To capture the generic phenomenology of laning – as it appears in a variety of
practical contexts [12, 29, 30, 31, 32, 33] – we adopt a simplified model that
is sufficient to capture the behavior of interest. Similar to previous works
[12, 31, 30, 32, 23, 33, 41, 24], the model captures the fact that two species
of particle are driven in opposite directions and that they feel strong
repulsive interactions that prevent them from overlapping. Their motion is
biased by the driving field but also includes a noise term, which can describe
Brownian motion in colloidal systems or other random influences from the
environment. Specific laning systems might include other interactions, such as
electrostatic forces in colloids or more complex decision rules for
pedestrians, but these are not generically required for laning to occur, so we
neglect them, consistent with [12, 31, 30, 32, 23, 33, 41, 24].

To this end, we consider a binary mixture of $N$ particles, with $N/2$ in each
species, occupying a two dimensional box of size $L_{\parallel}\times
L_{\perp}$, with periodic boundary conditions. The species of particle $i$ is
labelled as $\sigma_{i}=+1$ (blue) or $\sigma_{i}=-1$ (red). Particles are
driven by a force of strength $E$ along the $x$-axis and interact with each
other via WCA potentials [43],

$V(r)=4\epsilon[(d_{0}/r)^{12}-(d_{0}/r)^{6}+(1/4)]\Theta(2^{1/6}d_{0}-r),$
(1)

where $d_{0}$ is the diameter, $\epsilon$ is the repulsion strength, and
$\Theta$ is the Heaviside function. The position $\bm{r}_{i}=(x_{i},y_{i})$ of
particle $i$ follows Langevin dynamics, implemented in LAMMPS [44],

$m\ddot{\bm{r}}_{i}=-\gamma\dot{\bm{r}}_{i}-\nabla_{i}U+\sigma_{i}E\bm{\hat{x}}+\sqrt{2\gamma
k_{B}T}\,\bm{\eta}_{i},$ (2)

where $\gamma$ is a friction constant, $\bm{\hat{x}}$ is a unit vector in the
$x$-direction, $T$ is the temperature of the heat bath, $k_{B}$ is Boltzmann’s
constant, $U$ is the interaction energy, and $\bm{\eta}_{i}$ is a Gaussian
white noise with mean zero and
$\langle{\eta}_{i,\alpha}(t){\eta}_{j,\beta}(t^{\prime})\rangle=\delta_{ij}\delta_{\alpha\beta}\delta(t-t^{\prime})$,
where $\alpha,\beta$ indicate Cartesian components. The Brownian time is
$\tau_{\rm B}=d_{0}^{2}\gamma/(k_{B}T)$.

We non-dimensionalise the system using base units $d_{0},m,\epsilon$, see
Appendix A for full details. The dimensionless control parameters of the model
are the Peclet number $\text{Pe}=Ed_{0}/(k_{B}T)$, the area fraction
$\phi=N\pi d_{0}^{2}/(4L_{\parallel}L_{\perp})$, the reduced friction
$\tilde{\gamma}=\gamma d_{0}/\sqrt{m\epsilon}$, the reduced temperature
$\tilde{T}=k_{\rm B}T/\epsilon$ and the aspect ratio of the simulation box
$S=L_{\parallel}/L_{\perp}$. We take $\tilde{\gamma}=1000$ to mimic a high-
friction colloidal environment, and $\tilde{T}=1$ (results depend weakly on
this parameter). Since particles explore space much more quickly in the
$x$-direction, we take $S=5$, except where explicitly stated otherwise. We
focus on $\phi=0.35$ which is representative of a laning regime
$0.1\lesssim\phi\lesssim 0.4$, some results for other densities are presented
in the Appendices.

Figure 2: (a) $\psi_{\rm dem}$ against Pe at different system sizes. Dotted
and solid lines are obtained from disordered and demixed initial conditions,
respectively. (b) $h_{\sigma,\sigma}(\bm{r})$ plotted along the $x$-direction.
The solid line indicates $h_{\sigma,\sigma}(x,0)\sim x^{-3/2}$. (c) Dependence
of $h_{\rho,\rho}(\bm{r})$ on $\bm{r}=(x,y)$. In (b), (c),
$L_{\parallel}=150d_{0}$.

##  III Results

###  III.1 Steady states: Demixing is a smooth crossover

The steady state behaviour of this model is illustrated in Fig. 1(b), for a
range of Pe. We verified that all systems have converged to their steady
states by starting independent simulations from both homogeneous and fully de-
mixed initial conditions, which lead to the same results. This requires long
simulations, up to $4000\tau_{B}$ (see Appendix A for further details of
initialisation and equilibration).

To characterise the fully demixed state, we interpret $\sigma_{i}$ as a
dimensionless charge and evaluate the Fourier transform of the charge density
at $\bm{k}^{*}=(0,2\pi/L_{\perp})$ [45, 46, 39]:

$\psi_{\rm
dem}=\langle|\Psi({\bm{k}}^{*})|\rangle=\Big{\langle}\Big{|}\frac{1}{L_{\perp}L_{\parallel}}\sum_{i=1}^{N}\sigma_{i}e^{-i\bm{k}^{*}\cdot\bm{r}_{i}}\Big{|}\Big{\rangle},$
(3)

Here and throughout, angle brackets indicate averages in the steady state of
the dynamics.

We use $\psi_{\rm dem}$ as an order parameter for demixing, in a finite-size
scaling analysis, by increasing $N,L_{\parallel},L_{\perp}$, at fixed
$\mathrm{Pe},\phi$, and aspect ratio $S$. A demixing phase transition would
manifest as a critical Peclet number $\mathrm{Pe}^{*}$: For
$\mathrm{Pe}<\mathrm{Pe}^{*}$ the state would be homogeneous and
$\Psi({\bm{k}}^{*})$ would be small in modulus (of order $N^{-1/2}$ as
$N\to\infty$). For $\mathrm{Pe}>\mathrm{Pe}^{*}$ the state would be demixed:
the complex number $\Psi({\bm{k}}^{*})$ has a modulus of order $1$; its phase
is random and indicates the relative positions of the two domains. Hence,
demixing involves spontaneous symmetry breaking, via the phase of $\Psi$.

Figure 3: (a) $S_{\sigma}(\bm{k})$ at various $\rm{Pe}$ with
$L_{\parallel}=300$ and $\phi=0.35$, corresponding to $N=8100$. The red arrows
indicate the locations of the highest peaks. (b) $S_{\sigma}(0,k_{y})$ at
different $\rm{Pe}$ with $L_{\parallel}=300$. The peaks correspond to the
$k^{*}_{y}$ in (c). (c) The characteristic wavevector $k^{*}_{y}$ extracted
from $S_{\sigma}(\bm{k})$ against different $\rm{Pe}$. (d) The width of the
lobes extracted from $S_{\sigma}(\bm{k})$ against different $\rm{Pe}$.

The finite-size scaling results in Fig. 2(a) are _not_ consistent with any
such transition. For each system size, there is a crossover from homogeneous
to de-mixed states on increasing $\mathrm{Pe}$. However, increasing $N$ shows
that ever larger values of $\mathrm{Pe}$ are required to see demixing: there
is no sharp threshold $\mathrm{Pe}^{*}$ at which long-ranged order appears.
Figure 2(a) shows that these results are independent of whether the system is
initialised in a homogeneous (mixed) or fully demixed state.

Figure 2(a) suggests a correlation length which grows with $\mathrm{Pe}$ but
never diverges. We identify it with the width of the lanes (perpendicular to
the drive), which is $L_{\perp}/2$ in state IV, but much smaller in state III,
these correlations are discussed further below. A similar scenario was
proposed in [30] for correlations parallel to the drive; our focus here on
perpendicular correlations is essential for detecting the fully de-mixed state
IV, which resembles macroscopic phase separation at equilibrium. Comparing
with that case, our observation of demixing in large finite systems without
any phase transition is a distinctively non-equilibrium effect.

###  III.2 Spatial correlations

####  III.2.1 Real-space analysis

States II and III of Fig. 1 exhibit laning, but the particle density also
develops inhomogeneities, where right- and left-moving particles tend to block
each other. To investigate these correlations, define density-density and
“charge-charge” correlation functions:

$\displaystyle h_{\rho,\rho}(\bm{r})$
$\displaystyle=\langle\rho(\bm{{r}})\rho(\bm{{0}})\rangle/\overline{\rho}^{2}-1$
(4) $\displaystyle h_{\sigma,\sigma}(\bm{{r}})$
$\displaystyle=\langle\sigma(\bm{{r}})\sigma(\bm{{0}})\rangle/\overline{\rho}^{2}$

where $\rho(\bm{r})=\sum_{i=1}^{N}\delta(\bm{r}_{i}-\bm{r})$ is the empirical
particle density,
$\sigma(\bm{r})=\sum_{i=1}^{N}\sigma_{i}\delta(\bm{r}_{i}-\bm{r})$ is the
“charge” density, and $\bar{\rho}=N/(L_{\parallel}L_{\perp})$ is the mean
density.

Figure 4: $S_{\rho}(\bm{k})$ at various $\rm{Pe}$ with $L_{\parallel}=300$ and
$\phi=0.35$, corresponding to $N=8100$. The features near the origin
correspond to the emergence of TJCs.

The behaviour of these correlation functions on hydrodynamic length scales has
been analysed [47, 33, 41, 48] within a mean-field approximation, which is
valid when density fluctuations are not too large, for example at high density
with weakly interacting particles [33]. For the system considered here, the
particles have strongly repulsive cores and Fig. 1 shows that density
fluctuations are large. Hence mean-field theory is not expected to be
quantitative: it does not predict the demixing effect seen in state IV. Still,
the universal hydrodynamic correlations of Refs. [33] should emerge on large
scales, within homogeneous states. Figure 2(b) shows the charge
autocorrelation function measured parallel to the drive. For $\mathrm{Pe}=20$,
these results are consistent with Refs. [33], which predicted a decay as
$|x|^{-3/2}$. For larger Pe, the correlations grow smoothly with Pe, as the
lanes become more pronounced.

####  III.2.2 Correlations in Fourier space

It is useful to display these correlations in Fourier space. For the charge
density $\sigma$, the structure factor associated with two-point correlations
is

$S_{\sigma}(\bm{k})=\frac{1}{N}\Big{\langle}\Big{|}\sum^{N}_{j=1}\sigma_{j}\exp(-i\bm{k}\cdot\bm{r}_{j})\Big{|}^{2}\Big{\rangle},$
(5)

where $\bm{k}=(k_{\parallel},k_{\perp})$ is the wavevector.

Results are shown in Fig. 3. Figure 3(a) shows that these functions have a
characteristic shape, with strongly anisotropic correlations that were already
anticipated in Chapter 3.1 of Refs. [49], see also Sec. III.2.3 below. (The
system is not phase separated for these parameters and system sizes.)
$S_{\sigma}$ has two lobes that lie along the $k_{y}$-axis, which get stronger
as Pe increases. To characterise these lobes, we consider
$S_{\sigma}(0,k_{y})$ [Fig. 3(b)] and we denote its maxima by $\pm k_{y}^{*}$.

We identify $\pi/k_{y}^{*}$ as one half of the wavelength of the modulations
in $\sigma(\bm{r})$, which is the typical lane width. We also measure the
width of the lobes in $S_{\sigma}$, by plotting $S_{\sigma}(k_{x},k_{y}^{*})$
as a function of $k_{x}$ and computing its full width at half maximum. This is
denoted by $k^{*}_{\rm width}$ and represents the range of correlations _along
the lanes_.

Figure 3(c) shows how $k_{y}^{*}$ depends on Pe and system size (always at
fixed density). Similarly Fig. 3(d) shows $k^{*}_{\rm width}$. Several points
are notable: For small-to-moderate Pe, the length scale $\pi/k_{y}^{*}$ is a
few particle diameters, independent of system size. For the largest Pe and the
smaller system sizes, one sees a change in behaviour of $k_{y}^{*}$, which
corresponds to the onset of macroscopic phase separation. At the smallest Pe,
the correlations are weak but appear rather long-ranged, we return to this
point in Sec. III.2.3, below. We also find that $k_{\rm width}^{*}$ is much
smaller than $k_{y}^{*}$, consistent with the fact that correlations along the
lane have a range that is much longer than the lane width. For large Pe, the
width starts to depend on the system size, which naturally limits the range of
correlations (clearly, $k_{\rm width}^{*}\geq\pi/L_{\parallel}$).

Figure 5: (a) $S_{\sigma}(\bm{k})$ from both mean-field theory and simulation.
(b) $S_{\rho}(\bm{k})$ from both mean-field theory and simulation. Both
$S_{\sigma}(\bm{k})$ and $S_{\rho}(\bm{k})$ are shown with $\rm{Pe}=10$ and
$\phi=0.35$. For simulations, $L_{\parallel}=300$, corresponding to $N=8100$.

To complement the analysis of charge correlations, we show the structure
factor for the density field:

$S_{\rho}(\bm{k})=\frac{1}{N}\Big{\langle}\Big{|}\sum^{N}_{j=1}\exp(-i\bm{k}\cdot\bm{r}_{j})\Big{|}^{2}\Big{\rangle}\;,$
(6)

which is proportional to the Fourier transform of $h_{\rho,\rho}$, see Fig.
2(c). This function is shown in Fig. 4. On increasing Pe, a significant peak
appears at small $k$, signalling the formation of TJCs and associated large
density fluctuations. (The compressibility $\tilde{\chi}$ can be expressed in
terms of an orientational integral over $S_{\rho}$.) However, this peak also
has a significant orientation dependence, which is discussed in the next
subsection by comparing with results of mean-field theory.

####  III.2.3 Comparison with mean-field theory

It is instructive to compare these results in Fourier space with those of
mean-field theory [33]. Based on Eq. (9) from Refs. [33], one can obtain the
structure factors $S_{\sigma}(\bm{k})$ and $S_{\rho}(\bm{k})$:

$\displaystyle S_{\sigma}(\bm{k})$
$\displaystyle\sim\frac{k^{4}+A_{\sigma}k^{2}_{\parallel}\rm{Pe}^{2}}{k^{4}+B_{\sigma}k^{2}_{\parallel}\rm{Pe}^{2}},$
(7) $\displaystyle S_{\rho}(\bm{k})$
$\displaystyle\sim\frac{A_{\rho}k^{2}_{\parallel}\rm{Pe}^{2}}{k^{4}+B_{\rho}k^{2}_{\parallel}\rm{Pe}^{2}},$

where $A_{\sigma},B_{\sigma},A_{\rho},B_{\rho}$ are positive constants that
depend on the system of interest via a single parameter $\tilde{v}_{0}$, which
is the Fourier transform of the interaction potential evaluated at wavevector
$\bm{k}=0$. For the purposes of this paper, it only matters that these are
positive constants, of order unity.

Figure. 5 shows representative results from simulations above, compared with
results of the theory, for the representative parameter value
$\tilde{v}_{0}=1$, which leads to $A_{\sigma},B_{\sigma},A_{\rho},B_{\rho}$ of
order unity. The theory captures the anisotropic correlations at small $k$.
However it does not capture the peaks that appear in the structure factor, for
example it does not make predictions for $k_{y}^{*}$, nor does it capture the
peak in $S_{\rho}$ (which we attributed previously to TJCs). An interesting
prediction of Eqs. (7) is that small Pe leads to weak correlations with large
correlation lengths that scale as $\xi\sim 1/({\rm Pe}\sqrt{C})$ where $C$ is
one of $A_{\sigma},A_{\rho},B_{\sigma},B_{\rho}$. It would be interesting to
investigate these length scales further, but the differences between theory
and simulation are strong enough that we do not attempt this here. (It is not
clear to us if this theoretical prediction is linked to the increase of
$k_{y}^{*}$ with Pe, observed for small Pe in Fig. 4.)

Figure 6:  (a) $\tilde{\chi}-\tilde{\chi}^{\text{eq}}$ against Pe, for
different system sizes, with cut-off radius $R=5d_{0}$. The mean-field theory
predicts that this quantity vanishes for $R\to\infty$. (b) Probability
distribution for the number of particles in a randomly-chosen circular probe
region of radius $r_{0}=1.75d_{0}$, for various Pe. (c) Probability
distributions for the $x$-displacement $\Delta x$ of each species of particles
within time $\tau=0.25\tau_{B}$, at $\rm{Pe}=50$. The dashed line indicates
$P(\Delta x)\propto e^{-(\Delta x_{\text{free}}-\Delta x)/\ell}$. For (b),
(c), $L_{\parallel}=300d_{0}$.  Figure 7: (a) The characteristic trap time
$\tau_{\rm TJC}$ in units of $\tau_{B}$ at different $\rm Pe$ number when TJCs
are present with $L_{\parallel}=300d_{0}$. (b) The characteristic trap time
$\tau_{\rm TJC}$ against measurement window $\tau$ of system as in (a). The
characteristic trap time $\tau_{\rm TJC}$ converges with sufficiently large
measurement window $\tau$.

###  III.3 Transient jammed clusters

####  III.3.1 Characterisation of TJCs

States II and III in Fig. 1 show strikingly large density fluctuations
involving clusters of particles with opposite charge, which block each others’
motion, leading to a transient jamming phenomenon. To analyse this, Fig. 2(c)
shows density correlations for $\mathrm{Pe}=(0,50,100)$, corresponding to
states I, III, and IV in Fig. 1. For the equilibrium fluid state
($\mathrm{Pe}=0$), the correlations are isotropic and short-ranged. For
$\mathrm{Pe}=50$, the correlations are stronger (as indicated by the darker
shading) and they are also longer ranged – this is a signature of the TJCs.
For $\mathrm{Pe}=100$, the correlations resemble those for $\mathrm{Pe}=0$.
This re-entrant behaviour occurs because the demixed domains are similar to
equilibrium systems that are being advected at a constant velocity. This
advection also affects the usual nearest-neighbour correlations in
$h_{\rho,\rho}$, leading to to bands in $h_{\rho,\rho}$, extended along the
driving direction.

To further quantify the TJCs, we integrate the correlation function over space
to obtain an analog of the compressibility of an equilibrium system [50]:

$\tilde{\chi}=1+\overline{\rho}\int_{|\bm{r}|<R}h_{\rho,\rho}(\bm{r})d\bm{r}\;.$
(8)

The large-distance cutoff $R$ reduces statistical uncertainty, Appendix C
shows that these results depend weakly on $R$.

Figure 6(a) shows $\tilde{\chi}$, measured relative to its equilibrium value.
It grows rapidly as Pe increases from zero, indicating the formation of TJCs.
There is a plateau at intermediate Pe, showing that these clusters persist
across broad range of parameters. At large Pe, $\tilde{\chi}$ decreases, due
to the re-entrant demixing, recall Fig. 2(c). The mean-field theory of Refs.
[33] does not capture these TJCs, because of its assumption of small density
fluctuations. For example, the correlations predicted by that theory lead to
$\tilde{\chi}\to 0$ as $R\to\infty$ (see Appendix C): this prediction is not
consistent with the data presented here.

As an alternative characterisation of TJCs, Fig. 6(b) shows the probability
distribution of the number of particles $N_{\rm in}$ in a randomly chosen
circular probe region [51, 52, 53]. The variance of this distribution
increases dramatically with $\mathrm{Pe}$, with tails at large $N_{\rm in}$
reflecting the existence of TJCs. These clusters also have a dynamical
signature: Fig. 6(c) shows the distribution of particle displacements parallel
to the driving field. Each particle displacement is measured as

$\Delta x_{i}(\tau)=x_{i}(t_{0}+\tau)-x_{i}(t_{0}),$ (9)

where $t_{0}$ is an arbitrary initial time. Then the results are averaged
separately over the two species. The distributions are skewed, with
significant weight near zero displacement, reflecting particles which remain
almost stationary in TJCs. Similar slowing down has been observed in $3$D [54,
55, 56].

Figure 8: (a)-(c) Time evolution of the charge density profile, starting from
a demixed state, for $L_{\parallel}=100d_{0}$ and $S=1.25$. Symbols are
numerical results and dotted lines are fits to the diffusion equation, see
Appendix D for details of the fitting procedure. (d) Time series, highlighting
a TJC that nucleates new lanes. We show a subsystem from large simulation with
$L_{\parallel}=300d_{0}$ and ${\rm Pe}=78$.

####  III.3.2 Residence time for particles in TJCs

The characteristic trapping time of particles in TJCs can be estimated from
the distribution of longitudinal displacement in Fig. 6(c), using the
following simple model. The measurement window for the displacement is $\tau$
and we consider the species with positive charge, so a free particle (outside
of the TJC) travels an average distance approximately $\Delta
x_{\text{free}}=\tau E/\gamma$ within this time. (The average velocity of a
free particle is $v_{\text{free}}=E/\gamma$.)

From Fig. 6(c), the distribution of displacements $P(\Delta x)$ has an
exponential tail for small $\Delta x$, we write

$P(\Delta x)\approx Ae^{-(\Delta x_{\text{free}}-\Delta x)/\ell}$ (10)

where $A$ is a constant, $\ell$ is a characteristic length scale, which is an
estimate of the length by which a typical particle in a TJC is “held back” by
its temporary immobilisation. (The value of $A$ presumably depends on how many
particles participate in TJCs.) As a simple method to obtain a characteristic
trapping time we divide $\ell$ by the characteristic velocity to obtain
$\tau_{\rm TJC}=\ell/v_{\text{free}}$.

To better interpret this time scale, suppose that particles in TJCs stay there
for an exponentially distributed time $t_{\rm TJC}$ with mean $\tau_{\rm
TJC}$, during which time they do not move at all parallel to the field. Their
resulting displacement is then $\Delta x_{\rm free}-v_{\text{free}}t_{\rm
TJC}$. Assuming also that $P(t_{\rm TJC})\propto{\rm e}^{-t_{\rm
TJC}/\tau_{\rm TJC}}$, the tail of the displacement distribution is readily
seen to be Eq. (10) with $\ell=\tau_{\rm TJC}v_{\text{free}}=\tau_{\rm
TJC}d_{0}{\rm Pe}/\tau_{\rm B}$.

However, this simple analysis neglects effects of the measurement time window
$\tau$. In particular, the distribution of $\Delta x$ is not sensitive to
trapping times longer than $\tau$. As a result, the above argument tends to
underestimate the trapping time. Figure 7(a) shows the resulting $\tau_{\rm
TJC}$ in the regime $\rm Pe\in[30,70]$. These results indicate that large
measurement times $\tau$ can be used to obtain robust estimates for the
trapping time. Specifically, taking $\text{Pe}=50$ and $\tau=\tau_{B}$, we
obtained $\ell=5.4d_{0}$ from the exponential fit and therefore have
$\tau_{\rm TJC}=0.11\tau_{B}$. While this is smaller than the Brownian time,
it is significantly larger than the the time for a particle in a lane to move
its own diameter, which is $\tau_{B}/{\rm Pe}\approx 0.020\tau_{B}$. However,
smaller $\tau$ results in underestimates of the true $\tau_{\rm TJC}$, as
expected, see Fig. 7(b). Nevertheless, the resulting estimates of $\tau_{\rm
TJC}$ are all of a similar order of magnitude; they also depend weakly on Pe,
within the range $30<{\rm Pe}<70$ where TJCs play a significant role.

###  III.4 Dynamics of mixing – the role of TJCs

We have shown that these driven systems are always mixed in the limit of large
systems and that the mixed state supports large TJCs. We now establish a
connection between these two results.

We choose parameters such that the system’s steady state is mixed but we
initialise it in a demixed state, so that we can follow the mixing process. To
collect good statistics while maintaining a sufficiently large $L_{\perp}$, we
reduce the aspect ratio to $S=1.25$. Figure 8 shows time-dependent composition
profiles perpendicular to the drive, for various Pe. In the absence of driving
($\mathrm{Pe}=0$), mixing happens by diffusion. To illustrate this, we fit the
results to the (analytical) solution of the diffusion equation
$\partial_{t}\langle\sigma\rangle=D_{y}\partial_{y}^{2}\langle\sigma\rangle$,
using $D_{y}$ as a single fitting parameter (see Appendix D for details of the
fitting procedure). The fit is excellent. For $\mathrm{Pe}=20$, the system
still mixes readily: the diffusion equation still fits the data adequately but
the value of $D_{y}$ is larger, reflecting the ELD [32, 39].

However, for $\mathrm{Pe}=60$, the lifetime of the demixed state is longer:
the interface between the two domains is stable on short time scales, because
collisions between oppositely charged particles tend to deflect them back into
their own domains. On longer time scales, mixing does occur, but the mechanism
is not diffusive, as signalled by systematic deviations from the diffusive
theory in Fig. 8(c). In particular, the composition inside the demixed domains
relaxes more quickly than diffusive mixing would predict, indicating that a
collective mechanism is disrupting the domains.

This mixing mechanism is illustrated in Fig. 8(d) and the Supplemental
Material movies 4,5 [42]: collisions near the interface occasionally form TJCs
which disrupt the flow and mix the counterpropagating species. This initiates
a period of chaotic flow, which involves further TJCs and – eventually – the
creation of new stable lanes. This mechanism further emphasizes the rich
phenomenology of these non-equilibrium systems.

To understand the absence of demixed states in large systems, observe that
increasing $L_{\parallel}$ provides more locations where TJCs can form,
initiating the instability. This mechanism is intrinsic to the interface
between the oppositely-moving domains: Appendix D shows that increasing
$L_{\perp}$ (at fixed $L_{\parallel}$) does not change the local behaviour
near the interface. In other words, the absence of de-mixed states for large
$L_{\parallel}$ depends weakly on the aspect ratio of the system, since any
interface between counterpropagating domains is prone to fragmentation via
TJCs. Interfaces between oppositely-driven domains have been recently studied
in Refs. [57, 58], and nucleation of lanes is also reminiscent of other
complex flow phenomena [59].

Figure 9:  Time series, highlighting a TJC that causes coarsening of two
lanes. We show a subsystem from large simulation with $L_{\parallel}=150d_{0}$
and ${\rm Pe}=80$.

If the system is small enough that the steady state is fully demixed, the
mechanism of Fig. 8(d) can also operate in reverse: four lanes can merge into
two, via TJC formation at the interface, see Fig. 9. However, for large
systems, the mixing (lane creation) events predominate.

##  IV Outlook

Lane formation for oppositely driven particles relies on ELD, which tends to
cause demixing, perpendicular to the drive. We showed that this demixing is
not a phase transition – sufficiently large systems always remain mixed. They
also support large density fluctuations in the form of TJCs, which exist over
a wide range of Pe, and play an important role in mixing, by destabilising
interfaces between domains with oppositely driven particles. In this sense,
the complex patterns shown in Fig. 1 reflect a competition between ordering by
ELD and mixing via TJCs. We look forward to future work on such patterns. For
example, the possibility of system-spanning TJCs may be relevant for dynamical
arrest in small systems [60]. Similar phenomena may also be relevant for
pattern formation in other non-equilibrium systems, such as self-propelled
particles or sheared fluids [61, 62, 63, 52, 64, 65, 66, 67, 68, 69], and the
control of such patterns offers possibilities for the design of directed self-
assembly [70, 40]. At higher densities, we also anticipate an interesting
interplay between TJCs and crowding (jamming or glassy) behaviour [71, 72].

###### Acknowledgements.

We thank K. Thijssen, D. Frenkel, L. Berthier, O. Bénichou, M. Cates, and N.
Wilding for helpful discussions.

##  Appendix A Model and simulation details

We non-dimensionalise the Langevin Eq. (1) using base units
$d_{0},m,\epsilon$. The natural unit of time is
$t_{0}=d_{0}\sqrt{m/\epsilon}$. Then we define reduced quantities which are:
the particle position $\tilde{\bm{r}}_{i}=\bm{r}_{i}/d_{0}$, the Peclet number
$\text{Pe}=Ed_{0}/(k_{\rm B}T)$, the volume fraction $\phi=N\pi
d_{0}^{2}/(4L_{\parallel}L_{\perp})$, the reduced time $\tilde{t}=t/t_{0}$,
the reduced interaction energy $\tilde{U}=U/\epsilon$, the reduced friction
constant $\tilde{\gamma}=\gamma t_{0}/m$, the reduced temperature
$\tilde{T}=k_{\rm B}T/\epsilon$, and the aspect ratio of the simulation box
$L_{\parallel}/L_{\perp}$. Recall that the Brownian time is $\tau_{\rm
B}=d_{0}^{2}\gamma/(k_{B}T)=t_{0}\tilde{\gamma}/\tilde{T}$: this evaluates to
$10^{3}t_{0}$ for the parameters used in this paper.

Figure 10: (a) $\tilde{\psi}_{\rm{dem}}$ against $\rm{Pe}$ at different
packing fraction $\phi$ for $N=2026$. (b) Diagram showing the regime in which
phase separation is observed ($\tilde{\psi}_{\rm{dem}}\gtrsim 0.6$) for
$N=2026$. (c)-(e) $\tilde{\psi}_{\rm{dem}}$ against $\rm{Pe}$ at different
system sizes for $\phi=0.25,0.31,0.39$ respectively.

The non-dimensionalised Langevin equation is then

$\frac{d^{2}}{d\tilde{t}^{2}}\tilde{\bm{r}}_{i}=-\tilde{\gamma}\frac{d}{d\tilde{t}}\tilde{\bm{r}}_{i}-\tilde{\nabla}_{i}\tilde{U}+\sigma_{i}\textrm{Pe}\,\tilde{T}\hat{\bm{x}}+\sqrt{2\tilde{\gamma}\tilde{T}}\,\tilde{\bm{\eta}}_{i}(\tilde{t}),$
(11)

where $\tilde{\bm{\eta}_{i}}$ is non-dimensionalised Gaussian white noise with
mean zero and
$\langle\tilde{{\eta}}_{i,\alpha}(\tilde{t})\tilde{{\eta}}_{j,\beta}(\tilde{t}^{\prime})\rangle=\delta_{ij}\delta_{\alpha\beta}\delta(\tilde{t}-\tilde{t}^{\prime})$.

Simulations are performed in LAMMPS [44] with time step $\delta t=5\times
10^{-4}t_{0}$. For initialisation, we either use disordered or demixed initial
conditions. For the disordered case, we initialise $N/2$ red and $N/2$ blue
particles with random positions and simulate the equilibrium system ($E=0$)
for $10^{2}\tau_{B}$, after which the force $E$ is introduced. For the demixed
case, we initialise $N$ particles of a single type with random positions and
simulate the equilibrium system for $10^{2}\tau_{B}$. Then we assign types to
particles according to their $y$-coordinates: for $y_{i}\in[0,L_{\perp}/2]$ we
take $\sigma_{i}=-1$ and for $y_{i}\in[L_{\perp}/2,L_{\perp}]$ we take
$\sigma_{i}=1$. Then we turn on the force $E$. As discussed in the main text,
the steady state of the system is independent of whether disordered or demixed
initialisation is used, and all measurements are taken in the steady state.

For Figs. 1-7 and 10-12 the steady state measurements are performed after
simulation for at least $10^{3}\tau_{B}$, so that the system has relaxed to
its steady state. We confirm the steady state behavior by comparing the
simulations performed with both disordered and demixed initial conditions.
When the system phase separate, the coarsening process is extremely slow as
the interfaces are smooth and stable. When $\text{Pe}>90$, it takes
approximately $4\times 10^{3}\tau_{B}$ for the system to reach complete phase
separation, and results are collected after this time.

##  Appendix B Density dependence of $\psi_{\mathrm{dem}}$

We defined $\psi_{\rm{dem}}$ to characterise the steady state behaviour of the
system, perpendicular to the drive. In this section, we investigate the
density dependence of the steady state behaviour. We do this by fixing the
number of particles $N$ and the aspect ratio of the system $S$, while varying
the system size. We re-scale the order parameter so that results are
comparable between systems at different densities, specifically

$\tilde{\psi}_{\rm{dem}}=\frac{L_{\parallel}L_{\perp}}{N}\psi_{\rm{dem}}\;.$
(12)

To justify this choice, note that a perfectly phase separated state has
“charge” density $\sigma(\bm{r})=+\overline{\rho}$ when $y\in[0,L_{\perp}/2)$,
and $\sigma(\bm{r})=-\overline{\rho}$ when $y\in[L_{\perp}/2,L_{\perp})$, one
has

$\tilde{\psi}_{\rm{dem}}\approx\Big{\langle}\Big{|}\frac{1}{N}\int\sigma(\bm{r})e^{-i\bm{k}^{*}\cdot\bm{r}}d\bm{r}\Big{|}\Big{\rangle}=\frac{2}{\pi}\approx
0.64,$ (13)

which is independent of the total density.

Simulation results are shown in Fig. 10. Panels (a) and (b) show data for a
system of $N=2026$ particles. We identify phase-separated states as those with
$\tilde{\psi}_{\rm{dem}}>0.6$, and we determine the phase boundary (black
line) in Fig. 10(b) by identifying $\rm Pe$ closest to
$\tilde{\psi}_{\rm{dem}}=0.6$ . We emphasize that this “phase boundary” was
computed by varying the box size with $N=2026$; results will be different for
systems with different numbers of particles. One sees that phase separation
occurs across a range of concentrations $\phi$. At low concentrations
$\phi\lesssim 0.05$, particle collisions are rarer, which reduces the ELD
effect, and suppresses laning and phase separation. The tendency towards
laning is non-monotonic in $\phi$ is consistent with previous work [35].
Figures 10(c)- 10(e) show finite-size scaling results for the order parameter,
at three different densities within the laning regime. In all cases, the Pe
required for phase separation increases with $N$, indicating that sufficiently
large systems will always remain mixed, for any given Pe.

Figure 11: (a) Density fluctuation $\tilde{\chi}-\tilde{\chi}^{\text{eq}}$
against Pe at $N=2026$ with $\tilde{R}=5$ for various packing fractions. (b)
Density fluctuation $\tilde{\chi}-\tilde{\chi}^{\text{eq}}$ against Pe at
$L_{\parallel}=150d_{0}$, $L_{\parallel}=300d_{0}$ with $\tilde{R}=5,10$. The
mean-field curve is shown with $\tilde{v}_{0}=1$ (see Appendix C.3).

For larger concentrations $\phi>0.4$, we observe different qualitative
behaviour, where neither laning nor phase separation takes place. Instead the
system forms macroscopic domains of the two species whose sizes and shapes
fluctuate in time, as they flow past each other. This regime brings even
further computational challenges associated with the time required to converge
to a steady state, and the ability to control finite size effects, so we do
not address it here. (These aspects were already challenging for the results
presented here at lower density, but our extensive simulations have allowed us
to control them.)

##  Appendix C Generalised compressibility $\tilde{\chi}$

###  C.1 Dependence on density

In the main text, we showed the behavior of generalised compressibility
$\tilde{\chi}$ for $\phi=0.35$. In Fig. 11(a) we show in the laning regime the
behavior of $\tilde{\chi}$ for different packing fraction is the same. It
grows at low $\rm{Pe}$ due to the formation of TJCs and it plateaus at
intermediate $\rm{Pe}$, corresponding to the coexsitence of TJCs and lanes at
intermediate $\rm{Pe}$. At large $\rm{Pe}$, $\tilde{\chi}$ drops at large
$\rm{Pe}$, corresponding to the demixing of the particles and disappearance of
TJCs.

From Fig. 11(a), we also observe that higher packing fraction corresponds to
higher $\tilde{\chi}$ across all $\rm{Pe}$ regimes when TJCs are present. The
behavior is intuitive: at higher packing fractions particles become more
packed when moving in opposite directions and form larger TJCs.

We also note the drop in $\tilde{\chi}$ at high $\rm{Pe}$ corresponding to
demixing is density dependent: higher packing fraction requires higher
$\rm{Pe}$ to demix. The behavior is consistent with the density dependence of
the order parameter as shown in Fig. 10(a).

###  C.2 Dependence on cutoff $R$

In the main text, we defined $\tilde{\chi}$ in terms of the density
fluctuations of the system. In equilibrium systems, this $\tilde{\chi}$ is
related to the compressibility (defined in terms of dependence of the volume
on the applied pressure). In these non-equilibrium systems there is no
relation between $\tilde{\chi}$ and responses to applied pressure, so we refer
to it as a generalised compressibility.

To parameterise the cutoff in Eq. (8), we non-dimensionalise as
$\tilde{R}=R/d_{0}$. For the main text, we take $\tilde{R}=5$ to reduce
numerical uncertainties. In Fig. 11(b), we show that the behavior of
$\tilde{\chi}$ depends weakly on this cutoff. There is a small increase on
increasing $\tilde{R}$, from the tails of $h_{\rho,\rho}$ but the robust
qualitative feature is that $\tilde{\chi}$ increases strongly with $\rm Pe$,
saturates at a plateau, and then decreases.

###  C.3 Mean-field analysis

We describe here the mean-field prediction for density-density correlations,
see also Fig. 6(a) of the main text. From Ref. [33], one obtains the
hydrodynamic (or singular) part of the density-density correlation function
which we denote as $h^{s}_{\rho,\rho}(x,y)$. This hydrodynamic correlation is
valid on large length scales but it misses effects of particle packing, for
example one has $h^{s}=0$ for the equilibrium case $\mathrm{Pe}=0$ (although
of course $h\neq 0$ in this case). Hence we subtract the equilibrium part of
$\tilde{\chi}$ when comparing the numerical results with the theory.

The singular part of the correlation function in two spatial dimensions can be
obtained from Ref. [33] as

$h^{s}_{\rho,\rho}(x,y)=\mathrm{Pe}^{1/2}\frac{H_{\rho,\rho}}{|x|^{3/2}}\left(\frac{\text{Pe}\,y^{2}}{{\cal
D}|x|}-1\right)\exp\left(-\frac{\mathrm{Pe}\,y^{2}}{{2\cal D}|x|}\right),$
(14)

where $H_{\rho,\rho}$ and ${\cal D}$ are constants that are fully determined
by the parameter $\tilde{v}_{0}$ that was discussed in Sec. III.2.3. We always
have ${\cal D}>0$; for $\tilde{v}_{0}\lesssim 1.3$ then $H_{\rho,\rho}>0$
while for larger $\tilde{v}_{0}$ we have $H_{\rho,\rho}<0$.

We integrate the singular part of the density correlation function to get the
singular part of the generalised compressibility,
$\tilde{\chi}^{s}_{\text{mf}}=\overline{\rho}\int_{|{\bm{r}}|<R}h^{s}_{\rho,\rho}({\bm{r}})d{\bm{r}}$.
We obtain

$\tilde{\chi}^{s}_{\text{mf}}=\overline{\rho}H_{\rho,\rho}\sqrt{\text{Pe}}\int^{R}_{0}\int^{\sqrt{R^{2}-x^{2}}}_{0}\bigg{[}\frac{1}{x^{3/2}}\left(\frac{\text{Pe}\,y^{2}}{2{\cal
D}x}-1\right)\\\ \times\exp\left({-\frac{\text{Pe}\,y^{2}}{2{\cal
D}x}}\right)\bigg{]}dydx\;.$ (15)

Integrate $y$ by parts, and make the change of variable $u=x/R$ to simply the
integral:

$\tilde{\chi}^{s}_{\text{mf}}=-\overline{\rho}H_{\rho,\rho}\sqrt{\text{Pe}/R}\\\
\times\int^{1}_{0}\frac{\sqrt{1-u^{2}}}{u^{3/2}}\exp\left[{\frac{\text{Pe}\,R}{2{\cal
D}}\left(u-\frac{1}{u}\right)}\right]du\;.$ (16)

This final integral is easily computed numerically. Moreover, the argument of
the exponential is always negative so one easily sees that
$\tilde{\chi}^{s}_{\text{mf}}\to 0$ as the cutoff $R\to\infty$. Figure 11(b)
shows the $R$-dependence of results from simulation, showing that
$\tilde{\chi}$ increases weakly with $R$, contrary to this prediction.
(Agreement is not expected here, because the mean-field theory does not
capture large density fluctuations such as TJCs.)

Figure 12: (a) Dependence of $h_{\rho,\rho}(\bm{r})$ on $\bm{r}=(x,y)$. (b)
The density-density correlation function $h_{\rho,\rho}(0,y)$ and
$h_{\rho,\rho}(x,0)$. The inset is the difference
$h_{\rho,\rho}(0,y)-h_{\rho,\rho}(x,0)$. The sign of
$h_{\rho,\rho}(0,y)-h_{\rho,\rho}(x,0)$ suggests $H_{\rho,\rho}>0$ in the
mean-field theory. For both (a) and (b), $L_{\parallel}=150d_{0}$ and $\rm
Pe=20$.

For a more detailed analysis, Fig. 12 shows the behaviour of $h_{\rho,\rho}$
for ${\rm Pe}=20$. In particular, we show the behaviour along the axes
$\bm{r}=(x,0)$ and $\bm{r}=(0,y)$. For the latter case, Eq. (14) predicts that
the hydrodynamic part of the correlation is exponentially small but the
numerical results show a small but positive contribution that decays slowly at
large scales. (We attribute this to TJCs that are extended along the
$y$-direction, which are not accounted for by mean-field theory.) Since
$h_{\rho,\rho}(x,0)<0$ at large distances while $h_{\rho,\rho}(0,y)>0$, the
best chance for agreement with Eq. (14) is to take $H_{\rho,\rho}>0$,
indicating $\tilde{v}_{0}<1.3$. The mean-field results are obtained with
$\tilde{v}_{0}=1$, however, this leads to $\tilde{\chi}^{s}_{\text{mf}}<0$,
see Fig. 11(b) (recall also that this $\tilde{\chi}^{s}_{\text{mf}}\to 0$ for
large $R$).

##  Appendix D Evolution of Interface

In this section we describe how we fit the time-evolution of the interface
within demixed state described in the main text. We provide evidence that the
non-diffusive behavior of the time-evolution of the interface at the high Pe
regime depends weakly on the shape of the simulation box.

Figure 13: The charge density profiles for systems with different aspect
ratios near the interface at $t=25\tau_{B},75\tau_{B},125\tau_{B}$ for
$\text{Pe}=60$ with $L_{\parallel}=100d_{0}$.  Figure 14: The charge density
profiles with different aspect ratios, varying Pe with
$L_{\parallel}=100d_{0}$.

###  D.1 Fitting for diffusive interface profile

For small Pe, we expect the mixing for domains of red and blue particles to be
described by the diffusion equation:

$\partial_{t}\langle\sigma(y,t)\rangle=D_{y}\partial_{y}^{2}\langle\sigma(y,t)\rangle,$
(17)

where $D_{y}$ is a lateral diffusion constant. The initial condition is (for
$0\leq y<L_{y}$)

$\langle\sigma(y,0)\rangle=\begin{cases}1&0\leq y<L_{\perp}/2\\\
-1&L_{\perp}/2\leq y<L_{\perp}\\\ \end{cases},$ (18)

Since we consider a system with periodic boundary conditions, we extend the
function $\langle\sigma(y,t)\rangle$ to include all values of $y$ and we seek
periodic solutions, such that
$\langle\sigma(y+L_{\perp},t)\rangle=\langle\sigma(y,t)\rangle$ for all $y,t$.
The corresponding solution to (17) can be represented as

$\langle\sigma(y,t)\rangle=\sum_{k=-\infty}^{\infty}\left[\operatorname{erf}\left(\frac{y-kL_{\perp}}{\sqrt{4D_{y}t}}\right)-\operatorname{erf}\left(\frac{y-\frac{2k-1}{2}L_{\perp}}{\sqrt{4D_{y}t}}\right)\right]$
(19)

For $0<y<L_{\perp}$ and moderate values of $t$, the contributions from large
values of $k$ are negligible, so it is practical to truncate the sums. In
fact, for the relatively short times considered here, it is sufficient to
approximate

$\langle\sigma(y,t)\rangle=\operatorname{erf}\left(\frac{y}{\sqrt{4D_{y}t}}\right)-\operatorname{erf}\left(\frac{y-\tfrac{1}{2}L_{\perp}}{\sqrt{4D_{y}t}}\right)\\\
+\operatorname{erf}\left(\frac{y-L_{\perp}}{\sqrt{4D_{y}t}}\right).$ (20)

since all other contributions are small for $0\leq y<L_{\perp}$. This form
gives the fits shown in Fig. 8.

###  D.2 Interface mixing at different box shapes

We explain in the main text that TJCs form at the interface between demixed
red and blue domains, which causes the demixed state to break down in large
systems. As additional supporting evidence for this mechanism, we simulated
systems with different aspect ratios: $L_{\parallel}=100d_{0}$ and
$L_{\perp}=40d_{0},60d_{0},80d_{0},100d_{0}$ at $\text{Pe}=50,60,70$. This is
the regime where diffusive mixing does not apply, and TJCs play an important
role.

In Fig. 13 we show the time-evolution of the interface profile at
$\text{Pe}=60$ and the mechanism described is consistent at different time
steps. Each interface profile is averaged over at least 20 simulations and the
error bar for the density profile is taken to be the standard error of the
average.

The mixing effect depends very weakly on $L_{\perp}$, except for the smallest
value ($L_{\perp}=40d_{0}$) which is small enough that a single TJC can affect
the whole domain, and is no longer localised near an interface. For the larger
systems, this indicates that the mixing can be interpreted as an instability
of a single interface between red and blue, initiated by a TJC.

Finally, Fig. 14 shows snapshots of the interface profile at
$\text{Pe}=50,60,70$: the non-diffusive mixing mechanism is consistently
observed across this range.

## References

  * Couzin and Krause [2003] I. D. Couzin and J. Krause, Adv. Study Behav. 32, 10 (2003).
  * Ballerini  _et al._ [2008] M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, A. Orlandi, G. Parisi, A. Procaccini, M. Viale, _et al._ , Animal Behaviour 76, 201 (2008).
  * Bialek  _et al._ [2012] W. Bialek, A. Cavagna, I. Giardina, T. Mora, E. Silvestri, M. Viale, and A. M. Walczak, Proc. Natl. Acad. Sci. USA 109, 4786 (2012).
  * Sanchez  _et al._ [2012] T. Sanchez, D. T. Chen, S. J. DeCamp, M. Heymann, and Z. Dogic, Nature 491, 431 (2012).
  * Doostmohammadi  _et al._ [2018] A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans, and F. Sagués, Nat. Commun. 9, 3246 (2018).
  * Tailleur and Cates [2008] J. Tailleur and M. Cates, Phys. Rev. Lett. 100, 218103 (2008).
  * Thompson  _et al._ [2011] A. G. Thompson, J. Tailleur, M. E. Cates, and R. A. Blythe, J. Stat. Mech.: Theory Exp. 2011 (02), P02029.
  * Cates and Tailleur [2015] M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys. 6, 219 (2015).
  * Turing [1952] A. M. Turing, Proc. R. Soc. London B 237, 37 (1952).
  * Gray and Scott [1983] P. Gray and S. Scott, Chem. Eng. Sci. 38, 29 (1983).
  * Burger  _et al._ [2016] M. Burger, S. Hittmeir, H. Ranetbauer, and M.-T. Wolfram, SIAM J. Math. Anal. 48, 981 (2016).
  * Dzubiella  _et al._ [2002] J. Dzubiella, G. Hoffmann, and H. Löwen, Phys. Rev. E 65, 021402 (2002).
  * Leunissen  _et al._ [2005] M. E. Leunissen, C. G. Christova, A.-P. Hynninen, C. P. Royall, A. I. Campbell, A. Imhof, M. Dijkstra, R. Van Roij, and A. Van Blaaderen, Nature 437, 235 (2005).
  * Vissers  _et al._ [2011a] T. Vissers, A. Wysocki, M. Rex, H. Löwen, C. P. Royall, A. Imhof, and A. van Blaaderen, Soft Matter 7, 2352 (2011a).
  * Sütterlin  _et al._ [2009] K. Sütterlin, A. Wysocki, A. Ivlev, C. Räth, H. Thomas, M. Rubin-Zuzic, W. Goedheer, V. Fortov, A. Lipaev, V. Molotkov, _et al._ , Phys. Rev. Lett. 102, 085003 (2009).
  * Sutterlin  _et al._ [2009] K. R. Sutterlin, H. M. Thomas, A. V. Ivlev, G. E. Morfill, V. E. Fortov, A. M. Lipaev, V. I. Molotkov, O. F. Petrov, A. Wysocki, and H. Lowen, IEEE Trans. Plasma Sci. 38, 861 (2009).
  * Du  _et al._ [2012] C. Du, K. Sütterlin, K. Jiang, C. Räth, A. Ivlev, S. Khrapak, M. Schwabe, H. Thomas, V. Fortov, A. Lipaev, _et al._ , New J. Phys. 14, 073058 (2012).
  * Sarma  _et al._ [2020] U. Sarma, S. Baruah, and R. Ganesh, Phys. Plasmas 27, 012106 (2020).
  * Lee  _et al._ [2014] A. A. Lee, S. Kondrat, G. Oshanin, and A. A. Kornyshev, Nanotechnology 25, 315401 (2014).
  * Helbing and Molnar [1995] D. Helbing and P. Molnar, Phys. Rev. E 51, 4282 (1995).
  * Isobe  _et al._ [2004] M. Isobe, T. Adachi, and T. Nagatani, Physica A 336, 638 (2004).
  * Karamouzas  _et al._ [2014] I. Karamouzas, B. Skinner, and S. J. Guy, Phys. Rev. Lett. 113, 238701 (2014).
  * Oliveira  _et al._ [2016] C. L. Oliveira, A. P. Vieira, D. Helbing, J. S. Andrade Jr, and H. J. Herrmann, Phys. Rev. X 6, 011003 (2016).
  * Bacik  _et al._ [2023] K. A. Bacik, B. S. Bacik, and T. Rogers, Science 379, 923 (2023).
  * Couzin and Franks [2003] I. D. Couzin and N. R. Franks, Proc. R. Soc. London B 270, 139 (2003).
  * Wysocki and Löwen [2009] A. Wysocki and H. Löwen, Phys. Rev. E 79, 041408 (2009).
  * Vissers  _et al._ [2011b] T. Vissers, A. van Blaaderen, and A. Imhof, Phys. Rev. Lett. 106, 228303 (2011b).
  * Li  _et al._ [2021] B. Li, Y.-L. Wang, G. Shi, Y. Gao, X. Shi, C. E. Woodward, and J. Forsman, ACS Nano 15, 2363 (2021).
  * Rex and Löwen [2008] M. Rex and H. Löwen, The European Physical Journal E 26, 143 (2008).
  * Glanz and Löwen [2012] T. Glanz and H. Löwen, J. Phys.: Condens. Matter 24, 464114 (2012).
  * Kohl  _et al._ [2012] M. Kohl, A. V. Ivlev, P. Brandt, G. E. Morfill, and H. Löwen, J. Phys.: Condens. Matter 24, 464115 (2012).
  * Klymko  _et al._ [2016] K. Klymko, P. L. Geissler, and S. Whitelam, Phys. Rev. E 94, 022608 (2016).
  * Poncet  _et al._ [2017] A. Poncet, O. Bénichou, V. Démery, and G. Oshanin, Phys. Rev. Lett. 118, 118002 (2017).
  * Chakrabarti  _et al._ [2003] J. Chakrabarti, J. Dzubiella, and H. Löwen, Europhys. Lett. 61, 415 (2003).
  * Chakrabarti  _et al._ [2004] J. Chakrabarti, J. Dzubiella, and H. Löwen, Phys. Rev. E 70, 012401 (2004).
  * Reichhardt and Reichhardt [2006] C. Reichhardt and C. J. O. Reichhardt, Phys. Rev. E 74, 011403 (2006).
  * Vasilyev  _et al._ [2017] O. A. Vasilyev, O. Bénichou, C. Mejía-Monasterio, E. R. Weeks, and G. Oshanin, Soft Matter 13, 7617 (2017).
  * Schimansky-Geier  _et al._ [2021] L. Schimansky-Geier, B. Lindner, S. Milster, and A. B. Neiman, Phys. Rev. E 103, 022113 (2021).
  * Yu  _et al._ [2022] H. Yu, K. Thijssen, and R. L. Jack, Phys. Rev. E 106, 024129 (2022).
  * VanSaders and Vitelli [2023] B. VanSaders and V. Vitelli, arXiv preprint arXiv:2302.07402 (2023).
  * Geigenfeind  _et al._ [2020] T. Geigenfeind, D. de las Heras, and M. Schmidt, Commun. Phys. 3, 1 (2020).
  * [42] see supplemental material for the supplementary movies of the steady state behaviors at different peclet number regimes and the mixing between the counterpropagating red and blue domains triggered by tjcs.
  * Weeks  _et al._ [1971] J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chem. Phys. 54, 5237 (1971).
  * Thompson  _et al._ [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, Comp. Phys. Comm. 271, 108171 (2022).
  * Korniss  _et al._ [1995] G. Korniss, B. Schmittmann, and R. Zia, Europhys. Lett. 32, 49 (1995).
  * Korniss  _et al._ [1997] G. Korniss, B. Schmittmann, and R. Zia, J. Stat. Phys. 86, 721 (1997).
  * Démery and Dean [2016] V. Démery and D. S. Dean, J. Stat. Mech.: Theory Exp. 2016 (2), 023106.
  * Frusawa [2022] H. Frusawa, Entropy 24, 500 (2022).
  * Schmittmann and Zia [1995] B. Schmittmann and R. Zia, in  _Statistical Mechanics of Driven Diffusive System_ , Phase Transitions and Critical Phenomena, Vol. 17, edited by B. Schmittmann and R. Zia (Academic Press, 1995) pp. 3–214.
  * Hansen and McDonald [2013] J.-P. Hansen and I. R. McDonald,  _Theory of simple liquids: with applications to soft matter_ (Academic press, 2013).
  * Crooks and Chandler [1997] G. E. Crooks and D. Chandler, Phys. Rev. E 56, 4217 (1997).
  * Del Junco  _et al._ [2018] C. Del Junco, L. Tociu, and S. Vaikuntanathan, Proc. Natl. Acad. Sci. USA 115, 3569 (2018).
  * Omar  _et al._ [2021] A. K. Omar, K. Klymko, T. GrandPre, and P. L. Geissler, Phys. Rev. Lett. 126, 188002 (2021).
  * Dutta and Chakrabarti [2016] S. Dutta and J. Chakrabarti, Europhys. Lett. 116, 38001 (2016).
  * Dutta and Chakrabarti [2018] S. Dutta and J. Chakrabarti, Soft Matter 14, 4477 (2018).
  * Dutta and Chakrabarti [2020] S. Dutta and J. Chakrabarti, Phys. Chem. Chem. Phys. 22, 17731 (2020).
  * Del Junco and Vaikuntanathan [2019] C. Del Junco and S. Vaikuntanathan, J. Chem. Phys. 150, 094708 (2019).
  * Dean  _et al._ [2020] D. S. Dean, P. Gersberg, and P. C. Holdsworth, J. Stat. Mech.: Theory Exp. 2020 (3), 033206.
  * Bouchet  _et al._ [2019] F. Bouchet, J. Rolland, and E. Simonnet, Phys. Rev. Lett. 122, 074502 (2019).
  * Helbing  _et al._ [2000] D. Helbing, I. J. Farkas, and T. Vicsek, Phys. Rev. Lett. 84, 1240 (2000).
  * Wensink and Löwen [2012] H. Wensink and H. Löwen, J. Phys.: Condens. Matter 24, 464130 (2012).
  * Farrell  _et al._ [2012] F. D. C. Farrell, M. C. Marchetti, D. Marenduzzo, and J. Tailleur, Phys. Rev. Lett. 108, 248101 (2012).
  * Bain and Bartolo [2017] N. Bain and D. Bartolo, Nat. Commun. 8, 15969 (2017).
  * Bär  _et al._ [2020] M. Bär, R. Großmann, S. Heidenreich, and F. Peruani, Annu. Rev. Condens. Matter Phys. 11, 441 (2020).
  * Stopper and Roth [2018] D. Stopper and R. Roth, Phys. Rev. E 97, 062602 (2018).
  * Solon  _et al._ [2022] A. Solon, H. Chaté, J. Toner, and J. Tailleur, Phys. Rev. Lett. 128, 208004 (2022).
  * Anderson  _et al._ [2023] C. Anderson, G. Goldsztein, and A. Fernandez-Nieves, Sci. Adv. 9, eadd0635 (2023).
  * Zhang and Fodor [2023] Y. Zhang and É. Fodor, Phys. Rev. Lett. 131, 238302 (2023).
  * [69] Y.-E. Keta, J. U. Klamser, R. L. Jack, and L. Berthier, ArXiv:2306.07172.
  * Chennakesavalu and Rotskoff [2021] S. Chennakesavalu and G. M. Rotskoff, J. Chem. Phys. 155, 194114 (2021).
  * Berthier  _et al._ [2019] L. Berthier, E. Flenner, and G. Szamel, J. Chem. Phys. 150, 200901 (2019).
  * Keta  _et al._ [2022] Y.-E. Keta, R. L. Jack, and L. Berthier, Phys. Rev. Lett. 129, 048002 (2022).