Collective Motion of Repulsive Brownian Particles in Single-File Diffusion with and without Overtaking

Subdiffusion is commonly observed in liquids with high density or in restricted geometries, as the particles are constantly pushed back by their neighbors. Since this “cage effect” emerges from many-body dynamics involving spatiotemporally correlated motions, the slow diffusion should be understood not simply as a one-body problem but as a part of collective dynamics, described in terms of space–time correlations. Such collective dynamics are illustrated here by calculations of the two-particle displacement correlation in a system of repulsive Brownian particles confined in a (quasi-)one-dimensional channel, whose subdiffusive behavior is known as the single-file diffusion (SFD). The analytical calculation is formulated in terms of the Lagrangian correlation of density fluctuations. In addition, numerical solutions to the Langevin equation with large but finite interaction potential are studied to clarify the effect of overtaking. In the limiting case of the ideal SFD without overtaking, correlated motion with a diffusively growing length scale is observed. By allowing the particles to overtake each other, the short-range correlation is destroyed, but the long-range weak correlation remains almost intact. These results describe nested space–time structure of cages, whereby smaller cages are enclosed in larger cages with longer lifetimes.


Introduction
Particles in dense liquids are hindered from free motion, being constantly pushed back by their neighbors. This is often described as a "cage" that confines each particle. The cage effect makes the motion subdiffusive and, in certain cases, leads to the glass transition [1,2].
To be specific, let us consider a system consisting of Brownian particles with a nearly hardcore interaction. The position vector of the i-th particle, r i = r i (t), is governed by the Langevin equation with m and µ denoting the mass and the drag coefficient of the particle, µf i (t) representing the thermal fluctuating force, and the interaction being prescribed as U = U(r 1 , r 2 , . . .) = ∑ (j,k) in terms of the pair potential V jk . Among the fundamental statistical quantities characterizing this system is the mean square displacement (MSD), i.e., the second moment of the displacement R i = r i (t) − r i (0). If the interaction through U is negligible, each particle diffuses freely so that R 2 i grows in proportion to t (for timescales longer than m/µ). This occurs when the colloidal fluid modeled by Equation (1) is dilute enough. In contrast, particles in a denser colloidal fluid are hindered from free motion by the cage effect, so that the growth of R 2 i is much slower. As is illustrated in Figure 1a, every particle in such a system is almost arrested in a "cage" consisting of its neighbors. In extreme cases, the system ceases to be fluid and becomes a kind of amorphous solid, referred to as colloidal glass [1].
Although it is true that the cage effect suppresses the growth of MSD on the whole, the details are rather complicated [2]. The behavior of the MSD in dense liquids reflects at least three aspects of caged dynamics: nearly free motion within the cage for a short time, possible drift of the cage enclosing the particle at a longer timescale, and hopping of the particle out of the cage as a rare event. Proper characterization of these processes requires space-time description, typically in terms of some four-point space-time correlation [2][3][4], as the cage effect actually emerges from many-body dynamics involving collective motions of numerous particles correlated both spatially and temporally. In search of insight into the theoretical treatment of such collective motions, here, we take note of the one-dimensional (1D) system illustrated in Figure 1b, following several authors who studied it as a simplified model of the cage effect [4][5][6][7][8][9][10]. The slow dynamics of such a 1D system are known by the name of single-file diffusion (SFD). In what we call the ideal SFD, every particle is eternally trapped within the "cage" formed by their neighbors. The MSD in the ideal SFD is known to grow subdiffusively as R 2 i ∝ √ t [11][12][13][14] (for 1D systems, we write R i instead of R i ). The subdiffusion in SFD emerges from collective motion of particles [5,10,15] and is also detected as a negative longtime tail in the velocity autocorrelation [16,17], indicating that the particle is pushed back by its neighbors. The importance of the collective motion is understood by considering the origin of the effective stochastic equation for a single particle in SFD: the one-body equation (yielding the negative velocity autocorrelation) is actually based on the collective dynamics described in terms of the fluctuating density field [18].
Focusing on the collective motions in SFD, the group of the present authors has noticed the usefulness of the displacement correlation R i R j [4,10,19]. It is a kind of four-point space-time correlation that probes both the time scale and length scale; the definition of the displacement includes t, while the spatial scale is included as the mean distance between the two particles (i, j). In the ideal SFD in which the particles are forbidden to overtake each other, the displacement correlation has been calculated both analytically and numerically [4,10,20]. The calculated displacement correlation revealed collective motions behind the slow diffusion in SFD, in contrast to free diffusion in which R i R j vanishes (unless i = j). The formalism for analytical calculation of the displacement correlation can be extended to the case of two-dimensional (2D) colloidal liquids [10,21], which reproduces some numerical findings in 2D systems, such as vortical cooperative motion, with negative velocity autocorrelation being a manifestation of the cage effect. One of the delicate points in this extension is that the cage effect in 2D liquids cannot be infinitely strong, in the sense that eventually the particles can escape from the 2D cage. The escape from the cage is an important process, which requires further investigation.
Methodologically, it should be noted that one of the most powerful approaches to SFD, employing the relation between the position of the tagged particle and the density fluctuation [10,12,17,22], has been formulated in reliance on the assumption that overtaking is completely forbidden. To extend this formalism to the case of "non-ideal SFD"-in which overtaking is allowed-is a challenging problem, which is the main objective of the present work. Most of the existing works on SFD with overtaking have reported numerical simulations [23][24][25][26][27], while analytical results are quite rare. In the exceptional case of lattice SFD, the method of vacancy dynamics [15,16] has been applied to quasi-1D geometries allowing some kind of overtaking [28,29]. The analysis of lattice SFD with overtaking, however, is not readily extensible to the cases with a continuous space coordinate.
In the present work, we discuss how the analytical results on the displacement correlation in SFD [4,10,20] are modified, if the particles are allowed to overtake each other and thereby escape from the quasi-1D cage as a rare event. The probability of the escape is regulated by the height of the potential barrier, denoted by V max , so that V max → +∞ and V max → 0 correspond to the ideal SFD and the free diffusion, respectively. Some numerical solutions for finite V max were included in our previous work [19], but analytical calculations were limited to the ideal case without overtaking, for the very reason that overtaking was difficult to take into account in the theoretical framework based on the density fluctuation. Now, the effect of a non-zero overtaking rate on the displacement correlation will be shown analytically as a main result of the present work.
Logical presentation of the main results in Section 4 requires a considerable amount of review in preparatory sections. For this reason, the paper is organized as follows: In Section 2, after the governing equation of the 1D system is specified and the collective motion is illustrated in a space-time diagram, we define some basic concepts and variables, such as the displacement correlation, overtaking, the fluctuating density field ρ(x, t), and the label variable ξ, clarifying their background. In particular, the kinematics of overtaking are discussed in Subsection 2.5. The usage of the label variable is a keystone for the analytical calculation of the displacement correlation [4,10,30], as is reviewed in Section 3 in the case of SFD without overtaking. The displacement correlation in this case is expressed in terms of a similarity variable, implying a nested space-time structure of cages. Subsequently, we proceed to the main topic in Section 4, in which we incorporate the effect of overtaking into the calculation of the displacement correlation by considering the dynamics of overtaking in terms of the Ξ i (t) prepared in Subsection 2.5. It is shown analytically and confirmed numerically that the infrequent overtaking events destroy the short-range correlation, while the long-range weak correlation remains almost intact. The final section is allotted for discussion and concluding remarks.

Specification of the System
We consider a 1D system of Brownian particles with short-range repulsive interaction, confined in a narrow channel, as is depicted in Figure 1b. With the position of the i-th particle denoted by X i = X i (t), the system is governed by the 1D Langevin equation: where m and µ represent the mass and the drag coefficient of the particle, respectively. The system contains N particles within the length L. Posing the periodic boundary condition, X i+N = X i + L, we consider the limit of N → ∞ with the density ρ 0 def = N/L kept constant. The thermal fluctuating force, µ f i (t), is characterized by the variance, with T denoting the temperature of the medium. The whole system is assumed to be at thermal equilibrium, which implies spacial homogeneity and temporal steadiness. The interaction between the particles is expressed by the pairwise potential, V(r). We could choose any family of V(r) that interpolates between the limiting case of V = 0 and the opposite limit of the hardcore potential, with the diameter σ. Here, we choose which is parametrized by the barrier height V max . We also tested some other potentials [19], only to find that the basic behavior of the 1D system is qualitatively unaffected by different choices of V(r).
Preference was given to Equation (5) merely because its hard sphere limit (V max k B T) has been studied systematically [31] in the context of 3D glassy dynamics.
In regard to the system governed by Equation (2), we refer to the case of V max → +∞ as the ideal SFD, in which every particle is eternally caged by its neighbors. Large but finite values of V max allow the particle to exchange positions with one of its neighbors as a rare event which we call overtaking (borrowing the word from traffic flow). The ideal SFD means SFD without overtaking, and we may say "non-ideal SFD" referring to the case of finite V max /k B T. Note that the description of non-ideal SFD with the 1D equation (2) can be interpreted as modeling a quasi-1D system [19,[23][24][25][26][27] in which, typically, the position vector r i = (X i , Y i ) is governed by Equation (1a) with the potential term where V ex = V ex (y) denotes the external confinement potential such that V ex (±∞) → +∞. In this description, V(X k − X j ) in Equation (2) represents the free energy of the subsystem consisting of the neighboring particles j and k.
The specification of the system by Equations (2), (3) and (5), supplemented with the periodic boundary condition, involves some dimensional constants. As the basic scales of the length and the time, we take the particle diameter (σ) and the corresponding diffusive time (σ 2 /D), where D = k B T/µ is the diffusion constant of a free Brownian particle. A finite value of mass, such that m/µ : σ 2 /D = 1 : 1, is specified for computational ease, unless specified otherwise. The system size, L, must be infinitely large; though, in numerical computations, we must specify some finite values for it. For later convenience, we introduce 0 def = L/N = 1/ρ 0 which has the dimension of length. The nondimensional barrier height, V max /k B T, has an effect on the dynamics through the overtaking frequency, as will be discussed later.
In numerical simulations, the system is equilibrated by a preparatory run started at t = −T w . Subsequently, for a reason clarified in the next subsection, the particles are renumbered consecutively in the sense that at t = 0. It should be noted that T w must be longer than max t for sufficient equilibration [4].

Spatiotemporally Correlated Motion in SFD
As a graphic depiction of collective motions in SFD, let us examine Figure 2, in which a numerical solution of Equation (2) in the case of the ideal SFD is represented as worldlines in the (x, t)-plane.
To visualize the correlation of the worldlines, we measured the displacement for each particle i, for the time interval from 0 to t = 2 n × 10 σ 2 /D (with n = 1, 2, . . .), in accordance with Ref. [4]. If R i (t) > 5σ, the position of the particle is marked with a filled circle (•); if R i (t) < −5σ, it is marked with an open square ( ). As the time difference (t) increases, a string of the same kind of symbol is formed, expressing a cluster of particles moving together in the same direction.  , is plotted as worldlines in the (x, t)-plane (note that the t-axis is on a logarithmic scale). The symbols • and mark particles displaced (by more than 5σ) rightward and leftward, respectively.
The formation of clusters, visually shown in Figure 2, is quantified by calculating the displacement correlation R i R j . Note that the consecutive numbering in Equation (7) is needed to make R i R j meaningful as a function of j − i ( = ∆) and t. Since R i and R j have the same sign within the same cluster, their product must be positive if the distance in the numbering, ∆ = j − i, is small, while R i R j for distant particles (with 0 ∆ greater than the correlation length) is expected to vanish. The average, denoted by , is taken over the initial condition and the Langevin noise. Computationally, the displacement correlation is calculated as

Continuum Description
Our theoretical approach to the displacement correlation, R i R j , is based on a continuum description of the dynamics of Brownian particles. As fluctuating hydrodynamic fields describing the temporally coarse-grained dynamics of {X i } i=1,2,... for timescales longer than m/µ, one may take the density fields and their fluxes Note that the delta function in Equation (10) should be regarded as a blunted one, as a result of the coarse-graining (see §II-B in Ref. [21] and references therein). The density field, ρ(x, t), is governed by the Dean-Kawasaki equation [32][33][34][35][36][37], which can be presented as a set of equations of the following form: The effective potential, V eff in Equation (12c), is determined by the condition that the density fluctuation, described by Equations (12), should be consistent with the static structure factor, determined directly from Equation (2). More specifically, with the Fourier mode of the density and its correlation defined aŝ evidently, F(k, t = 0) must be equal to S(k), and the initial decay of F(k, t) is known to be exponential, as is shown in Ref. [38] as Equation (4.144). In our notation, it reads with D c = D c (k) = D/S(k) referred to as the (short-time) collective diffusion coefficient [38,39]. Once V eff is determined so as to make the linearized dynamics of Equations (12) consistent with Equation (14), we may redefine S(k) by the ratio of D to D c (k).

Label Variable
Now let us outline the main idea that makes it possible to calculate R i R j analytically on the ground of Equations (12) [4,10,19,30]. The point is to introduce a new variable, ξ, referred to as the label variable, to incorporate the notion of particle tracking into the mathematical formalism that links R i R j with Equations (12).
The necessity of the new variable for particle tracking is understood by noticing how a more straightforward approach, based on ρ i (x, t) in Equation (10), is confronted by a difficulty. In principle, R i R j can be obtained from the Fourier mode of ρ i (x, t), because the definition If the correlation on the left side is successfully evaluated by nonlinear analysis of the equation governing ρ i and ρ j , analogous to Equation (12b) and later shown as Equation (48), then R i R j can be obtained from the power series on the right side. This is insurmountable, unfortunately, as is evident from the complication encountered in the apparently easier problem of evaluating ρ i (k, t)ρ i (−k, 0) in SFD [8,9]. The difficulty originates from the choice of the standard field representation with the independent variables (x, t), referred to as the Eulerian description (according to the terminology of fluid mechanics [40,41]). Since R i R j is treated as a kind of four-point space-time correlation [42,43], it implies a four-body correlation in the Eulerian description, as is seen in Equation (16). However, this difficulty can be avoided by switching to another way referred to as the Lagrangian description [30,40,41], which means to include particle tracking mechanism in the definition of the fields and their correlations. In this way, the displacement correlation can be treated simply as a two-body Lagrangian correlation of the field [4,10,19].
Deferring a concrete calculation of R i R j until the next section, here we only define the label variable ξ to lay the foundation for it. Instead of (x, t) for the standard space-time coordinate system, we introduce a stretchable coordinate system (ξ, t), requiring ξ = ξ(x, t) to satisfy the convective equation, which states that ξ should be convected with the velocity u = Q/ρ. To satisfy Equation (17), we define ξ as a solution to This is solvable because of the continuity equation (12a), with the solution determined uniquely by some initial condition, such as ξ(X 0 (t 0 ), t 0 ) = 0. Subsequently, by inverting the mapping (x, t) → ξ, we obtain the coordinate system with the independent variables (ξ, t) [4,10,19,21,30]. It should be emphasized that we take Equation (18), not Equation (17), as the definition of the mapping between ξ and x. In other words, we do not solve Equation (17) in the usual sense of the word; rather, Equation (17) is satisfied as a consequence of Equation (18). In this way, we avoid complication of the attempt to define ξ = ξ(x, t) directly with Equation (17), which would require specification of the initial condition We also remind the readers that the delta function in the definition of ρ is a blunted one, as has been noted immediately after Equations (10) and (11).
Using the label variable ξ, defined in this way, we can calculate R i R j analytically. The (ξ, t) coordinate system has an advantage of making it easy to trace the worldlines, such as the ones plotted in Figure 2, because ξ is expected to keep the same value if we follow the identical particle. In order to see it, we define as a function of the particle number i and the time t, and consider its t-derivative [30,44]. Provided that the conventional chain rule is valid, the time-derivative of Ξ i (t) is where the definition of ξ in Equation (18) is taken into account. The expression on the right side of Equation (20) vanishes unless the i-th particle overlaps with another. In the ideal SFD, in which the particles never overlap, Ξ i (t) is none other than the numbering i; this is the key ingredient for the analytical calculation of R i R j in the ideal SFD.
In the absence of overtaking, the particles can move only as a result of changes in inter-particle spaces, which is illustrated schematically in Figure 3a as a migrating "vacancy" causing correlated motion of particles. The overtaking allows another kind of motion, illustrated in Figure 3b, which does not require migration of a vacancy. Before proceeding to a concrete calculation of R i R j , let us discuss how to describe an overtaking event within the framework of the (ξ, t) coordinate system.

Kinematics of Overtaking Events
In the presence of overtaking, Ξ i (t) = ξ(X i (t), t) is not necessarily equal to i. In this case, we must distinguish between them and discuss correspondence among three variables: the numbering (i ∈ Z), the label variable (ξ ∈ R), and the position (x).
We still define ξ by Equation (18) without regard to overtaking. The constant of integration is chosen appropriately so that Ξ i (t) = ξ(X i (t), t) has an integer value unless the i-th particle overlaps with another. On this premise, the mapping from i ∈ Z to Ξ i (t) = ξ(X i (t), t) ∈ Z is injective, which means, so to speak, a kind of exclusion principle whereby different particles must carry different values of ξ.
The overtaking process is a transition from one such mapping to another. Since the effect of overtaking on the mapping i → Ξ i (t) is local, as is readily seen from Equation (18), we may identify the overtaking event with a process in which two particles exchange the value of their label. In the case of a pair of particles (i, j), with t 1 and t 2 denoting points in time immediately before and after the exchange, this process is described as To justify Equation (21), we consider what occurs in the time interval from t 1 to t 2 . Since the particles may overlap in the course of overtaking, the value of Ξ i (t) is not limited to Z but rather belongs to R. Using ρ and ρ i given in Equation (10), we rewrite Equation (20) as which implies that Ξ i (t) is locally conserved. Note that only overlapping particles contribute to the integral. If no particle overlaps the i-th one, then the expression on the right side of Equation (20) vanishes so that Ξ i remains constant. If only one particle, say the j-th one, overlaps particle i, Equation (22) gives By integrating Equation (23) over the time interval from t 1 to t 2 and taking the conditions such as so that only two cases are possible: one is the case of unsuccessful or retracted overtaking, with the same value of (Ξ i , Ξ j ) restored after the interaction, and the other case corresponds to Equation (21) that describes (successful) overtaking. Note that the other possibilities are eliminated by the "exclusion principle" for the configurations before and after the overtaking process. As a rare case, three-body or four-body interactions might be possible, but it suffices to approximate such a case with a sequence of two-body exchange processes. Changes of Ξ i (t) = ξ(X i (t), t) and X i (t) = x(Ξ i (t), t) are exemplified in Figure 3. Since ξ = ξ(x, t) satisfies ∂ξ/∂x = ρ by definition in Equation (18), {Ξ i } is always spatially consecutive, as is shown in Figure 3 with small digits. The important point is that this spatial consecutiveness of {Ξ i } holds true even if the numbering is inverted by overtaking. In other words, ξ labels the order in the file and not the particles themselves. In the case of Figure 3b, the particles are initially numbered consecutively so that Ξ i = i for all i, until particles 4 and 5 start to overlap. During overtaking, the values of Ξ 4 and Ξ 5 evolve according to Equation (23) with (i, j) = (4, 5), while the labels for other particles, such as Ξ 2 , Ξ 3 and Ξ 6 , remain unchanged. Finally, after overtaking, particle 4 carries the label Ξ 4 = 5, while label 4 is now carried by particle 5 so that Ξ 5 = 4.
In the next section, we start with the case of the ideal SFD, in which the requirement of (d/dt)Ξ i (t) = 0 is fulfilled so that ξ simply plays the role of particle numbering. Later, we consider the temporal change of Ξ i (t) to allow for overtaking in Section 4.

Analytical Calculation of Displacement Correlation
Here we review the analytical calculation of R i R j in the ideal SFD [4,10,19], in which Ξ i is independent of t.
Out of the three equations composing the Dean-Kawasaki equation (12), the continuity equation (12a) is already included in ξ = ξ(x, t), or its inverse mapping x = x(ξ, t), through Equation (18). To relate x = x(ξ, t) to the remaining two equations, we solve Equation (12b) for u = Q/ρ, which should be equal to u = ∂ t x(ξ, t). Noticing that u(ξ, t) is the (negative) flux of 1/ρ = ∂ ξ x(ξ, t), we introduce [30] so that 1/ρ = 0 (1 + ψ). The fluctuating field ψ can be interpreted as a continuum representation of migrating vacancies [16] or elongation of a chain [45]. As is illustrated schematically in Figure 3a, migration of a vacancy can give rise to correlated motion of particles.
Once the dynamics of ψ are known from Equations (12), the remaining task for the calculation of R i R j is only a problem of kinematics-given the field ψ, how can we find the displacement? This is readily solved in a Fourier representation, allowing us to express x = x(ξ, t) as an antiderivative of ψ: where X G (t) corresponds to the center-of-mass motion which should be negligible in the limit of N → +∞ [10]. The positions of the i-th or j-th particle are then obtained by substituting Ξ i or Ξ j into ξ in Equation (26), which readily yields (27) in the absence of overtaking (i.e., dΞ i /dt = 0). To calculate R i R j , we multiply Equation (27) by its duplicate with (i, k) replaced by (j, −k). Subsequently, taking it into account that ψ (k, t)ψ(k , t ) generally vanishes unless k + k = 0 (due to the space-translation symmetry of the system), we find where we have defined and used Ξ j − Ξ i = j − i = ∆. We refer to Equation (28) as the Alexander-Pincus formula [10,12], which relates the displacement correlation to C ψ . Since C ψ is a two-body correlation, it is much more tractable than the four-body correlation in Equation (16).
To allow concrete calculation of C ψ , Equations (12b) and (12c) are rewritten as an equation foř ψ(k, t) in the following form: The statistics of the random force term are specified as where D * = D/ 2 0 . The coefficient of the linear term is which also gives S = D/D c = D * /D c * , and the coefficients V pq k in the nonlinear term are found in Refs. [4,30].
Within the linear approximation to Equation (30), C ψ is readily calculated as Although a nonlinear theory producing a correction term to be added to Equation (33) is also possible [4], here we ignore this correction, giving priority to the simpler expression in Equation (33). The displacement correlation is obtained by substituting Equation (33) into the Alexander-Pincus Formula (28) and evaluating the integral. As the contribution from the longwave modes is dominant, S = S(k) and D c * = D * /S(k) can be replaced by their limiting values for k → +0. Thus, R i R j is obtained explicitly as a function of ∆ and t, expressible in terms of a similarity variable as [4,10,19] The dynamical correlation length λ(t) grows diffusively, which seems to be consistent with the observation of growing clusters in Figure 2.
Note that infrared and ultraviolet cutoffs are not necessary in Equation (28), as the integrand is regular around k = 0 and decays for k → ±∞ (algebraically, but fast enough). This should not be confused with the infrared divergence of x(ξ, t) in Equation (26) in the limit of L → ∞.

Particle Simulation in the Absence of Overtaking
The theoretical prediction in Equation (35) is tested in Figure 4a, by plotting the values of the displacement correlation, R i R j , obtained from numerical simulation of the ideal SFD. The computed values are plotted in terms of rescaled variables; according to Equation (35), a plot of R i R j /(K √ t) against the similarity variable θ should give a single master curve for all values of the time interval t. The ideal SFD was simulated by solving Equation (2) numerically with a very high barrier (we chose V max = 50k B T). The system size and the density were specified as N = 10 4 and ρ 0 = N/L = 0.2 σ −1 so that L = 5Nσ.
The three kinds of symbols in Figure 4a correspond to three different values of t. The plots for all these values of t are seen to be reducible to a single master curve given by ϕ( · ) in Equation (35), supporting the prediction of the analytical calculation.
To be precise, a small but finite discrepancy is found at θ = 0 for t = 10 σ 2 /D. One may be tempted to explain this short-time discrepancy simply as indicating a lack of time to establish correlated motion, because a particle, on average, takes time on the order of 1/D c * to encounter its neighbors. Unfortunately, this argument appears too simple to explain the numerical result in Figure 4a in which a positive correlation for θ = 0 (i.e., |∆| ≥ 1) has already been established at t = 10 σ 2 /D. We note, on the other hand, that the discrepancy can be ascribed to the nonlinear term in Equation (30) which was ignored in the previous subsection. It is shown that inclusion of the nonlinear term gives a correction to Equation (35), which is significant only for |θ| 1 and D c * t < 1 [4], and the sign of the correction for the MSD is negative [4,10]. Thus, the short-distance correlations are affected by nonlinear coupling ofψ, while correlations over a long distance seem to be tractable with a linear theory.
It should be also noted that the theoretical predictions discussed here are based on the Dean-Kawasaki equation in which inertia is completely ignored, while the particle-based simulation is performed with finite m/µ. To check for consistency between the numerical simulation and the theoretical predictions, three cases with different values of m/µ are compared in Figure 4b. It is seen that the change in m/µ does not make any remarkable difference. This is in agreement with the general expectation that the Langevin dynamics on time scales longer than m/µ are basically independent of the inertia, because the momentum can be eliminated by temporal coarse-graining [46][47][48].
In regard to correlations over long distances, one might be tempted to suppose that the static structure factor, S(k) in Equation (13), could be helpful in the detection of such long-ranged correlations. This point was discussed in Ref. [19], leading to the conclusion that the structure of the collective motion is not properly captured by the static structure factor. It is for this reason that we focused on R i R j rather than S(k).

Effects of Overtaking on Displacement Correlation
Having reviewed the analytical calculation of the displacement correlation R i R j , which leads to Equation (35) for the ideal SFD, let us now consider effects of overtaking that were ignored in the previous section. We start with numerical observation, noticing how the behavior of R i R j deviates from Equation (35) due to overtaking. This deviation is then compared with a modified theory in which overtaking is allowed for.

Particle Simulation of SFD with Overtaking
The governing equations of the system, namely Equations (2), (3) and (5), contain a parameter V max representing the barrier height. This parameter, V max , regulates the frequency of overtaking, if the other parameters are kept unchanged.
Two extreme cases are already known theoretically: the case of V max → ∞ implying the ideal SFD in which overtaking is completely forbidden, and V max = 0 corresponding to free diffusion in which overtaking is always allowed. In the ideal SFD, there is a positive correlation between displacements of two particles, as is shown in Equation (35), while in free diffusion, displacements of two particles are totally uncorrelated, as is easily seen by proving for V max = 0. Between these two extreme cases, there are cases of finite values of V max , allowing overtaking with some probability. Three such cases are shown in Figure 5, where the computed values of R i R j at t = 200σ 2 /D are plotted against ∆ = j − i. In the case of the lowest barrier, V max = k B T, the plot is similar to Equation (36) in that R i R j almost vanishes for i =j; instead, the MSD (i = j) is greater than in the other two cases, indicating that the particles are diffusing almost freely. The case of the highest barrier with V max = 5k B T resembles the ideal SFD, although a close inspection reveals a slight deviation from Equation (35) as a result of overtaking that occurs at a very small rate.
The intermediate case with V max = 3k B T is interesting. At large distances, the same correlation is observed in the case of V max = 3k B T as in the case of V max = 5k B T (and as in the ideal SFD). In contrast, at ∆ = ±1 and ±2, the correlation in the case of V max = 3k B T is remarkably smaller than that for V max = 5k B T. The decrease in the displacement correlation and the increase in MSD must be attributed to overtaking.
The numerical observation on the effects of a finite V max may be understood intuitively, if R i R j is regarded as representing a nested structure of cages with different radii. From this point of view, the plot for V max = 3k B T in Figure 5 can be interpreted as describing the breakdown of inner cages, while the outer ones persist (at least until t = 200σ 2 /D). To elevate this pictorial idea to a quantitative theory on collective motion in SFD with overtaking, we raise the question: How can we modify Equation (35) allowing for overtaking?

Theory of Displacement Correlation in SFD with Overtaking
To find out how Equation (35) should be modified by overtaking, let us re-examine its derivation process. A crucial step is found in Equation (27) where R i is obtained on the assumption that Ξ i is independent of t. This is the point at which the "no overtaking" rule was enforced [30,44].
Overtaking is incorporated into the theory though temporal changes of Ξ i (t) [44]. The kinematics of overtaking were discussed in Subsection 2.5: as is illustrated in Figure 3, two particles exchange their labels according to Equation (21).
In order to recalculate the displacement correlation, now we need to consider the dynamics of overtaking. In principle, the stochastic dynamics of Ξ i (t) should be determined from "first principles" through Equation (22). As a crude approximation, however, here we adopt a phenomenological description characterized by the parameter ν α , referred to as the overtaking frequency or the hopping rate.
The dynamics of Ξ i (t) is modeled as a simple Markovian process in which the labels of neighboring particles are exchanged at a rate of ν α . More precisely speaking, for every pair (i, j) such that Ξ j (t 1 ) = Ξ i (t 1 ) ± 1, their labels are exchanged according to Equation (21) with a probability of 1 − ν α ∆t in every short time interval, ∆t = t 2 − t 1 . This process is essentially equivalent to what is known as "Amida-kuji" (Amitabha's lottery) or the "ladder lottery" [49,50]. In regard to the dynamics of a single tagged particle, the Amida-kuji process is merely a diffusion on a 1D lattice [49,51] such that This is readily shown by solving the master equation for the one-tag diffusion propagator [51], where P(a, t) = P(0, 0; a, t) denotes the probability that the tagged particle, the 0-th one such that Ξ 0 (0) = 0, carries the label Ξ 0 (t) = a at the time t.
In order to calculate the contribution of overtaking to the displacement correlation R i R j , we need to know the two-body diffusion propagator, P(i, j, 0; a, b, t), which represents the probability that the tagged particles, initially labelled with Ξ i (0) = i and Ξ j (0) = j, are found to carry Ξ i (t) = a and Ξ j (t) = b at the time t. Assuming j − i = ∆ > 0 without loss of generality, we write the master equation for P a,b = P(i, i + ∆, 0; a, b, t) with a = b as ∂ t P a,b = ν α (P a+1,b + P a−1,b + P a,b+1 + P a,b−1 − 4P a,b ) + (δ a+1,b + δ a,b+1 )ν α (P b,a + P a,b ) (39) and P a,a = 0. This is more complicated than Equation (38) but still solvable in a Fourier representation (as is shown in Appendix A), so that two-body correlations of Ξ i (t) are given in terms of the modified Bessel function in the limit of N → ∞. Now, we are prepared for calculation of R i R j . From Equation (26) and where we have defined for the sake of brevity. Following the same line of argument as in the derivation of Equation (28), on the assumption thatψ and Ξ i are uncorrelated, we have The summand can be expanded as again, with the assumption thatψ and Ξ i are uncorrelated. The terms including the exponentials of δΞ i can be evaluated analytically on the basis of solutions to the master equations (38) and (39).
In this way, after some calculation, we obtain the displacement correlation. For i = j, we have where D * = D c * + ν α . Note that the last term, 2ν α t, originates from Equation (37). In the case of i = j, using ϕ( · ) defined in Equation (35), we obtain The last term (∆ = j − i ≥ 1 without loss of generality) needs to be evaluated with the two-body diffusion propagator in the form of Equation (A13) in Appendix A, which yields with I n ( · ) denoting the n-th modified Bessel function. It is easy to verify that, in the limit of ν α → 0, Equations (43) and (44) are reduced to the ideal case in Equation (35). The theoretical prediction in Equations (43)-(45) is compared with a result of our particle simulation in Figure 6. With the barrier height and the density chosen as V max = 3k B T and ρ 0 = N/L = 0.2 σ −1 (N = 10 4 and L = 5Nσ), we calculated R i R j and delineated it for three different values of t. The same rescaled variables were used as in Figure 4a: namely, R i R j /(K √ t) is plotted against θ, with K given immediately below Equation (35). The hopping rate was evaluated numerically and estimated to be ν α = 0.0057 D/σ 2 (see Appendix B), which was used to plot Equations (43)- (45) as theoretical curves in Figure 6. The prediction for the ideal SFD (ν α = 0) in Equation (35) is also included with a broken line.
In regard to the difference between the particle simulation and theory for the ideal SFD, Figure 6 exhibits qualitatively the same behavior as was observed in Figure 5-the difference due to overtaking is remarkable only for small θ and occurs in such a way that, except for the self part (i = j, i.e., the MSD), the numerical values of the displacement correlation are smaller than the prediction for the ideal SFD in Equation (35). This means that the effect of overtaking should manifest itself as a negative correction to R i R j for i = j. In this sense, the present theory modifies Equation (35) in the right direction, as the theoretical curve in Figure 6 predicts smaller values of R i R j for i = j in comparison to Equation (35 Figure 6. Comparison of the theoretical predictions of R i R j with particle-based numerical calculations of SFD with overtaking (V max = 3k B T). The three symbols have the same meaning as in Figure 4 with the axes rescaled also in the same way: The thinnest red line, the thin blue line, and the solid black line represent Equation (44) that allows for overtaking (ν α = 0.0057 D/σ 2 ), evaluated at the values of t corresponding to the three kinds of symbols, namely, at t/(σ 2 /D) = 10, 70, and 200, respectively. The theory without overtaking [ν α = 0, i.e., Equation (35)] is shown by a broken line.

Discussion and Concluding Remarks
We have presented calculations of two-particle displacement correlation, R i R j , in a 1D or quasi-1D system of Brownian particles with repulsive interaction, i.e., in SFD with or without overtaking, as an illustrative model of collective dynamics associated with the cage effect. In the ideal SFD (without overtaking), correlated motion with a diffusively growing length scale, λ = 2 √ D c t, was observed. Subsequently, we studied how this result is modified by overtaking; it was shown both numerically and analytically that the overtaking processes destroy short-range correlations alone, leaving long-range correlations nearly intact. This behavior of R i R j , evidenced in Figures 5 and 6, suggests a spatiotemporally nested structure of cages, such that smaller cages are enclosed in larger cages with longer lifetimes.
The main objective of the present work was to shed theoretical light on overtaking in SFD, by extending an analytical theory of SFD to the case of non-ideal SFD in which overtaking is allowed. The analytical theory is based on the method of the label variable ξ. The Lagrangian correlation of the field ψ links R i R j to the Dean-Kawasaki equation (12) that describes the fluctuating density field, while overtaking is taken into account through δΞ i (t). The main analytical result is represented by Equations (43) and (44). This result is reasonably consistent with the numerical behavior of R i R j in Figure 6. A linear solution to the transformed Dean-Kawasaki equation (30) seems to suffice for the description of the outer cages (long scales). Contrastively, the inner cages are not only affected by overtaking via δΞ i (t) but also subject to nonlinear coupling ofψ, as is suggested by the deviation from the linear theory in Figure 4a.
In spite of the reasonable agreement between the analytical and numerical results, however, there are at least two issues that need to be discussed and probably improved in the future. Firstly, the hopping rate (ν α ) could be obtained from βV max and other parameters in a more first-principle-oriented manner, as opposed to the numerical fitting adopted here. Secondly, the validity of the decoupling approximation in Equation (42) is questionable.
In regard to the first issue, we note that the numerical fit for ν α in Equation (A14) in Appendix B is not simply given by an Arrhenius-like expression, e −βV max [46], but includes a prefactor that depends on ρ 0 in a nontrivial way. Although accurate computation of overtaking is difficult and may require improvement in the numerical scheme, which is outside the scope of the present work, the numerical prefactor is interesting enough to motivate theoretical attempts to explain it. While V max in Equation (2) corresponds to the Helmholtz free energy barrier in quasi-1D systems governed by Equations (1a) and (6) [10,27], the hopping rate is rather related to a barrier in the Gibbs free energy [52]. The ρ 0 -dependent prefactor is also reminiscent of the escape probability in predator-prey problems on a 2D lattice [53]. Theoretical evaluation of the hopping rate (ν α ) in the present system will be quite suggestive in a wider context, such as that of 2D colloidal liquids.
In fact, correspondence between the 2D dynamics of colloids and the non-ideal SFD may deserve serious consideration. Displacement correlations in 2D colloidal liquids [54] were recently calculated analytically with the method of the label variable [21]. A linear analysis of the transformed Dean-Kawasaki equation was found to suffice to explain displacement correlations at larger scales, while a phenomenological correction was needed for behavior at smaller scales. The 2D dynamics involve dilatational and rotational modes; the former modes change the local density and correspond to ψ(ξ, t) in SFD, while the latter (or, to be precise, their short-wave components) allow the particles to escape from the cage and correspond to δΞ i (t) in this sense. More intuitively speaking, the overtaking event in Figure 3b can be understood as a small vortex involving two particles, namely 4 and 5.
The second issue concerns the approximation of treatingψ and δΞ i in Equation (42) as uncorrelated, which seems to have made the behavior of Equations (43) and (44) quantitatively incorrect for shorter distances. While a quantitative test of Equation (44) is already given in Figure 6, the MSD in Equation (43) requires more discussion. The validity of Equation (43) could be checked by way of which can be computed numerically and compared with analytical predictions. It seems, however, that the ν α -dependence of D α is in dispute. Hahn and Kärger [23] asserted D α ∝ ν α , while Mon and Percus [24] claimed D α ∝ ν 1/2 α . If Equation (43) is taken literally, it predicts that the longtime behavior of MSD is dominated by the second term on the right side and makes essentially the same prediction as Hahn and Kärger [23]; unfortunately, it contradicts the numerical results of Mon and Percus [24], at least for a certain range of ν α . On the other hand, the reasoning about the origin of D α ∝ ν 1/2 α by Mon and Percus [24] is unsatisfactory, as it seems to lack connection with the collective motion.
As a possible scenario for reconciliation, we may conjecture that the correlation betweenψ and δΞ i , ignored in the present theory, makes a difference to the short-range behavior of R i R j and in the MSD as its limiting case. It seems physically plausible that the collapse of smaller cages may be influenced by density fluctuations with long wavelengths which correlateψ and δΞ i . If this correlation modifies the term containing S in Equation (43) so as to give with some constant c, then it predicts D α to be a sum of terms proportional to √ ν α and ν α . This form might be consistent at once with Hahn and Kärger [23] and with Mon and Percus [24], depending on the values of parameters.
All these issues originate from the phenomenological treatment of overtaking. Improvement upon the present analysis, going beyond the "Amida-kuji" or random exchange approximation, will need to be grounded on the integral in Equation (22) that gives dΞ i /dt. Its systematic treatment will allow calculation of ν α , and it will also make it possible to take the correlation betweenψ and Ξ i into account. (22) is given, in principle, by a solution of the Dean equation [32] for the single particle density, ρ i = ρ i (x, t), and its flux, Q i :

The integrand in Equation
Technically difficult though it might be, this strategy seems quite natural, as it is consistent with the treatment of ψ based on the Dean-Kawasaki equation (12). Development of this strategy for systematic treatment of overtaking effects on SFD will provide useful insights into cage-breaking events in 2D and 3D colloidal glasses. in particular, for p = 2 we obtain Equation (37). The two-body problem in Equation (39) is more complicated. As a counterpart of ∆ = j − i, Although a general solution could be sought in Fourier representation as P a,b = ∑ lPl,d (t)e iaλ , here we concentrate on the mode with l = 0, introducing Q d = Q(∆, d, t) def = ∑ a P(i, i + ∆, 0; a, a + d, t); (A3) this is sufficient for the present purpose, as the knowledge of Q d will allow us to calculate In terms of Q d , the master equation (39) is rewritten as while Q 0 vanishes identically. The initial condition can be expressed as with ε def = e iµ . Subsequently, evaluation of Equation (A5) for d = ±1 gives which also implies A m + A −m = B m + B −m .
Although the initial condition (A6) seems to suggest A m = B m = ε −∆ , this naive choice does not satisfy Equations (A9). Instead, we assume where C = C(ε) is a polynomial of ε (consisting of terms with non-negative powers only) so that ∑ m C(ε)ε d = 0 for d > 0. Substituting Equation (A10) into Equations (A9) yields which is satisfied if we choose Note that this is indeed a polynomial, since ∆ ≥ 1. Thus Q d is obtained in the form of Equation (A7), with A m and B m given by Equations (A10) and (A12), where ε = e iµ = e 2πim/N . In the limit of N → ∞, the summations in Equation (A7) are evaluated as integrals, expressible in terms of Φ n (s) in Equation (A1). The result reads This result allows us to calculate the moments of δΞ j − δΞ i in Equation (A4). In particular, using the second moment, we can also calculate δΞ i δΞ j as where δΞ 2 i = δΞ 2 j = 2ν α t according to Equation (37). After some rearrangement with the recurrence formula for the modified Bessel function, we arrive at Equation (45).

Appendix B. Numerical Evaluation of the Hopping Rate
Within the approximation of the present work, the hopping rate ν α is regarded as a constant whose value depends on the control parameters of the system. The problem of how to specify ν α might be trivial if the lattice dynamics were adopted instead of the present system with the continuous space, because the probability of overtaking, i.e., the ratio of ν α to the collision frequency, would be given as a part of dynamical rules governing the lattice system. More consideration is required in the present case, because overtaking is not given as an elementary process but determined as a result of the Langevin dynamics regulated by the interaction potential V(r) with r ∈ R.
For the present, we have determined the function ν α = ν α (ρ 0 , V max ) by solving Equation (2) numerically, with V(r) given in Equation (5); for numerical ease, the mass is chosen to be finite so that m/µ : σ 2 /D = 1 : 1. For each run with some specific values of (ρ 0 , V max ), we counted overtaking events that had occurred since the time 0, and plotted the cumulative number of overtaking, divided by N, against t. Thus we obtain a plot analogous to Figure 4(a) in Ref. [25] that can be fitted with a straight line, whose slope gives ν α . The values of ν α , collected in this way as a function of (ρ 0 , V max ), are fitted with an Arrhenius-like expression, with a prefactor that depends both on V max and ρ 0 [44]: where a 0 ≈ 1/2 and a 1 ≈ 1/6. The numerical values of ν α reported here, however, should not be regarded as decisive. We must be cautious about difficulty in precise numerical calculation of rare events such as overtaking. In comparison to other quantities, computation of ν α may be sensitive to the choice of the numerical scheme, and in the present calculation, based on a Verlet-like scheme, the computed values of ν α seem to depend on the ratio m/µ:σ 2 /D. A systematic study of this issue will be expected in the future.