Molecular hydrodynamics: Vortex formation and sound wave propagation Free

In principle, the hydrodynamics behavior of a fluid is determined by the dynamics of the molecules constituting the considered fluid. In practice, however, direct molecular dynamics (MD) simulations of hydrodynamic phenomena are severely limited in time and become prohibitive as the system size goes beyond the micrometer scale. The difficulty originates from the large scale difference between MD and hydrodynamics. To resolve the collisions between fluid molecules, the time step and the shortest length scales of MD simulations must resolve the mean collision time and the size of a fluid molecule, respectively. On the other hand, hydrodynamics describes the collective motion of fluid molecules, which corresponds to the zero-wavenumber limit. In this limit, density fields of the globally conserved quantities such as mass, energy, and momentum are only considered.^1–4 Nevertheless, the MD simulation technique^5–7 represents an indispensable tool to investigate the molecular origins of the hydrodynamics theory from first principles. Powerful computing capacities available now are decisive in bridging the gap between the continuum-based analysis of transport phenomena and their modeling on the microscopic level.^8–10 

There have been various attempts to understand the hydrodynamic description from the molecular viewpoint by means of the MD simulation technique. The first approach was to compute transport coefficients appearing in the phenomenological hydrodynamic description. It is based on the Green–Kubo relations, where each transport coefficient is expressed as the time integral of the autocorrelation function of a corresponding dynamical variable. While this approach provides a practical methodology for estimating transport coefficients,^5–7 it assumes the validity of the phenomenological hydrodynamic description and the complete scale separation between MD and hydrodynamics.

One question raised in various MD studies is whether the hydrodynamic description still holds down to the molecular scale, where the validity of the continuum description becomes questionable. To this end, some well-known results from hydrodynamics have been tested in a molecular setting. For example, the applicability of the Stokes–Einstein relation for a molecular-sized tracer particle was extensively investigated. While overall good agreement of MD simulation results with the Stokes–Einstein relation was reported,^11–13 this does not give a definite answer to the original question due to the subtlety existing in determining the particle-solvent boundary condition and the hydrodynamic radius.

Moreover, specifically nonlinear aspects of fluid dynamics have been studied by means of MD simulations such as the formation of eddies behind an obstacle,^14–16 structure formation and flow instabilities,^17–21 flow of immiscible fluids,^22,23 and turbulent mixing,²⁴ to name but a few. The emerging flow patterns suggest that even at the considered molecular length scales, the respective systems can be treated as continuum fluids. Points of instability, such as bifurcation points, however, are prone to large fluctuations seen in the MD results rendering a quantitative comparison with the hydrodynamic description difficult.^8,25

Another direction of investigation is to understand the time correlation functions of a molecular fluid at equilibrium by using the hydrodynamic theory. Since these quantities are directly related to the transport properties, one can investigate the continuous transition from the molecular time scale to the hydrodynamic time scale. The best-known and most extensively studied quantity is the velocity autocorrelation function (VACF) of a tagged particle in a d-dimensional simple fluid (d = 2, 3). Alder and Wainwright observed an algebraic decay proportional to t^−d/2 and showed that this slow decay is caused by a vortex flow forming around the particle.²⁶ Several types of theoretical approaches, having a variety of theoretical perspectives and mathematical sophistications, have been applied and refined to account for this unexpected behavior. For example, the mode-coupling theory²⁷ relates the algebraic decay of the VACF to the vorticity diffusion based on the hydrodynamic description, whereas the kinetic theory²⁸ sees its origin as correlated binary collisions (i.e., ring collisions) in the microscopic dynamics. Nevertheless, all approaches give the identical expression for the long-time decay. Recent extensive MD studies have confirmed that the emergence of the long-time tail is universal for fluid systems.^29–32 

The validity of the assumptions inherent in the derivation of governing equations of hydrodynamics and their limitations has not been systematically scrutinized yet by MD simulations. In particular, no general rules are known indicating the actual ranges of validity of the common restriction to large temporal and spatial scales for specific characteristics, say, of a flow pattern. Also, in general, one does not know precisely when deviations from the hydrodynamic picture must be expected and of which nature and magnitude they would be. A further fundamental postulate of hydrodynamics is that of local equilibrium. It presupposes the separation of temporal and spatial scales.

The flux correlation functions determining the transport coefficients via the Green–Kubo relations, such as viscosity, heat conduction, and diffusion coefficients, typically display a rapid decay at short time scales characteristic of the molecular motion and a slowly decaying long-time tail resulting from relatively large scale spatial structures of hydrodynamic nature. While in three dimensions the total decay rate, including the hydrodynamically caused slow contribution, is sufficiently fast to yield finite, time-independent transport coefficients, the relevant hydrodynamic patterns in a two-dimensional fluid are more persistent. Their decay is so slow that, for example, the diffusion coefficient characterizing the spreading of a tagged particle diverges.

In this paper, we study the hydrodynamic description of collective modes in a molecular fluid system by examining the mode-coupling approach.³³ To this end, we performed an MD simulation study of a two-dimensional simple fluid by computing relevant hydrodynamic modes directly from MD simulations. In order to determine the flow pattern that is generated by the motion of a tagged particle, we select those initial configurations from the equilibrium distribution for which the tagged particle sits in the origin and moves in the positive x-direction. Letting all particles interact via short-range pairwise potentials, a flow pattern emerges that is constructed from the positions and velocities of all particles resulting from the MD simulation. The resulting flow pattern can be decomposed into a potential and a solenoidal flow by means of the Helmholtz decomposition. The latter flow corresponds to a vortex pattern as it was predicted by Alder and Wainwright²⁶ and confirmed in Refs. 34 and 35. The two flow fields are compared to the solutions of the linearized Navier–Stokes (LNS) equations containing time-dependent transport coefficients resulting from Green–Kubo-like formulas with the actual time as the upper integration limit. The spatial resolution of the MD-based flow patterns is at 10% of the fluid particle diameter, considerably finer than the 70% resolution used in the recent work.³⁴ Such high resolution is needed in order to obtain a reliable separation into the potential and the solenoidal flow patterns. This though requires a large ensemble of MD data in order to keep the statistical errors sufficiently small. For efficient computation, we employed not only the ensemble average (over independent samples) but also an average over all particles and a large number of instants of time.

The rest of this paper is organized as follows. In Sec. II, the necessary theory is reviewed. In Sec. III, our MD model and the numerical procedures are introduced. In Sec. IV, the MD simulation results are presented and analyzed. Section V concludes the paper with a summary and an outlook.

From the viewpoint of the mode-coupling approach, to obtain the VACF of a tagged particle, two field variables, the fluid velocity and the distribution of the tagged particle, are essential. While a similar argument can be found elsewhere (e.g., Refs. 30 and 36), we derive a refined expression (17) for the long-time VACF by assuming time-dependent transport coefficients. This leads to a substantially improved agreement with MD simulation results for two-dimensional fluids. On the other hand, with this refined approach, one recovers the well-known expression³³ for the case of constant coefficients. In Sec. II A, we derive relations between the VACF and the average velocity of the tagged particle conditioned on its initial velocity. In Sec. II B, we obtain an approximated form of the latter in terms of a particular fluid velocity field and the tagged particle distribution. In Sec. II C, we obtain the VACF with the help of the hydrodynamic forms of these fields.

We consider a tagged particle suspended in an isotropic fluid in d dimensions; for a detailed description on the molecular setting, see Sec. III A. By denoting the velocity of the tagged particle at the time t by $v$(t), we can express the VACF, C(t) = ⟨$v$(0)·$v$(t)⟩, in terms of the average velocity at time t, ⟨$v$(t)|$v$₀⟩, conditioned on the initial velocity $v$(0) ≡ $v0$⁠, as

whereby the outer average is determined by the Maxwell–Boltzmann velocity distribution Φ(⁠$v0$⁠).

Because of the isotropy of the velocity in equilibrium, the conditional average points into the direction of the initial velocity, i.e., $⟨v(t)|v0⟩=v0|v0|f(t,|v0|)$⁠, where f(t, |$v0$|) is a scalar function of time and of the absolute value of the initial velocity. For sufficiently small initial velocities, this scalar function becomes linearly proportional to |$v0$|, and hence one finds

In order to relate the velocity $v$(t) of a tagged particle to the velocity field of the fluid, we assume that the velocity of the tagged particle agrees with the velocity of the fluid at its actual position x(t). In other words, we assume that $v(t;v0)=du(x(t),t)$⁠, where $u(x,t)$ is the fluctuating velocity field of the fluid. Here, $v$(t; $v0$⁠) denotes a realization of the velocity of a tagged particle with initial value $v$(0; $v0$⁠) = $v0$ and $=d$ indicates equality in distribution. The conditional velocity average ⟨$v$(t)|$v0$⟩ can be expressed in terms of the fluctuating velocity field as

and further transformed as

On going to the second line, we assumed that the velocity field and the tagged particle density are uncorrelated with each other. Their respective averages are denoted as

In order to avoid a very clumsy notation, we suppressed the dependence of the average fields u(x, t; $v0$⁠) and n^tag(x, t; $v0$⁠) on the initial position x(0) ≡ x₀ of the tagged particle. Without restriction, we may choose x₀ = 0 by using a proper coordinate system. Further we note that as a scalar density of a tagged particle in an otherwise isotropic fluid in equilibrium, n^tag(x, t; $v0$⁠) can only depend on the scalar quantities |x|², |$v0$|², and $x⋅v02$⁠. For sufficiently small velocities $v0$⁠, a possible dependence on the last two invariants can be neglected because both are quadratic in the velocity and hence

Combining (1), (4), and (7), one obtains, for the VACF,

For sufficiently small velocities $v0$⁠, the absolute value of the velocity field u(x, t; $v0$⁠) is proportional to the absolute value $v0$ ≡ |$v0$|. This allows one to express the velocity field by its average over all absolute values of the initial velocity as

where the bar, $•¯≡∫dv0φ(v0)•$⁠, indicates a thermal average over the absolute values $v0$ with respect to the Maxwell–Boltzmann distribution $φ(v0)=mv0/(kBT)exp(−mv02/(2kBT))$ for d = 2 with m as the mass of the tagged particle, k_B as the Boltzmann constant, and T as the temperature. Note that $u(x,t;v0)¯=u(x,t;v0ev0)¯$ only depends on the direction of $v0$ denoted by $ev0≡v0/v0$ but not on $v0$⁠. Therefore, upon replacing the velocity u(x, t; $v0$⁠) in (8) by the right-hand side of (9), one can perform the integral over the absolute value of the initial velocity and find

Here, due to the isotropy of the fluid, $ev0⋅u(x,t;v0ev0)¯$ is independent of the orientation $ev0$⁠. In Sec. IV, we set $ev0$ to the standard unit vector e_x along the x-axis and compute $u(x,t)≡u(x,t;v0ex)¯$ as well as $ntag(x,t)≡ntag(x,t;v0ex)¯$ from MD simulations.

Within the hydrodynamic description, the spreading of the tagged particle density is described by a diffusion equation of the form

where, without loss of generality, we choose the origin for the starting point of the tagged particle, i.e., x₀ = 0. As already mentioned above, in order to properly account for the effects of the long-time tails, we allow for a time-dependent diffusion coefficient, which is directly related to the VACF by $D(t)=1d∫0tC(t′)dt′$⁠. This approach will also allow us to consider the dynamics at short molecular time scales. Molecular expressions for D(t) and other transport coefficients are presented in  Appendix B. The solution of the diffusion equation on a square with side length L and periodic boundary conditions can be conveniently given for the Fourier transform of the tagged particle density, ñ^tag(k, t) = ∫ dxn^tag(x, t)e^−ik·x, where k = 2πn/L with n = (n₁, n₂, …, n_d) for n₁, …, n_d = 0, ±1, ±2, …. It becomes

According to the hydrodynamic description, the velocity field u(x, t; $v0$⁠) is given by the solution of the LNS equations (see  Appendix A) subject to the initial condition

where $n¯$ is the mean number density of the molecular fluid. As for the MD simulations, we assumed that the tagged particle exciting the velocity field is initially located at x₀ = 0 and moves with velocity $v0$ relative to the resting fluid at this moment. The resulting solution can be split into the divergence-free component u^⊥ and the rotation-free component u^∥ according to the Helmholtz decomposition,

The rotation-free contribution u^∥ can be represented as a superposition of the sound- and density-modes of the LNS (see  Appendix A) and hence propagates as an attenuated wave with sound velocity. The heat-mode, which is also a longitudinal mode of the LNS, does not contribute because of the uniformity of the temperature field. On the other hand, the divergence-free field u^⊥ coincides with the transverse mode of the LNS. It spreads only diffusively,

The solution evolving from the initial condition (13) is readily obtained in the k-space as

Because the divergence-free part spreads diffusively, it decays much slower than the rotation-free part, which, as just mentioned, propagates with the sound velocity. Hence, the latter can be neglected as a contribution to expression (8) for the VACF at large times.^1,27 Finally, using Parseval’s theorem, one may express the spatial integral on the right-hand side of (8) by means of a sum over all k values, which can be approximated by an integral in the limit of large systems. One then obtains the expression

To summarize, this result is based on the following assumptions: validity of the LNS equations with time-dependent transport coefficients and also large times and large system size. Note that for time-independent transport coefficients D and ν, the well-known algebraic decay expression for the VACF³³ is recovered.

We consider a standard model of a two-dimensional molecular fluid. It consists of N identical fluid particles with mass m in a square of side length L with periodic boundary conditions. The fluid particles interact pairwise via a potential function V(r) of the inter-particle distance r. Hence, the Hamiltonian of the system is given as

where x_i and p_i = $mvi$ are the position and momentum of the ith fluid particle, respectively. For the pair potential V(r), we employ the Weeks–Chandler–Andersen potential, which is a purely repulsive potential of Lennard-Jones type,

Here, σ is the diameter of a fluid particle and ε is the interaction strength parameter. Based on the molecular parameters in (18) and (19), we use reduced (dimensionless) MD units. That is, the units of mass, length, and energy are set to m, σ, and ε, respectively, and the units of any other quantities are determined from them (e.g., $σm/ε$ for time and ε/k_B for temperature).

We simulate a fluid having number density $n¯=0.6$ and temperature $T¯=1$⁠. This state is chosen so that the algebraic decay of the VACF can be readily observed. While this long-time behavior is universal, the time scale on which it is clearly seen depends largely on the density of the fluid.^29,30 At lower densities, the algebraic decay emerges at later times. At higher densities, on the other hand, the VACF displays negative values at short times (due to backscattering effects) before the algebraic decay pattern develops. Moreover, since it takes less time for sound waves to travel across the domain, the disturbance of the long-time tail occurs at earlier times. To identify finite-size effects, we simulate three system sizes having N = 512, 1024, and 2048 fluid particles, which correspond to domain sizes $L=N/n¯=29.2$⁠, 41.3, and 58.4, respectively. The mean free path and the mean collision time are roughly estimated to be 0.2 and 0.2, respectively, from a corresponding hard-disk system having the same number density and temperature.³⁷ 

We perform NV E simulations⁷ using the standard velocity Verlet algorithm with time step size Δt_MD = 10⁻³. Equilibrium samples are prepared as follows. For each initial configuration where the positions and momenta of the fluid particles are randomly chosen with zero total momentum, velocity scaling with the target temperature $T¯=1$ is performed 10 times every 10³ steps and then equilibration is performed for additional 10⁵ steps. While performing a single run for these simulations is not computationally expensive at all, we generate a large NV E-ensemble of $N=1024$ MD trajectories in order to obtain smooth flow patterns.

The number density n(x, t) and the velocity u(x, t) of a fluid, generated by the motion of a tagged particle, together with the resulting tagged particle density n^tag(x, t) are expressed as averages over MD trajectories of the individual fluid particles in terms of the formal expressions

Here, particle 1 is considered as the tagged particle. The double brackets ⟨⟨ ⟩⟩ denote the thermal equilibrium average conditioned on the specific initial position and initial velocity-direction of the tagged particle at x₁(0) = 0 and $v1$(0)/|$v1$(0)| = e_x, respectively.

In order to comply with this initial condition in those cases when the actual initial data of particle 1 deviate from the required values, we used the homogeneity of the configuration space and also assumed its homogeneity by applying to each particle the translation T that shifts particle 1 to the origin (Tx₁(0) = 0) and a rotation R that brings its velocity into the positive x-direction (⁠$Rv1$(0) = |$v1$(0)|e_x), yielding transformed positions (⁠$x^i(t)=RTxi(t)$⁠) and velocities (⁠$v^i(t)=Rvi(t)$⁠) of all particles. Strictly speaking, the configuration space underlying the MD simulations, which is a two-dimensional torus due to the periodic boundary conditions, does not have rotations as symmetry operations because they fail to be bijective. However, as long as the flow pattern has not covered the full torus but is rather concentrated in a region around the origin, one may rotate this pattern as if it were defined on the entire Euclidean plane.

Due to the permutation symmetry with respect to the numbering of particles, any other particle j can also be assigned as the tagged particle and the resulting set of trajectories can be used to perform the averages in (20). Moreover, the time-translational invariance of the trajectories in thermal equilibrium allows one to take any time point lΔt as initial time, yielding new trajectories

where k = 0, …, N₂ labels the lth trajectory with l = 0, …, N₁ − N₂, and N₂ ≪ N₁. Here, Δt is an integer multiple of the time step Δt_MD and the translation T^(j,l) and rotation R^(j,l) are such that the initial tagged particle position x_j(lΔt) is shifted to the origin and its velocity $vj$(lΔt) is rotated into the positive x-direction. Hence, the brackets in (20) denote an average over a large ensemble generated by $N$ MD runs with all particles at any time generating an initial value of trajectories as specified in (21).

In order to estimate the averaged fields (20), we discretized the state space by introducing square cells of side length Δx. The Dirac delta functions were approximated by means of the indicator function $I(C)$⁠, which is 1 if the condition $C$ is true and 0 otherwise. The three fields in (20) were estimated in the following way:

where $⟨•⟩N$ denotes the arithmetic average over the results of $N$ independent MD runs. The average fields were determined using Δx = 0.1, Δt = 500Δt_MD = 0.5, N₁ = 4020 (for N = 512 and 1024) and 2020 (for N = 2048), and N₂ = 20. Thus, the field quantities were computed up to the time N₂Δt = 10 from MD trajectories of length N₁Δt ≈ 2 × 10⁶Δt_MD (for N = 512 and 1024) and 10⁶Δt_MD (for N = 2048). A total of $N=1024$ MD samples were generated.

We determined the longitudinal and transversal components of the velocity field u(x, t) obtained from the MD simulations using the Helmholtz decomposition. In the two-dimensional case, this decomposition is obtained from two scalar fields Φ(x, t) and A(x, t) as³⁸ 

where −$∇$Φ(x, t) is the rotation-free component corresponding to the longitudinal velocity u^∥(x, t) and J$∇$A(x, t) is the divergence-free component corresponding to the transverse velocity u^⊥(x, t). Here, J denotes the counterclockwise rotation by π/2.

Both scalar fields are solutions of the Poisson equation. The source terms of Φ(x, t) and A(x, t) are determined by the divergence and the rotation of the velocity field, respectively, i.e.,

and

where n is the outward normal to the boundary. Because the velocity field does not strictly obey periodic boundary conditions due to the construction of the ensemble on which (22) is based, we used Neumann boundary conditions.

Below we compare the velocity fields and their longitudinal as well as transversal components obtained from the MD simulations with solutions of the LNS equations, which are presented in  Appendix A. The initial condition for the velocity field complies with (13), whereas the fluid density and temperature fields initially are assumed as uniform,

where $n¯$ and $T¯$ are the mean number density and temperature of the system, respectively.

All parameters needed to solve the LNS equations are determined by the MD simulations. As already mentioned, the time-dependent transport coefficients entering the LNS equations are obtained from Green–Kubo-like expressions that are presented in  Appendix B. Thermodynamic variables, such as the adiabatic speed of sound c_s, the ratio of specific heats γ, and the thermal expansion coefficient α, are computed by the method of pressure derivatives;³⁹ see also  Appendix C. As numerical values for the present two-dimensional system, we obtained c_s = 4.43, γ = 1.83, and α = 0.31.

In Sec. IV A, we give a brief description of the field variables that were computed from the MD simulations. In Sec. IV B, we compare these results with the hydrodynamic description. In Sec. IV C, we examine the validity of the assumptions implied by the hydrodynamic description of the long-time VACF and discuss the applicability and limitations of the description.

The velocity fields u(x, t), at time t = 4, obtained from the MD simulations for three system sizes (N = 512, 1024, 2048), are displayed in Fig. 1. In all cases, the velocity field is caused by a tagged particle initially moving in the positive x-direction. The magnitude of the field is obtained from an average over all absolute values of the initial velocity of the tagged particle according to the Maxwell–Boltzmann distribution as introduced in (9). For the small and middle sized systems, the resulting flow patterns cover the entire square, while the flow field of the largest system is still localized around the initial position of the tagged particle. The distribution of the tagged particle n^tag(x, t) is indicated by three contours within which the particle is found with probabilities 0.5, 0.9, and 0.99 when going from the inner most to the outer circle. This distribution is well described by a Gaussian, which is slightly shifted to the right due to the directional bias of the initial velocity of the tagged particle; see also Fig. 8. As expected from the fact that the mass diffusion is much slower than the momentum transport in a dense fluid, n^tag(x, t) spreads slower than u(x, t). The fluid density n(x, t) (not shown) has a background value $n¯$ with an atomistic structure which has the same origin as the peaks found in the pair correlation function.

Whereas the velocity field u(x, t) contains both vorticity and sound-wave contributions, a clear separation into u^⊥(x, t) and u^∥(x, t) is achieved by the Helmholtz decomposition as presented in the upper rows of Fig. 2 (at t = 0.5) and Fig. 3 (at t = 3). The sizes of the patterns represented by the perpendicular and the longitudinal fields clearly reflect the fast wave-like expansion of the latter and the slow diffusional motion of the former. Figure 2 exemplifies the flow patterns at time t = 0.5 corresponding to the two- to threefold of the mean collision time. At this very short time, the velocity field has still very large values close to the center. Yet panel (b) displays a vortex pair that has already developed above and below the center. The largest velocity contributions to u(x, 0.5) though come from the parallel component presented in panel (c). At the somewhat later time t = 3, these large velocity contributions have been already dissipated leaving a vortex and a wave-like contribution as exemplified by Figs. 3(a)–3(c).

When the sound wave has propagated by the distance L/2, the flow pattern obtained from rotated MD configurations starts failing to obey the periodic boundary conditions. The upper limit of time before this happens is given by t_max = L/(2c_s), which corresponds to 3.23, 4.66, and 6.60 for N = 512, 1024, and 2048, respectively. The finite-size effect becomes visible at larger times. Due to the diffusive character, the shear mode propagation and the mass diffusion are slower than the sound wave propagation and finite-size effects on u^⊥(x, t) and n^tag(x, t) only appear at later times.

Even though the identification of the vortex pairs in panels (b) of Figs. 1–3 appears straightforward from an intuitive point of view, we corroborated their existence by two objective criteria; see  Appendix D.

In the second rows of Figs. 2 and 3, the solutions of the LNS equations (A1), together with the corresponding transverse and longitudinal components, are displayed for t = 0.5 and t = 3. While at very short time t = 0.5 large deviations between the MD (depicted in the upper row) and the LNS velocity fields exist, the transverse components (displayed in the middle column) already agree surprisingly well. After the still rather short time t = 3, the agreement between the MD and the LNS velocity fields is very good. The largest deviations can be seen in the parallel field components. The deviations between MD and LNS have two distinct origins. One reason is the continuum nature of the Navier–Stokes equations, which is in manifest contrast to the atomistic structure underlying the MD simulations. An additional reason for the observed differences is the approximation of the nonlinear Navier–Stokes equations by the linearized equations (A1), which is not justified for the considered singular initial condition at least at short times.

In order to get a better understanding to which extent the atomistic structure of the fluid plays a role, we pass from the x-space to the k-space and compare magnitudes of the Fourier transformed fields $ũ•(k,t)=∫dxu•(x,t)e−ik⋅x$ for $•$ = ⊥, ∥ resulting from the MD simulations with the corresponding LNS fields, which are explicitly given by (A2). In Fig. 4, the absolute values $|ũ•(k,t)|=|ũx•|2+|ũy•|21/2$⁠, $•$ = ⊥, ∥, at time t = 0.5 are compared. While, as in direct space, the agreement between the MD and LNS results for the perpendicular contribution is almost perfect, the MD result for the parallel component displays distinct structures at large k values, while the differences at small k values are minor. The most pronounced difference is a ridge along a circle with the radius k = 2π/r₁ corresponding to the nearest-neighbor distance r₁ = 1.096 in the fluid. The appearance of a maximum at this particular k value clearly reflects the influence of the atomistic structure of the fluid on the MD velocity field. Beyond this radius, the absolute velocity $|ũ∥(k,0.5)|$ first decreases and then develops a broad and shallow maximum at even larger k values. The amount of quasi-momentum transferred to small k values is restricted to two regions around 2πn/Le_x with n = ±1 according to both MD and LNS methods. Deviations between the MD and LNS results already at slightly larger k values are clearly visible. Whether they are caused by the atomistic structure of the fluid or by the linearization of the Navier–Stokes equations is not clear.

Accordingly atomistic structures are also visible in direct space at sufficiently small scales. Figure 5 presents the x-component of u^∥(x, t) in the approximately threefold magnified central region at two times t = 0.5 and t = 1. At the smaller time, the MD result displays a shell-like structure, which is determined by concentric circles with radii r_i coinciding with the locations of the first three maxima of the radial distribution function at r₁ = 1.096, r₂ = 2.226, and r₃ = 3.390. In the central circle with radius r₁ and also in the two rings bordered by neighboring radii r_i, there are roughly croissant-shaped regions of forward motion surrounded by regions with backward motion. In total, the amount of backward motion is less pronounced in accordance with the fact that it is generated by backscattering events. These alternating structures of forward and backward motion are much less pronounced at the later time t = 1, when the initial excitation of the fluid by the tagged particle has already propagated away from the center. Of course, none of the atomistic structures are present in the LNS results. Also, the LNS velocity patterns displayed in Fig. 5 appear to be more symmetric with respect to the y-axis than the corresponding MD result.

According to the LNS equations, the centers of the two vortices, defined as the locations where u(x, t) vanishes, move in the directions perpendicular to the initial tagged particle velocity in proportion to $[∫ν(t′)dt′]1/2$⁠; see (A4). In Fig. 6, this prediction is compared with the movement of the respective locations of the perpendicular velocity component determined by MD simulations. The agreement is good as long as finite-size effects can be neglected. The LNS motion of the vortices is strictly confined to the y direction in contrast to the results obtained from the MD simulation according to which the vortices also move to the x-direction approaching the maximal distance x₀ ≈ 0.6. Apart from this minor difference, the agreement of the vortex pattern found in the MD simulations and obtained from LNS is excellent. Likewise, we compared the speed of the wave-like propagation in Fig. 7 and found that both front positions $xmax±$⁠, where |u^∥(x, t)| attains maxima on the x-axis, move according to the MD simulations with the sound speed in complete agreement with the LNS.

We finally compare the distribution of the tagged particle, n^tag(x, t), obtained from the MD simulations with the Gaussian distribution, n^Gauss(x, t), which is the solution of the diffusion equation (11). For a better testing, we correct for the small bias of n^tag(x, t) that emerges from the initial condition of the tagged particle having only velocities in the positive x-direction, by shifting the Gaussian density to the x-direction by the mean value ⟨x(t)⟩ of the tagged particle at the considered time t. The latter is well described by^40,41

This expression is obtained by integrating both sides of (2) up to time t and taking the average over the absolute values of the thermally distributed initial velocities. As can be seen from Fig. 8, the marginal densities n^tag(x, t) ≡ ∫ n^tag(x, y, t)dy according to the MD simulations and the shifted Gaussian density resulting from the diffusion equation (11) agree almost perfectly with each other. In conclusion, the tagged particle density is faithfully described by a Gaussian distribution with variance 2∫ D(t′)dt′ and mean values ⟨y(t)⟩ = 0 and ⟨x(t)⟩ given by (27).

Based on the previous results, we also examine the validity of the analytic expression (17) for the VACF C(t). Its derivation is mainly based on the following three assumptions:

• Assumption 1. According to (8) and (10), C(t) can be expressed in terms of the velocity field u(x, t), which is generated by the motion of the tagged particle relative to the fluid at the initial time, and the probability density n^tag(x, t) to find the tagged particle at a later time t at the position x.

• Assumption 2. The velocity field u(x, t) can be obtained as a solution of the LNS.

• Assumption 3. The velocity field u(x, t) can be substituted by its transverse component u^⊥(x, t).

We note here that Assumption 1 indicates that the tagged particle behaves as any other fluid particle. It still leaves open whether the velocity field u(x, t) is estimated from the positions and velocities of all fluid particles or whether it is obtained from a hydrodynamic consideration as postulated in Assumption 2. Assumption 3 introduces a further simplification which is expected to lead to errors at short times and for finite systems also at large times.

In Fig. 9(a), we compare the VACFs obtained by the MD simulations of systems with three different sizes. The expected t⁻¹ long-time tails emerge first around t = 1 and well describe the VACFs up to the appearance of a series of humps with maxima at times $nx2+ny2L/cs$ with n_x, n_y = 0, 1, 2, …. At these times, a signal propagating with the sound velocity c_s reappears when moving on a square with side length L and periodic boundary conditions. Hence, the humps can be interpreted as the “echoes” of the initial state. Actually, relatively insignificant deviations from the long-time tail behavior appear already before the first hump as was found from MD simulations, which are not presented here, of a much larger system with N = 6.5 × 10⁵ particles.

In the intermediate time interval 5 < t < L/c_s, the analytic expression (17) being indicated by the green line in Fig. 9(a) very well agrees with the MD result for the large system with N = 2048. Yet small but noticeable deviations are present in the time span 1 < t < 5. A much better conformity is achieved with the result calculated with the full solution u(x, t) of the LNS, see  Appendix A, and the Gaussian distribution of the tagged particle (blue line in Fig. 9). Hence, the deviations in the region 1 < t < 5 can be fully attributed to a violation of Assumption 3.

In order to assess the finite-size effects in particular in combination with further approximations, we also considered the finite-size version of (10), which is of the form

Here, a cutoff k_c = π/Δx_c at large k values must be introduced in order to avoid the divergence of C(t) at t = 0.^42,43 We choose as cutoff length Δx_c = 0.1, in agreement with the spatial resolution used for estimation (22) of the MD fields.

In Fig. 9(b), the VACF directly obtained from the MD simulations is compared with the result from (10) with the full velocity field obtained from MD simulations (blue circles) as well as the result from (28) with the LNS u(x, t) of the corresponding finite system with periodic boundary conditions (blue dashed line). The LNS equations are solved with time-dependent transport coefficients based on molecular expressions; see  Appendix B. Whether these transport coefficients and the thermodynamic parameters listed in  Appendix C are determined in finite-size systems or in the limit of an infinite-size system turns out to be insignificant.

While the result from the full MD velocity field perfectly agrees with the direct MD VACF from the shortest to the largest times, the LNS result coincides at short and intermediate times with the infinite system results hence revealing the aforementioned deviations at short times. At large times, the humps caused by echoes on the finite torus are perfectly reproduced by (28), provided that the full velocity u(x, t) is used, whereas with a transversal velocity field u^⊥(x, t), these humps in C(t) disappear, independent of whether the respective MD (green dots) or the LNS velocity field (green dashed line) is used.

In view of the above-formulated assumptions, we may conclude that Assumption 1 is surprisingly well satisfied even in the kinematic region with t < 1, where the LNS predictions fail mainly for the longitudinal components as evident from Figs. 2(c), 2(f), 4(c), 4(d), and 5. Therefore, Assumption 2 cannot be used for times t < 1, which corresponds to less than 5 mean collision times between fluid particles. For t ≳ 1 though, the full LNS velocity field yields a valid description of C(t) for all times, including the echo humps of finite systems; hence, Assumption 2 is justified from the time t = 1 onwards. Assumption 3 sets in to hold for times t ≳ 5. This time scale is with approximately 25 mean collision times still extremely short if considered at the hydrodynamic time scale. Finite-size effects are mainly suppressed with this assumption, which therefore well describes the thermodynamic limit.

In the present study, we investigated self-diffusion in a system of N pairwise interacting soft disks moving in a two-dimensional square with periodic boundary conditions and compared with the results predicted by the linearized Navier–Stokes equations. By computing hydrodynamic field variables directly from extensive MD simulations with high spatial and temporal resolution and by comparing with the solutions of the LNS, we performed a quantitative analysis indicating the regime of validity of the LNS approach both with respect to the length and time scales. In particular, in view of the fact that certain transport coefficients, such as the diffusion coefficient, do not exist for two-dimensional fluids, it was of significant importance to use time-dependent transport coefficients in the LNS. These time-dependent transport coefficients are determined by finite-time Green–Kubo relations. We expect that their use may also improve the agreement of the solutions of the LNS and results from MD simulations for three-dimensional systems at small spatial and short time scales and in this way may narrow the gap between the kinetic and the field theoretic descriptions of a fluid.

In the presently investigated two-dimensional case, we were able to completely close this gap and demonstrate that the hydrodynamic description is valid already after a few mean collision times. The predictions of the LNS set in to hold already at very short times for the perpendicular part of the velocity field, u^⊥(x, t), that contains the vortex pair, eventually being responsible for the long-time tail of the VACF. On the other hand, the longitudinal wave-like part of the velocity field, u^∥(x, t), resulting from the MD simulations contains at short times features at molecular scales that are absent in the respective LNS field. The layered spatial structure of the MD longitudinal velocity field u^∥(x, t) at short times is caused by the typical shell-like structure of the radial distribution function. The asymmetry of the forward and backward wavefronts can be explained by the backscattering events of the tagged particle necessary to build up the backward oriented wave front. However, by comparing the growth of flow patterns, we confirmed that outside the kinetic region the LNS solutions provide accurate quantitative description for both u^⊥(x, t) and u^∥(x, t).

Based on these observations, we revisited the long-time decay of the VACF, which was first observed and explained by Alder and Wainwright and subsequently analyzed by the mode-coupling theory. We investigated the implicit assumptions of the mode-coupling theory and specified the underlying assumptions as well as their applicability. As the central results, we found that the description of the VACF in terms of (i) a fluid velocity field u(x, t) that is conditioned on the initial position and velocity of the tagged particle and (ii) the probability density n^tag(x, t) of the tagged particle quantitatively agrees with the standard definition as the velocity-velocity correlation function of a selected fluid particle from very short times onwards, which are less than the mean collision time. The replacement of the named velocity field by the solution of the LNS equations with properly determined time-dependent transport coefficients is restricted to times that must be of the order of five mean collision times or larger and is hence already valid on quite short time scales compared to standard hydrodynamic times. The use of time-dependent transport coefficients is mandatory because of the divergence of transport coefficients in two-dimensional fluids but may also improve the solutions of three-dimensional fluids at short times.

The deviations of the actual flow patterns from those resulting from LNS at very short times can be traced back to the atomistic structure of the fluid but also might be influenced by the nonlinearity of the full Navier–Stokes equations, which at short times must not be neglected in view of the singular initial condition (13). Finally, the velocity field may be replaced by its transversal part after the relatively short time that it takes until the wave-like propagating part of the velocity does no longer overlap with the tagged particle probability density. The transversal velocity part contains in two dimensions a vortex pair and in three dimensions a vortex ring. These structures are responsible for the algebraically slow decay of the VACF as it already was predicted by Alder and Wainwright.

Expression (17) relating the VACF and the diffusion coefficient to each other in a nonlinear way was recently used by some of the present authors to find the large time behavior of these quantities resulting in proper corrections to the t⁻¹ decay law of C(t) and the conforming diffusion coefficient D(t) for a two-dimensional fluid.⁴⁴ 

A phenomenon closely related to self-diffusion to which the present approach might be applied is Brownian motion.^45,46 In the case of a Brownian particle, the considered particle has different, mostly much larger mass and also a larger volume than a fluid particle. While in the Brownian limit of heavy particles normal diffusion governed by a Markovian Langevin equation characterizes the dynamics well, the crossover behavior in a two-dimensional fluid from anomalous self-diffusion to normal diffusion presents an open problem.

This work was supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics Program under Contract No. DE-AC02-05CH11231 and Korea Advanced Institute of Science and Technology (KAIST), College of Natural Science, Research Enhancement Support Program under Grant No. A0702001005. An award of computer time was provided by the ASCR Leadership Computing Challenge (ALCC) program. This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract No. DE-AC05-00OR22725.

We present the form of the LNS equations considered in the paper and obtain analytic results (A4) for the transverse velocity in the two-dimensional case. For the number density n(x, t), the velocity u(x, t), and the temperature T(x, t), the linearized governing equations of a compressible fluid are written as^1,2,33

where $n¯$ is the mean number density, c_s is the adiabatic speed of sound, γ is the ratio of specific heats, ν is the kinematic viscosity, D_l is the kinematic longitudinal viscosity, α is the thermal expansion coefficient, and D_T is the thermal diffusivity.

The solution of (A1) obeying the initial conditions (13) and (26) is readily expressed in a closed form in the Fourier space. By assuming time-dependent coefficients ν(t), D_l(t), and D_T(t) and $v0$ = $v0$e_x, the longitudinal and transverse velocities $ũ•(k,t)=∫dxe−ik⋅xu•(x,t)$⁠, $•$ = ⊥, ∥, become

where k and $k^$ denote the magnitude and unit vector of k, respectively, and Γ_s(t) = D_l(t) + (γ − 1)D_T(t) is the sound attenuation coefficient.

In the two-dimensional case, from the inverse Fourier transform of (A2a), the transverse velocity has the following closed-form expression:

where r and θ denote the polar coordinates of x and $B(t)=∫0tν(t′)dt′$⁠. From the condition u^⊥ = 0, we can easily find the location of each vortex center (0, ±y₀), where

and ξ ≈ 1.26 is the positive solution of $(ξ+12)e−ξ=12$⁠. For convergent ν(t) (i.e., $limt→∞ν(t)=ν$⁠), we have $∫0tν(t′)dt′≈νt$ at large t and thus obtain the $νt$ time scale of y₀.

We consider the following five time-dependent transport coefficients which are defined by the Green–Kubo relationship in terms of correlation functions:

Then the respective flux is defined as

Here V and $P$ are the volume and the average pressure of the system, respectively, and e_i, f_ij, and r_ij are the energy of atom i, the inter-atomic force, and the position vector between atoms i and j, respectively. Figure 10 displays the MD estimation results of these flux correlation functions and the respective transport coefficients.

The sound attenuation coefficient, defined as

can be obtained from the heat conductivity λ(t), shear viscosity η(t), and bulk viscosity ζ(t), whereby the thermal diffusivity D_T(t) and the longitudinal diffusivity (kinematic longitudinal viscosity) D_l(t) are related to these three coefficients as

where C_P and γ are the isobaric heat capacity and the ratio of the specific heat of the system, respectively.

Here we briefly review Meier and Kabelac’s method³⁹ and apply this method for the calculation of the necessary thermodynamic parameters in two-dimensional fluid systems. Starting from the entropy S = k_B ln Ω, where Ω is the phase-space volume, the adiabatic speed of sound, for example, defined as

can be expressed in terms of derivatives of the phase-space volume with respect to the energy E and the volume V, yielding

where $ω=∂Ω∂E$ is the phase-space density and M and ρ are the total mass (i.e., $M=∑i=1Nmi$⁠) and mass density (i.e., ρ = M/V) of the system, respectively.

Defining the center of mass related quantity G,

we express the phase-space volume depending on the conserved quantities N, V, E, P, G as

where $P=∑i=1NPi$ is the total momentum, H = K + U is the Hamiltonian of the system, K is the kinetic energy of the system, U is the potential energy of the system, Θ is the Heaviside step function, and C_N is a normalization constant rendering the phase-space volume dimensionless.

We define the phase-space functions as $Ωmn=1ω∂m+nΩ∂Em∂Vn$⁠. Expressions of Ω_mn in terms of the kinetic energy K and the volume derivatives of the potential energy ∂ⁿU/∂Vⁿ can be derived by using derivatives of (C4). We list the phase-space functions Ω_mn (m + n ≤ 2) of the two-dimensional NV EPG ensemble in Table I.