Particle-conserving dynamics on the single-particle level

TL;DR

This paper extends the particle-conserving dynamics method to binary mixtures, applying it to one-dimensional hard rods, and compares its accuracy with Brownian dynamics and density functional theory, highlighting its strengths and limitations.

Contribution

The authors generalize particle-conserving dynamics to binary mixtures and analyze its effectiveness for one-dimensional hard rods, especially in tagged-particle dynamics.

Findings

Particle-conserving dynamics outperforms dynamical density functional theory at short and intermediate times.

The method accurately reproduces simulation data but has limitations at long times due to neglect of particle order.

Fundamental limitations of density-based approaches are highlighted in systems with particle caging.

Abstract

We generalize the particle-conserving dynamics method of de las Heras et al. [J. Phys. Condens. Matter: 28, 24404 (2016).] to binary mixtures and apply this to hard rods in one dimension. Considering the case of one species consisting of only one particle enables us to address the tagged-particle dynamics. The time-evolution of the species-labeled density profiles is compared to exact Brownian dynamics and (grand-canonical) dynamical density functional theory. The particle conserving dynamics yields improved results over the dynamical density functional theory and well reproduces the simulation data at short and intermediate times. However, the neglect of a strict particle order (due to the fundamental statistical assumption of ergodicity) leads to errors at long times for our one-dimensional setup. The isolated study of that error makes clear the fundamental limitations of (adiabatic) density-based theoretical approaches when applied to systems of any dimension for which particle caging is a dominant physical mechanism.

Equations

$$\Xi ({\mu }_1,{\mu }_2)=\sum _{N_1=0}^{\infty }\sum _{N_2=0}^{\infty }{e}^{\beta ({\mu }_1{N}_1+{\mu }_2{N}_2)}{Z}_{N_1{N}_2}.$$

$${\rho }^{(\nu )}(x)=\sum _{N_1=0}^{\infty }\sum _{N_2=0}^{\infty }{p}_{N_1{N}_2}({\mu }_1,{\mu }_2){\rho }_{N_1{N}_2}^{(\nu )}(x),$$

$$p_N({\mu }_1,{\mu }_2)=e^{\beta ({\mu }_1{N}_1+{\mu }_2{N}_2)}\frac{Z_{N_1{N}_2}}{\Xi ({\mu }_1,{\mu }_2)}.$$

$${\rho }_{N_1{N}_2}^{(\nu )}(x)=\sum _{({\mu }_1,{\mu }_2)}^M{\mathbf{P}}_{N_1{N}_2}({\mu }_1,{\mu }_2{)}^{-1}{\rho }^{(\nu )}(x).$$

$$\frac{\delta \Omega [\{{\rho }^{({\nu }^{\prime })}\}]}{\delta {\rho }^{(\nu )}(x)}=0.$$

$$\Omega [\{{\rho }^{(\nu )}\}]={\mathcal{F}}_{\text{id}}+{\mathcal{F}}_{\text{ex}}+\sum _{\nu =1}^2\int \mathrm{d}x{\rho }^{(\nu )}(x)\left(V_{\text{ext}}^{(\nu )}(x)-{\mu }_{\nu }\right).$$

$$\beta {\mathcal{F}}_{\text{id}}[\{{\rho }^{(\nu )}\}]=\sum _{\nu =1}^2\int \mathrm{d}x{\rho }^{(\nu )}(x)(\ln (\Lambda {\rho }^{(\nu )}(x))-1)$$

$$\beta {\mathcal{F}}_{\text{ex}}[\rho ]=-\int \mathrm{d}x{n}_0(x)\ln \left(1-n_1(x)\right).$$

$$\displaystyle u(x)=\left\{\begin{array}[]{c}0\ \ \ \ \,\text{if}\ |x|>\sigma\\ \infty\ \ \ \text{else}\hfill\end{array}\right.$$

$$n_i(x)$$

$${\omega }^{(0)}(x)$$

$${\omega }^{(1)}(x)$$

$$F_N^{\text{tot}}[\{{\rho }_N^{(\nu )}\}]=F_N[\{{\rho }_N^{(\nu )}\}]+\sum _{\nu =1}^2\int \mathrm{d}x{\rho }_N^{(\nu )}(x)V_{\text{ext}}^{(\nu )}(x)$$

$$\beta F_N[\{{\rho }_N^{(\nu )}\}]=-\ln Z_N-\sum _{\nu =1}^2\int \mathrm{d}x{\rho }_N^{(\nu )}(x)\beta V^{(\nu )}(x).$$

$$\beta V_n^{(\nu )}(x)=\beta V_{n-1}^{(\nu )}(x)-\ln {\rho }_N^{(\nu )}(x)+\ln {\rho }_{n-1}^{(\nu )}(x).$$

$$V_{\text{ext}}^{(\nu )}(x)-V^{(\nu )}(x)-\sum _{{\nu }^{\prime }=1}^2\int \mathrm{d}{x}^{\prime }{\rho }_N^{(\nu )}(x^{\prime })\frac{\delta V^{({\nu }^{\prime })}(x^{\prime })}{\delta {\rho }_N^{(\nu )}(x)}=0$$

$$Z_2=\frac{1}{2!{\Lambda }^2}{\int }_R^{L-R}\mathrm{d}{x}_1{\int }_R^{L-R}\mathrm{d}{x}_2{e}^{-\beta u(|x_2-x_1|)}$$

$$Z_{11}\equiv \frac{1}{{\Lambda }^2}{\int }_R^{L-3R}\mathrm{d}{x}_1{\int }_{x_1+2R}^{L-R}\mathrm{d}{x}_2\equiv {\mathcal{Z}}_{11}[1,2],$$

$$Z_{11}=\frac{1}{{\Lambda }^2}{\int }_R^{L-R}\mathrm{d}{x}_1{\int }_R^{L-R}\mathrm{d}{x}_2{e}^{-\beta o(x_2-x_1)},$$

$$\displaystyle o(x)=\left\{\begin{array}[]{c}0\ \ \ \ \,\text{if}\ x>\sigma\\ \infty\ \ \ \text{else}\hfill\end{array}\right.\,$$

$$Z_2=\frac{1}{2!}\left({\mathcal{Z}}_{11}[1,2]+{\mathcal{Z}}_{11}[2,1]\right),$$

$${\rho }_{11}^{(1)}(x)$$

$${\rho }_{11}^{(2)}(x)$$

$${\rho }_2^{(1)}(x)={\rho }_2^{(2)}(x)=\frac{1}{2}\left({\rho }_{11}^{(1)}(x)+{\rho }_{11}^{(2)}(x)\right)\equiv \frac{1}{2}{\rho }_2(x),$$

$$\frac{\partial {\rho }^{(\nu )}(x,t)}{\partial t}=D_0\frac{\partial }{\partial x}$$

$$+{\rho }^{(\nu )}(x,t)\frac{\partial \beta V_{\text{ext}}^{(\nu )}(x)}{\partial x}),\text{mbox}$$

$${\mathbf{f}}_{\text{ad}}^{(\nu )}(x,t)=-{\rho }^{(\nu )}(x,t)\frac{\partial }{\partial x}\frac{\delta {\mathcal{F}}_{\text{ex}}[\{{\rho }^{({\nu }^{\prime })}\}]}{\delta {\rho }^{(\nu )}(x,t)}$$

$$\beta {\mathbf{f}}_{N,\text{ad}}^{(\nu )}(x,t)={\rho }_N^{(\nu )}(x,t)\frac{\partial }{\partial x}\left(\beta V_{\text{ad}}^{(\nu )}(x,t)+\ln {\rho }_N^{(\nu )}(x,t)\right).$$

$$\frac{\partial {\rho }_N^{(\nu )}(x,t)}{\partial t}=\beta D_0\frac{\partial }{\partial x}\left({\rho }_N^{(\nu )}(x,t)\frac{\partial \left(V_{\text{ext}}^{(\nu )}(x)-V_{\text{ad}}^{(\nu )}(x,t)\right)}{\partial x}\right)$$

$${⟨f(x)⟩}_{\nu }\equiv \frac{1}{N_{\nu }}\int \mathrm{d}x{\rho }^{(\nu )}(x)f(x)$$

$$\beta {\mathcal{F}}_{\text{id}}[\rho ]=\int \mathrm{d}x\rho (x)(\ln [\Lambda \rho (x)]-1),$$

$${\mathbf{f}}_i(t)=-\frac{\partial }{\partial x_i}{V}_{\text{ext},0}^{(\nu )}(x_i)\Theta (-t)+R_i(t),$$

Full text

I Introduction

Dynamical density functional theory (DDFT) is a widely used and versatile tool for investigating the dynamics of bulk and inhomogeneous classical systems of interacting Brownian particles. By assuming that all pair and higher-order correlation functions equilibrate much faster than the one-body (or one-point) density (an adiabatic approximation) the DDFT exploits the formally exact statistical mechanical method of density functional theory (DFT) to approximately treat nonequilibrium situations marconi_tarazona .

DDFT has been applied with much success to study a variety of important physical phenomena, e.g. spinodal decomposition archer_evans , colloidal sedimentation royall and quasicrystal formation quasi . When combined with the test-particle method, whereby one of the particles is treated as an external field, DDFT can be used to calculate the self and distinct parts of the van Hove function dtp ; stopper_bulk and thus address the dynamics of equilibrium states. Extensions of the theory to treat driven systems can reproduce the phenomenology of colloidal system under external shear flow bk1 ; scacchi ; stopper_shear . More recently a general and exact variational framework, power functional theory, has been developed to treat nonequilibrium Brownian systems pft1 . Within this approach, DFT and DDFT emerge as equilibrium and adiabatic limits, respectively. Moreover, power functional theory enables the superadiabatic contributions to the dynamics to be approximated in a physically intuitive way pft2 .

An implicit drawback to standard implementations of DDFT is that interaction forces are generated from a grand-canonical free energy functional. For confined systems with small numbers of particles or, perhaps more generally, systems which exhibit strong density inhomogeneities, grand-canonical DFT can predict density profiles which differ significantly from canonical simulations at fixed particle number $N$. This issue was addressed in the late 1990’s by González et al. canonical1 ; canonical2 , who employed an expansion in inverse powers of the particle number to systematically approximate canonical density profiles using grand-canonical information as an input. The problem of calculating equilibrium density profiles in the canonical ensemble was revisited in 2014 by de las Heras and Schmidt canonical3 who showed how to obtain exact canonical information from grand-canonical DFT given an exact functional, performing explicit calculations for the one-dimensional hard-rod system. The method of de las Heras and Schmidt was further generalized in Ref. [canonical4, ] to generate a theory of particle-conserving dynamics (PCD) for the time-evolution of the one-point density of $N$ particles. This theory, while still relying on the adiabatic approximation, eliminates spurious effects arising from the grand-canonical ensemble and yields predictions for the time-dependent density in good agreement with Brownian dynamics (BD) simulation data.

In this paper we provide an intuitive generalization of the PCD approach of Ref. [canonical4, ] to binary mixtures and use this to examine the dynamics of a tagged particle in a one-component system. By tagging a particle we can investigate the physics of dynamical confinement or “caging” within the framework of an adiabatic time-evolution equation for the one-point density. For clarity in our terminology we only speak here of a “one-point density”, as opposed to two- or $n$-point densities (correlations), instead of using the equivalent term “one-body density”, which should not be confused with the particular (canonical) density for $N=1$ particle, i.e., a single-particle density profile. Our calculations are performed for hard-rods in one spatial dimension because (i) the exact grand-canonical density functional is known and (ii) the Tonks gas tonks provides an extreme case of a non-ergodic fluid, which allows us to illustrate most clearly the physics of localization and caging. Finally, (iii), the dynamic state of hard rods constitutes a fundamental model for single-file diffusion of particles that never swap positions, which can be solved in a closed form LA1 ; LA2 but not in the context of a variational framework.

The paper is organized as follows. In Sec. II we describe the transformation from grand-canonical to canonical information for mixtures and discuss the statistical background of hard rods in one dimension. Then we embed in Sec. III the canonical free energy functional in a dynamical framework by making an adiabatic approximation of DDFT. Then we apply this PCD approach to the relaxation dynamics of single-particle profiles of hard rods and discuss in detail the similarities and differences to BD results. Our findings are discussed in Sec. IV.

II Variational calculus for canonical mixtures

The variational character of DFT is fundamental to its usefulness in addressing the physics of the liquid state. However, for this to apply it is necessary to work in the grand-canonical ensemble. We thus adapt the recent method developed by de las Heras and Schmidt to use the grand-canonical data obtained from DFT, specifically the one-point density profiles and partition functions, to calculate the one-point density profile of one-component canonical systems, for which the particle number does not fluctuate canonical4 . Here we show that it is straightforward to generalize the grand-canonical–canonical transformation method to binary (or $\kappa$-component) mixtures.

II.1 Canonical information from a grand-canonical theory

Consider a two-component system containing $N_{\nu}$ particles of species $\nu\in\{1,2\}$. Suppose we have total grand-canonical information on such a system, i.e., for the given chemical potentials $\mu_{\nu}$, we know the grand partition function $\Xi(\mu_{1},\mu_{2})$ and the grand-canonical one-point density profiles $\rho^{(\nu)}(x)$ of each species $\nu$ at position $x$. Formally, these can be obtained by an infinite summation of the canonical partition function $Z_{N_{1}N_{2}}$ and one-point density $\rho_{N_{1}N_{2}}^{(\nu)}(x)$ according to

[TABLE]

and

[TABLE]

where $\beta=(k_{\text{B}}T)^{-1}$ denotes the inverse temperature with Boltzmann’s constant $k_{\text{B}}$ and the probability $p_{N_{1}N_{2}}(\mu_{1},\mu_{2})$ to find a state with $N_{1}$ and $N_{2}$ particles at a given pair of chemical potentials $(\mu_{1},\mu_{2})$ is given by

[TABLE]

In practice, the sums can be truncated at the maximum particle number $N_{\text{max}}$ in finite systems, i.e., all partition sums and probabilities with $N_{1}+N_{2}>N_{\text{max}}$ vanish. Otherwise, there is no difference between the two ensembles in the thermodynamic limit.

We calculate the grand-canonical density distributions and partition functions to set up a system of $M\equiv(N_{\text{max}}+1)(N_{\text{max}}+2)/2$ linear equations in the form of Eq. (1), with $M$ being the number of possible pairs of particle numbers, $(N_{1},N_{2})$, with $N_{1}+N_{2}\leq N_{\text{max}}$. We then solve this equation system for $Z_{N_{1}N_{2}}$ and calculate the probabilities $p_{N_{1}N_{2}}(\mu_{1},\mu_{2})$ via Eq. (3). By interpreting these probabilities as entries of an $M\times M$ matrix $\mathbf{P}_{N_{1}N_{2}}(\mu_{1},\mu_{2})$ and inverting this matrix we finally obtain the canonical densities via

[TABLE]

The sum on the right-hand side includes $M$ arbitrary pairs of chemical potentials $\mu_{1}$ and $\mu_{2}$ of the two components. The explicit choice does not matter given the available grand-canonical information is exact. For the sake of numerical robustness, a certain range of resulting average particle numbers should be covered. As detailed in Ref. [canonical3, ], we could reduce the number of equations by one upon removing the trivial case of zero particles in each species, which has been omitted in the present study due to the negligible effect on computation time.

The above methods can be easily generalized to systems containing any number $\kappa$ of different components, where the number of coupled linear equations to be solved grows exponentially with $\kappa$.

It is important to notice that the density profiles and partition sums depend explicitly on the numbers of particles of each species and not only on the total particle number $N=N_{1}+N_{2}$, as soon as the species are physically distinguishable (by their particle-particle interactions or by interactions with external potentials). However, we will consider systems with fixed values of $N_{1}$ and $N_{2}$ in the remainder of the paper and hence we will indicate quantities in canonical ensembles by $N$ to unclutter the notation. Quantities with two indices $N_{1}$ and $N_{2}$ will then indicate ordered ensembles (see Sec. II.4).

II.2 Classical density functional theory

Our starting point is classical DFT, providing, via a variational formalism, the grand partition function and the grand-canonical equilibrium density profiles of an arbitrary mixture of particles under the influence of any external field $V_{\text{ext}}^{(\nu)}(x)$ acting on each component $\nu$. Given a density functional $\Omega[\{\rho^{(\nu)}\}]$ of the grand potential, one obtains the grand-canonical equilibrium density profiles of each component from a variational minimization according to

[TABLE]

Substituting the resulting equilibrium density profile into the functional yields the equilibrium grand potential $\Omega$ and, therefore, the grand partition function $\Xi=\exp(-\beta\Omega)$. Canonical information is then accessible via inversion of Eqs. (1) and (2).

Here, we restrict ourselves to hard rods of length $\sigma$ in one spatial dimension, where the exact density functional reads

[TABLE]

with the contribution

[TABLE]

of an ideal non-interacting gas, where $\Lambda$ is the thermal wavelength. The excess free energy functional of hard-rods was derived by Percus percus ; percusmix and is given by

[TABLE]

This functional describes the contribution of the pair interaction potential

[TABLE]

and is a function of the two weighted densities

[TABLE]

where $i$ is not a species label but serves to enumerate the weighted densities. The corresponding weight functions are given by

[TABLE]

where $R=\sigma/2$ is half the length of a rod and $\delta(x)$ and $\Theta(x)$ denote the Dirac distribution and the Heaviside step function, respectively. Here we dropped all species labels, since we only consider identical particles (in general, the functional allows for the description of $\kappa$ actually different species of lengths $2R_{\nu}$). For the virtual mixture considered here, we further use the notation $\mathcal{F}_{\text{ex}}[\rho]$ to emphasize that Eq. (8) then depends only on the total density profile $\rho(x)=\rho^{(1)}(x)+\rho^{(2)}(x)$.

II.3 Canonical intrinsic free energy functional

Although there exists an explicit expression for a density functional $\mathcal{F}^{\text{tot}}=\Omega+\sum_{\nu}\int\mathrm{d}x\,\rho^{(\nu)}(x)\,\mu_{\nu}$ of the total Helmholtz free energy, cf. Eq. (6), its minimization under the constraint of fixed particle numbers $N_{\nu}=\int\mathrm{d}x\rho^{(\nu)}(x)$ of each component would still result in a grand-canonical profile, i.e., a superposition of canonical density profiles with different particle numbers that average to $N_{\nu}$. This constraint is exactly what requires one to introduce chemical potentials. Our objective is rather to provide a true canonical DFT of the form

[TABLE]

where the total Helmholtz free energy is formally written in its natural variables, which implies that all valid canonical “target” profiles ${\rho_{N}^{(\nu)}}(x)$ must integrate to $N_{\nu}$ and the functional must be minimal in canonical equilibrium. To this end we must perform an iterational search for the intrinsic Helmholtz free energy functional $F_{N}$ on the right-hand side of Eq. (15).

In generalization of one-component case described in Ref. [canonical4, ] we determine for a given pair ${\rho_{N}^{(\nu)}}(x)$ of “target” profiles the corresponding generating external potentials $V^{(1)}(x)$ and $V^{(2)}(x)$ acting on each species such that the target profiles would be equilibrated. Then the intrinsic Helmholtz free energy functional is obtained as

[TABLE]

To determine $V^{(\nu)}(x)$ for a given canonical target profile, we start with an initial guess for $V_{0}^{(\nu)}(x)$ and then employ the gradient-free iteration scheme canonical4

[TABLE]

In each iteration step $n$ we make use of the canonical equilibrium profile $\rho_{n-1}^{(\nu)}(x)$ of species $\nu$ in the external potential $V^{(\nu)}_{n-1}(x)$ of the previous step. It is found by minimizing Eq. (6) with $V^{(\nu)}_{n-1}(x)$ taking the role of $V_{\text{ext}}^{(\nu)}(x)$ and then inverting Eq. (2) as described in Sec. II.1. In practice, we ensure a proper convergence of Eq. (17) by introducing a damping factor close to unity, lowering the weight of the logarithmic terms, and adding a very small number to the density, avoiding a divergence of the logarithm in case of vanishing density. Note that the canonical partition function changes in each iteration step and generally differs from $Z_{N}$ entering in Eq. (16), which is calculated for the true external fields $V_{\text{ext}}^{(\nu)}(x)$ and that the obtained $V^{(\nu)}(x)$ are unique only up to an irrelevant constant related to the initial guess.

In the convention chosen here, the potentials $V^{(\nu)}(x)$ formally replace $V_{\text{ext}}^{(\nu)}(x)$ but still take into account the constraints given by these actual external fields, since, otherwise, Eq. (15) would be ill-defined. Hence, we can separate $V^{(\nu)}(x)$ into $V_{\text{ext}}^{(\nu)}(x)$ and a “nonequilibrium” correction, i.e., due to considering nonequilibrium target profiles ${\rho_{N}^{(\nu)}}(x)$. To illustrate this, we demonstrate how to minimize the canonical free energy functional, Eq. (15), by setting $\delta F^{\text{tot}}_{N}/\delta{\rho_{N}^{(\nu)}}(x)=0$. The resulting condition

[TABLE]

is satisfied by $V_{\text{ext}}^{(\nu)}(x)=V^{(\nu)}(x)$, which can be formally solved for $\rho_{N}^{(\nu)}(x)$ by an iteration similar to Eq. (17). In practice, such a calculation would just amount to determine the canonical equilibrium profile in the external field $V_{\text{ext}}^{(\nu)}(x)$ indirectly by one single grand-canonical-canonical transformation according to Eq. (4).

II.4 Distinguishability in one dimension

Given that the canonical transformation procedure has been demonstrated to be formally correct canonical4 it is interesting to proceed to investigate the limitations of the canonical ensemble for describing systems subject to non-ergodic dynamics. Due to the constraint that the particles remain ordered on the line, the one-dimensional hard rod model presents one of the simplest non-ergodic model systems. We will first highlight the inability of the canonical ensemble to correctly describe species-labeled density profiles in equilibrium, before proceeding in Sec. III to consider the PCD of tagged particles. We will then argue that our findings have strong implications for the ability of any approach based on ensemble-averaged density (adiabatic or superadiabatic) to describe non-ergodic behavior arising from particle caging.

If we employ species labeling simply as a formal device to track either individual particles or subsets of particles, then, within a canonical description, the equilibrium density profiles of a species holding $N_{\nu}$ particles are always given by $\rho_{N}^{(\nu)}=N_{\nu}\rho_{N}(x)/N$, where $\rho_{N}(x)$ is the total canonical equilibrium density profile of $N$ identical particles, irrespective of their species labeling. This is in contradiction with the real situation in systems with broken ergodicity, such as densely packed crystals, glasses or, in the example we lay out in the following, hard rods in one-dimension, where the single-particle profiles should reflect the spatial localization.

To illustrate the origin of localization in a statistical description, we analyze the simplest case with a non-trivial pair interaction: $N=2$ hard rods of length $2R$ confined to a slit of length $L$. The canonical partition function, $Z_{2}=(L-4R)^{2}/(2\Lambda^{2})$, of two particles can be calculated in two different ways. We stress that by the word canonical we always imply the ergodic assumption, i.e., the statistical average implies no particular particle order. The standard approach for completely indistinguishable particles is to calculate

[TABLE]

via the completely symmetric pair interaction potential $u(x)$, specified in Eq. (11). Alternatively, if we distinguish between the two particles and require that particle 1 is always to the left of particle 2, then the ordered partition function reads

[TABLE]

thus implying a broken ergodicity. In the last step we have introduced the formal arguments 1 and 2 to express the explicit dependence on the order of the two particles (the alternative functional notation $\mathcal{Z}_{11}[1,2]$ implies that particle 1 lies to the left of particle 2, which is not decisive for the mathematical value of $Z_{11}$ but serves to indicate the explicit particle order). Equivalently we can calculate

[TABLE]

where the ordering pair potential

[TABLE]

depends on the relative distance and not on its absolute value.

It is straightforward to see that both partition functions are mathematically equal (they evaluate to the same number) and we can formally write

[TABLE]

which is obviously equal to $Z_{11}$. The point we wish to make here is that, although the value of both partition functions $Z_{11}$ and $Z_{2}$ is the same, the corresponding one-point densities of labeled particles, i.e., the ensemble average of the density operators $\hat{\rho}_{\nu}=\delta(x-x_{\nu})$ with respect to the different probability distributions implied in Eq. (19) and Eq. (21), respectively, are different.

The proper calculation for mixtures via the ordered probability distribution, associated with $Z_{11}$, yields the density profiles

[TABLE]

shown as the dotted lines in Fig. 1. The (expected) difference between the two profiles arises from the physical distinction due to the imposed particle order, also present in BD simulations. In contrast, from the canonical distribution, associated with $Z_{2}$, we obtain identical results for each profile

[TABLE]

shown as dashed lines in Fig. 1. In a manner of speaking, we can say that the two species-labeled profiles in Eq. (27) follow from an additional average of those in Eq. (26), accompanied by a loss of microscopic information, whereas their sum equals the exact total canonical profile $\rho_{2}(x)$ for two particles in both cases.

Note that the ordered density profiles $\rho^{(\nu)}_{11}(x)$ need to be understood as correlated averaged quantities, so that the overlap of the profiles does not indicate an explicit crossing of their trajectories . The decreasing probability to find a particle in a specific region rather reflects the likely presence of the other particle according to the predefined order. There are still exclusion strict regions of exactly one particle length on one side for each profile, in addition to the confining wall, This correlation is lost for $\rho^{(\nu)}_{2}(x)$, so that each particle could be found at a position, where the other would not fit in any more. In general, to detect unphysical mixing, we can use the criterion of a nonzero probability for the center of a particle penetrating the minimal region which should remain available for the other particles to respect the predefined order.

The situation described above remains qualitatively the same if the rods are physically distinguishable, e.g., by their lengths. The only difference is the value of the partition functions, since the factor $1/2!$ in $Z_{2}$ has to be removed, so that we find for the only virtual mixture of physically indistinguishable particles $Z_{2}=\left(\mathcal{Z}_{11}[1,2]+\mathcal{Z}_{11}[2,1]\right)$, which equals $2Z_{11}$ in contrast to Eq. (25). This factor does, however, not affect the density profiles generalizing Eqs. (26) and (27) to true mixtures and, therefore, all conclusions drawn and illustrated in the following equally apply to this more general case.

III Particle-conserving dynamics

In this section we use the canonical intrinsic free energy functional, Eq. (16) to drive the dynamics of a mixture of $N=N_{1}+N_{2}$ hard rods in one dimension. In contrast to the standard (grand-canonical) DDFT, this PCD approach operates at fixed numbers of particles of each species, thus providing a more realistic representation of the BD of a system, which, by construction, resolves the positions of all particles at each time $t$ (measured in units of the Boltzmann time $t_{\text{B}}=\sigma^{2}/D_{0}$, where $D_{0}$ is the common diffusion coefficient). Our simulation setup and the averaging process to obtain the species-resolved density profiles $\rho^{(\nu)}_{N_{1}N_{2}}(x,t)$ from BD are described in the Appendix. By considering identical particles and choosing $N_{1}=1$, we can further resolve within PCD the time evolution of the probability density $\rho^{(1)}_{N}(x,t)$ associated with the (average) location of a single particle.

III.1 Adiabatic approximation

The crucial approximation, which allows one to employ any sort of equilibrium DFT in a nonequilibrium framework is to assume that the correlations, i.e., the $n$-particle densities for $n>1$, at each instant of time follow from the time-dependent one-point density in the same way as in equilibrium. This relation is provided by an equilibrium density functional and represents an adiabatic approximation, since we approximate the dynamics as a sequence of equilibrium states.

In DDFT, the one-point density of each species evolves in time according to marconi_tarazona ; archer_evans

[TABLE]

where we do not consider an explicit driving by time-dependent external fields or non-conservative forces. The nonequilibrium interaction force

[TABLE]

is approximately related to the excess free energy functional from Eq. (8). With the above choices, the whole expression in the brackets on the right-hand side of Eq. (29) could be represented, as in Eq. (30), by a functional derivative of the full grand-canonical free energy functional, i.e., a local chemical potential, thereby representing entropic, internal and external forces. We thus describe the adiabatic time evolution of the grand-canonical one-point density. To fix the average particle numbers of each component, we calculate the initial density profiles such that $N_{\nu}=\int\mathrm{d}x\rho^{(\nu)}(x)$ by accordingly choosing the chemical potentials.

Within the PCD framework, we calculate the adiabatic internal force density on species $\nu$ according to canonical4

[TABLE]

These expressions follow intuitively from the (adiabatic) balance with external and entropic forces on each species or, more formally and similar to Eq. (30), from the functional derivative of the canonical excess free energy $\beta F_{N}-\sum_{\nu}\int\mathrm{d}x\rho_{N}^{(\nu)}(\ln\Lambda\rho_{N}^{(\nu)}-1)$ with $F_{N}$ given by Eq. (16). In any case, it is necessary to determine at each time step the adiabatic potentials $V^{(\nu)}_{\text{ad}}(x,t)$ that would generate the instantaneous canonical density profiles $\rho_{N}^{(\nu)}(x,t)$ in equilibrium, as described in Sec. II.3. Having made the adiabatic approximation for the canonical system, we obtain from the analog to Eq. (29) the time evolution equations

[TABLE]

for a mixture in PCD. Here it becomes clear that we can interpret the driving forces of the dynamics as the counterforces to (the nonequilibrium part of) the forces arising from $V^{(\nu)}_{\text{ad}}(x,t)$, providing a physical meaning to the construction of the adiabatic potentials.

When compared to the exact BD, the DDFT approach has, in general, three drawbacks. (i) DDFT conserves only the average number of particles of each species. This has been corrected by our modified PCD approach. Combining the adiabatic canonical profiles according to Eq. (2) with $p_{N}$ independent of time, we could also provide a PCD for a grand-canonical system (different from DDFT canonical4 ), which will not be considered here. Moreover, (ii) superadiabatic forces are always neglected FortiniPRL , and (iii) an inexact canonical functional, obtained through inexact grand-canonical information, results in further deviations which are difficult to quantify. In the latter case, the time evolution will depend on the initial guess $V_{0}^{(\nu)}(x,t)$ in the iteration procedure of Eq. (17), since the vertical offset will formally renormalize the chemical potentials employed in each iteration step. While it is not necessary in the present study using the exact Percus functional, it might be possible to eliminate this problem by properly shifting the generating potentials in each iteration step or, equivalently, adapting the set of chemical potentials. In any case, it is convenient to start with the equilibrium result $V_{0}^{(\nu)}(x,t)=V_{\text{ad}}^{(\nu)}(x,t-\Delta t)$ of the previous time step, which we also do here. In other words, an exact functional, results in exact adiabatic PCD, as prescribed by the Percus functional of the total density profile of hard rods canonical4 .

In the following we use the generalization to mixtures of the Percus functional, which is exact in the grand-canonical sense, i.e., when the number of particles fluctuates and their order does not matter, but, as discussed in Sec. II.4, provides incorrect information on single localized particles. Hence, we expect that PCD for one-dimensional mixtures suffers from both (i) the adiabatic approximation and (ii) the lack of any exact density functional for a system with strict particle order (besides in the trivial case of an ideal gas, where particle trajectories may cross also in BD). To demonstrate the still significant advantages of PCD compared to ordinary DDFT when compared to BD, we discuss in the following the initial relaxation dynamics of hard rods in a slit.

III.2 Relaxation of the one-particle density

To illustrate the performance of the PCD approach for a two-component mixture with a tagged particle, $N_{1}=1$, we study in Figs. 2 and 3 the relaxation of the species-resolved density after bringing the system out of equilibrium by switching off confining external potentials at $t=0$, in which the particles initially were equilibrated with respect to the chosen approach. This initial condition ensures a high degree of localization in the density profiles, which means that, the configurations with the wrong particle order are suppressed in the canonical DFT underlying the PCD. To be more specific, such unphysical states are very unlikely, but not completely eliminated, so that there is already some small, albeit barely noticeable, difference at $t=0$ between PCD and BD with the same fixed particle numbers. We also compare the time evolution to DDFT, which provides grand-canonical states. Even though not physically meaningful we start from the same initial density profiles as for the PCD calculations for better comparison of the dynamics.

In DDFT the profile is interpreted as a grand-canonical one implying repulsion within each species. Thus the time evolution exhibits clear differences even at short times and it becomes apparent that PCD improves significantly over DDFT when comparing to the reference BD.

To learn more about the differences between PCD and BD, we must look a little closer. The PCD and BD profiles for two particles, i.e., $N_{1}=N_{2}=1$, in Fig. 2 are quantitatively nearly indistinguishable for small times. Hence, the neglected superadiabatic forces are insignificant in this case. The increasing deviations emerging at $t\gtrsim 0.6$ are thus of different origin. As the particles can not exchange positions in the one-dimensional trajectory-based BD simulation, the ordering property with particle 1 on the left-hand side of particle 2 is conserved throughout the time evolution. However, we clearly observe that the PCD does not conserve the order of the particles. While the unphysical states with particle 1 on the right-hand side of particle 2 are initially suppressed, the overlap of the density profiles of the two components grows when the system evolves in time. This confirms our expectations based on the model calculations in Sec. II.4 that the present version of DFT is not able to properly drive order-preserving dynamics, and not only PCD. In particular, the PCD profiles will ultimately approach a mixed state without distinction between the components, i.e., the canonical equilibrium profiles depicted in Fig. 1.

We also consider in Fig. 3 the dynamics for a system with $N_{1}=1$ and $N_{2}=2$ particles with the single particle located in the middle of the slit using the same methods. Here, for $N=3$, the role of the superadiabatic forces is a little more obvious than for $N=2$. Regarding the time evolution of species 1, we observe best that the PCD is a little faster than BD, which has also been observed for the one component case canonical4 and can be generally expected from numerical simulations FortiniPRL . At early times, the PCD profiles again show good agreement with the BD results. Moreover, it does not take as long as for two particles until the the onset of the unphysical mixing occurs at $t\gtrsim 0.2$.

To quantify our findings we introduce a general average with respect to the one-point species-resolved densities $\rho^{(\nu)}_{11}(x,t)$ for BD simulations, $\rho^{(\nu)}_{2}(x,t)$ for PCD and $\rho^{(\nu)}(x,t)$ for DDFT, generally defined by

[TABLE]

for any test function $f(x)$. In particular, we consider the values of the density $\rho_{\mathrm{m}}^{(\nu)}(t)\equiv\rho^{(\nu)}(\left\langle x\right\rangle_{\nu},t)$ evaluated at the mean $x$-coordinate $\left\langle x\right\rangle_{\nu}$ and the variance $\text{Var}^{(\nu)}(t)\equiv\left\langle(x-\left\langle x\right\rangle_{\nu})^{2}\right\rangle_{\nu}$ as functions of time.

The corresponding results $\rho_{\mathrm{m}}^{(1)}$ and $\text{Var}^{(1)}(t)$ for species $\nu=1$ are shown in Figs. 4 and 5 for the cases $N=2$ and $N=3$, respectively, as discussed before. In general, all observations on the different dynamics are qualitatively the same for the two systems. The values of DDFT differ significantly from the BD simulation and PCD, which shows that differences between canonical and grand-canonical ensembles are important in such a small system and that the PCD provides a reliable description of the early relaxation dynamics. Here, we also observe the one apparent difference between Figs. 4 and 5. For $N=3$, the PCD has both a clearly lower peak density and a greater variance than BD already at $t=0$. This reflects the importance of choosing some convenient initial conditions, which are less restrictive here than for $N=2$, since the two particles of species $\nu=2$ are not confined in a trap (cf. the caption of Fig. 3) so that they can more easily interchange with the particle in between. This is also the reason why the unphysical mixing becomes apparent at earlier times in Fig. 2 than in Fig. 3.

For larger times, the studied average quantities in the BD simulations eventually reach a plateau, whereas both PCD and DDFT continue to evolve, although the BD are initially slower due to superadiabatic forces. In particular, the value of $\rho_{\mathrm{m}}^{(1)}(t)$ decreases further while $\text{Var}^{(1)}(t)$ continues to grow. This illustrates that the density distributions become wider and flatter than in BD, which is consistent with the mixing of particles. It is clear that the overall timescale for spreading over the whole available system is larger than for the more localized particles in BD.

III.3 An inverse mixing paradox for the PCD of two hard rods

To better understand the underlying mechanism which leads to the unphysical mixing in the PCD and to emphasize that it does not solely occur as a consequence of having chosen improper, already slightly mixed, initial conditions deviating from those of BD, we consider now a more extreme example. We initiate the system in the true equilibrium state $\rho_{11}^{(\nu)}(x)$, as given by Eq. (26) and shown in Fig. 1, respecting the order of two identical hard rods in a slit. Of course, in the BD case, the density profiles will not change over time, since they are already equilibrated.

For the PCD, let us first consider the one-component case, where the (grand-canonical) density functional depends only on the total density $\rho(x)=\rho^{(1)}(x)+\rho^{(2)}(x)$ of two species and different components can be introduced on a formal level. In this case, the ideal free energy functional from Eq. (7) has to be replaced with

[TABLE]

whereas the excess free energy, Eq. (8) depends on the total density only, even for a mixture of identical particles. Using this functional as the basis for the PCD, it does not matter, at any time, which particle we associate with species 1 or 2. So, formally speaking, each pair of grand-canonical (canonical) density profiles for two species (particles), which sum up to the total density are valid adiabatic or equilibrium distributions. In particular, there is no difference between imposing either of the two pairs of one-particle profiles in Eq. (26) or Eq. (27), corresponding to the true canonical equilibrium result of two distinguishable particles or the equilibrium result of the canonical DFT from Sec. II.3, respectively. In both cases, the system remains in equilibrium and the PCD approach is correct.

Now we return to the description of mixtures of identical particles, where the ideal contribution to the free energy is given by Eq. (7). This functional depends explicitly on the individual density profiles and not only on their sum. Hence, there is only one “equilibrium” solution for the densities minimizing the functional. However, since there is no distinction between the two species in the excess free energy, Eq. (8), the resulting density profiles represent the most disordered state that is compatible with the total interaction. In the canonical case with $N=2$ particles, this means that the profiles Eq. (27), corresponding to the dashed lines in Fig. 1, result in a smaller value of the canonical density functional than in the physical equilibrium state. As a logical consequence, choosing Eq. (26) as the initial profiles, the PCD for mixtures spuriously drive the system out of the actual equilibrium, which we show in Fig. 6(a). This behavior illustrates clearly and already at early times the unphysical artifact observed in Sec. III.2 that the particles tend to mix. Intriguingly, the PCD for mixtures also spoils the time evolution of the total density profile, which becomes obvious from Fig. 6(b), where the system, seemingly initiated in equilibrium, exhibits a non-trivial dynamical behavior, just to finally return to a state with the same total density, but different single-particle profiles. This is a clear indication that the present form of the PCD do not reproduce the (correct) results of the one-component version canonical4 .

The situation described above is somehow reminiscent of the mixing paradox, which tells us that one should not assign a different entropy to a mixed and a demixed system (separated by a wall) of ideal particles if one is not able to measure or does not care about the physical difference between two species. In this case, no entropy change upon mixing or reseparation may occur. For the mixture of one-dimensional hard rods considered here, the particle interactions take the role of a wall inserted into the system. In inversion of the argumentation for the mixing paradox, we expect a higher entropy (or lower free energy) for the demixed state (true equilibrium), where the mixed state should even be entropically forbidden. This means that, in our theory, we care about a difference that is not reflected by the mathematical structure of the PCD for mixtures. Therefore, the entropy (or the canonical free energy) employed for a mixture is ill-defined.

To resolve this “inverse mixing paradox”, we continue the discussion from Sec. II.4. Both partition functions $Z_{2}$ from Eq. (19) and $Z_{11}$ from Eq. (21) imply a well-defined entropy. In the DFT language, the first corresponds to the intrinsic free energy functionals $\mathcal{F}_{\text{id}}[\rho]$ (all particles distinguishable) and $\mathcal{F}_{\text{ex}}[\rho]$ [derived for the symmetric pair potential $u(|x|)$] of the total density of all particles (species). The form of $Z_{11}$ [or generally $Z_{N_{1}N_{2}}$] should be represented by $\mathcal{F}_{\text{id}}[\{\rho^{(\nu)}\}]$ (two indistinguishable species) and an excess term $\mathcal{F}_{\text{ex}}[\{\rho^{(\nu)}\}]$, explicitly depending on the species-labeled profiles. Such a (yet unspecified) functional should be based on the asymmetric interaction potential $o(x)\neq o(-x)$ from Eq. (24), which preserves the order of the particles (species), i.e., the non-ergodicity of the system, and thus allows for a physical distinction. In contrast, the functional introduced in Sec. II.2 as the starting point of PCD, corresponds to an increased partition function $2Z_{2}$ [or generally $Z_{N}N!/(N_{1}!N_{2}!)$], which means that the ideal free energy $\mathcal{F}_{\text{id}}[\{\rho^{(\nu)}\}]$ implies the combinatorics of two species, whereas the excess free energy $\mathcal{F}_{\text{ex}}[\rho]$ is the same as for (orderwise) indistinguishable particles. More generally, this means that $\mathcal{F}_{\text{ex}}$ is built from symmetric pair potentials $u(|x|)$ of possibly physically distinguishable particles. Statically, this overcounting of states does not change the canonical equilibrium density profiles $\rho^{(\nu)}_{2}(x)$, ensuring the correct result for the total density $\rho_{2}(x)$ according to Eq. (27). In the dynamical case, however, the inconsistency between the entropic force [first term in brackets in Eq. (29), related to $\mathcal{F}_{\text{id}}$] and the interaction force [second term in brackets in Eq. (29), related to $\mathcal{F}_{\text{ex}}$] is the ultimate reason for the wrong time evolution of $\rho_{2}(x,t)$.

To illustrate the consequences of applying a theory built from symmetric interactions to a non-ergodic system, let us consider the time evolution of the adiabatic potential in Fig. 7(a), corresponding to the density profiles from Fig. 6(a) of two particles. Before the first time step, the potential is infinitely steep at the points where the density becomes zero. This shows that the initial confinement is not intrinsically described by the interaction functional but has to be artificially generated. Thus, there is a net force due to $V^{(\nu)}_{\text{ad}}$ in Eq. (32) that drives the dynamics of each particle into the physically forbidden region. At later times, the generating external fields become less and less restrictive on the interpenetration of the particles and become equal (up to a constant) to the external potentials $V_{\text{ext}}^{(\nu)}$ when the equilibrium state of the underlying functional is reached. It must be the goal to describe the intrinsic interactions in a way that $V^{(\nu)}_{\text{ad}}(x)\equiv V_{\text{ext}}^{(\nu)}(x)$ in the true equilibrium state.

Finally, we show in Fig. 7(b) that, at early times and in the low-density regions, the PCD of one particle in a two-particle mixture is remarkably similar to the proper adiabatic dynamics of a single particle ($N=1$) with the initial condition being identical to one of the species-labeled density profiles for $N=2$. Mathematically, this can be easily explained by the similarity of the adiabatic potentials in the regions where the density of the corresponding species is small, compare Fig. 7(a), so that also in the two-particle system the main contribution to the free energy stems from the interaction with the generating adiabatic potential. This behavior suggests that the speed with which the two profiles mix, which would represent a “hopping rate” on the particle level, is independent of the (local) density.

IV Discussion

In conclusion, we have shown that the presented generalization of PCD to mixtures provides a very good description of the early relaxation dynamics of individual Brownian particles in an interacting system, in particular, it clearly improves on the grand-canonical DDFT results. Only at later times, the PCD exhibits an unphysical mixing behavior in one dimension, which dominates any deviations arising from the neglect of superadiabatic forces. Our study of the somewhat artificial one-dimensional case both provides valuable insights into more realistic systems in higher spatial dimensions and is of fundamental theoretical interest in its own right.

We stress that our approach is also relevant for systems much larger than those considered here and even in bulk. By choosing appropriate external potentials it is possible to isolate a single particle in the initial state. Therefore, the difference between the canonical and grand-canonical ensemble remains significant for the (single-file) dynamics even though the differences decrease rapidly for the static properties and joint density of all particles.

To properly describe the non-ergodic Tonks gas, statistical mechanics only constitute a workable approach if an asymmetric interaction potential is considered which does not only depend on the relative distance between two particles. For the present variational approach to work out, one would have to construct a DFT based on such an interaction potential between members of different species. Such a mixture would then additionally be non-additive, a case in which even for a symmetric potential in one dimension only an approximate hard-rod functional can be derived schmidtNAM .

Having an accurate ergodicity-breaking theory in one dimension would be very instructive, since one could gain more explicit insight into superadiabatic contributions on the single-particle level. A more general and formally exact variational approach, based on a mixture of DFT methods and statistics intrinsically respecting the particle order, will be presented in a future publication. Another promising route would be to analyze the exact results for hard-rod dynamics provided by Lizana and Ambjörnsson LA1 ; LA2 to see if their lengthy analytical expressions obtained by a Bethe ansatz can be reformulated within the context of a variational approach. It will also be insightful to provide a simulation setup that can reproduce the mixed canonical density profiles in equilibrium, which will be a challenging task for both BD and (dynamic) Monte Carlo. This would allow to compare the adiabatic-superadiabatic splitting that applies to the artificial case prescribed by PCD with the (known FortiniPRL ) splitting that applies to the true BD case with preserved particle order.

In higher dimensions, the problem described above is seemingly resolved. Most importantly, these systems are, in general, ergodic, so that the final equilibrium state of PCD is the same as in BD and the pair potential must also be symmetric. We showed, however, that this does not guarantee a correct description of the transient dynamical regime [cf. Fig. 6(b)]. By contrast it is very likely that there are still some artifacts in the PCD due to its insensitivity to dynamical caging effects. For example, one would expect in BD that the mixing of the particle-resolved density profiles becomes increasingly slow with increasing total density, which is not reflected by the conclusions drawn in one dimension from Fig. 7(b). More specifically we hold these artifacts accountant for the fact, that density-based theories in general overestimate long-time diffusion constants stopper_long , as it adds unphysical particle exchange dynamics to the physical circuiting of particles, which is slower in dense suspensions of any dimension. In this sense, we have performed a minimalistic but extreme test (with infinite circuiting time) for the caging scenario in higher dimensions. It will be challenging to study in detail how significantly this shortcoming of PCD would influence the adiabatic dynamics, since there exists no exact grand-canonical functional for interacting systems in higher dimensions. Moreover, to describe caging statistically, one would require a complex many-body interaction (and not only an asymmetric pair potential), which leaves not much hope for a theoretical implementation. Despite these caveats, the available approximate forms of DFT in three dimension are very accurate FMT ; Tar99 ; HaRo06 , so that we expect that the corresponding PCD will provide a pretty good account for the early (adiabatic) dynamics of single-particle profiles, especially at a low overall density.

Exceptions to the above are presented by glassy, jammed or otherwise arrested systems. In such cases the interparticle coupling is so strong (due to, e.g., high density or strong attractions) that the phase space can not be fully explored and canonical averaging is not appropriate. In particular, for high-density systems with purely repulsive interactions the cage effect is a dominant physical mechanism. The failure of the ‘three-rod-caging test’ in Fig. 3 ultimately points to the impossibility of describing such glassy states using theories for which the density is the only variable. The standard observable used to quantify the dynamics of arrested states is the van Hove function, the self part of which describes the dynamics of a single tagged particle. DDFT has been employed in the test particle limit to approximate this self van Hove function dtp ; stopper_bulk . In the light of our present study one can safely conclude that the dynamic arrest which has been observed in DDFT calculations is an artifact arising from an approximate free energy functional dtp , rather than a true indication of vitrification.

To provide a more sophisticated description of caging effects on the level of symmetric pair potentials in any dimension, it seems unavoidable that one must extend the variational approach beyond the one-point density alone. For example, variational approaches based on two-point correlations might be better able to cope with caging. Alternatively, to obtain improved results using the one-point density alone in a non-variational framework, it may be possible to incorporate superadiabatic effects by relaxing the requirement of time locality, i.e., incorporate memory functions. By far the most natural extension of DDFT (in the authors opinion) lies in the framework of power functional theory pft1 ; pft2 , based on a functional of both density and current, which is nonlocal in both space and time. Indeed, it seems clear that a vector field is necessary in order to describe the motion of a fluid. So far, workable approximations to the power functional have remained time-local pft2 . It remains, however, unclear, whether implementing superadiabatic effects will automatically provide a better description of caging, which, as laid out in this paper, can be understood as an adiabatic many-body effect, for the description of which we require an accurate treatment of the static interactions.

Acknowledgements.

T.S. was supported by the Deutsche Forschungsgemeinschaft as part of the Forschergruppe GPSRS under grant ME1361/13-2.

Appendix A Brownian dynamics simulations

For calculation of trajectories in our BD simulations we employ a hybrid algorithm. It treats random forces and soft external potentials in discrete time steps according to the standard BD approach. The interactions between particles or between particles and walls, however, are treated as instantaneous collisions to account for the hard core repulsion. We chose this approach over the usual approximation of hard core interactions with steep soft potentials to avoid smoothing effects and instead obtain exact hard core density profiles.

In each time step of length $\Delta t=5\times 10^{-5}\tau_{\mathrm{B}}$ we calculate the instantaneous velocity $v_{i}(t)$ of a particle $i=1,...,N$ (of species $\nu$), which, in overdamped Langevin dynamics, directly adjusts to the force $D_{0}\beta{\mathbf{f}}_{i}(t)=v_{i}(t)$ acting on the particle. The force reads

[TABLE]

where $R_{i}(t)$ is a random force drawn from Gaussian white noise with variance $2/(D\beta^{2}\Delta t)$ and zero mean, the particular harmonic trap potentials $V_{\text{ext},0}^{(\nu)}$ are specified in the caption of Figs. 2 and 3 and the Heaviside step function represents the switch-off of the harmonic trap at $t=0$. From the old position $x_{i}(t)$ and the velocity we can compute the preliminary new position via

[TABLE]

If particles overlap with each other or the walls, then we resolve these overlaps in billiard-like collisions. For each created overlap we first calculate the collision time,

[TABLE]

for a collision of particles $i$ and $i+1$, or

[TABLE]

for a particle-wall collision, where $x_{\mathrm{w}}=R,L-R$ is the minimum or maximum position of a particle in the slit. Then, these collisions are scheduled and resolved chronologically to obtain the final positions of the time step via

[TABLE]

for particle particle collisions or

[TABLE]

for particle-wall collisions. (If hereby new overlaps are created, then the collisions, again, have to be scheduled and resolved.) When all collisions are resolved, the algorithm resumes with calculating the velocities for the next time step.

To sample the time dependent density profiles we generate an equilibrated configuration with harmonic traps switched on and an initial equilibration time of $1\tau_{\mathrm{B}}$. Between each sample we equilibrate the initial configuration for another $10^{-3}\tau_{\mathrm{B}}$ to decorrelate the samples from each other. A sample is then taken by switching off the harmonic traps and relaxing the system to its final state. The density profiles of the species shown in Sec. III.2 are obtained by averaging over $4\times 10^{6}$ individual samples via

[TABLE]

where the sum runs over all particles belonging to the respective species.