Hydrodynamic energy flux in a many-particle system

Turbulent flows exhibit multi-scale energy cascade, which has been derived using Navier-Stokes and Euler equations Kolmogorov (1941a, b); Frisch (1995). However, a hydrodynamic fluid is a collection of random molecules. Hence, it is imperative to derive energy transfers in a multi-particle system. In this letter, we capture such energy transfer in a multi-particle system using high-resolution molecular dynamics (MD) in two dimensions (2D). We simulate a large hydrodynamic vortex superposed on a thermal background in an isolated system with total energy conserved. The vortex in the system generates an energy cascade from large to small scales, akin to that in hydrodynamic turbulence Kolmogorov (1941b).

According to kinetic theory  (Bellomo and Schiavo, 1997) and statistical mechanics (Huang, 2008; Reif, 1965; Goldstein, 2001), multiparticle systems thermalize or reach equilibrium after a sufficiently long time, whose behavior is described by thermodynamic laws (Grad, 1958; Loeb, 2004). A popular multiparticle system is Lenard-Jones (LJ) gas, where particles interact through the short-range LJ potential. In this letter, we demonstrate that an LJ gas with vortex exhibits energy flux and thermalization after the hydrodynamic energy is fully converted to thermal energy. This scenario resembles thermalization in three-dimensional (3D) Euler turbulence Cichowlas et al. (2005). This is the first ab initio or first-principle MD simulations exhibiting turbulent energy flux.

Euler turbulence is studied using an inviscid, incompressible fluid equation given as:

where ${\bf u},p$ are the velocity and pressure fields respectively. An Euler flow conserves total kinetic energy ($\int d{\bf x}u^{2}/2$). In addition, 2D and 3D Euler equations conserve, respectively, the total enstrophy ($\int d{\bf x}\omega^{2}/2$) and total kinetic helicity ($\int d{\bf x}({\bf u}\cdot\boldsymbol{\omega})$), where $\boldsymbol{\omega}=\nabla\times{\bf u}$. Lee (Lee, 1952) and Kraichnan (Kraichnan, 1973a) showed that Euler turbulence admits equilibrium solutions that exhibit approximately $k^{2}$ energy spectrum for 3D and $k$ energy spectrum for 2D (also see Onsager (1949)). For the 3D Euler turbulence, Cichowlas et al. (Cichowlas et al., 2005) showed the Taylor-Green vortex evolves from $k^{-5/3}$ spectrum with positive energy flux at the small and intermediate wavenumbers and $k^{2}$ spectrum with zero energy flux at large wavenumbers. Eventually, the system thermalizes and exhibits $k^{2}$ spectrum for all wavenumbers. Note that $k^{-5/3}$ spectrum and constant energy flux are the predictions by Kolmogorov Kolmogorov (1941a, b); Frisch (1995) for fully-developed hydrodynamic turbulence.
For two-dimensional Euler turbulence, some researchers advocate it reaching equilibrium asymptotically  (Robert and Sommeria, 1991; Kraichnan, 1973b; Lee, 1952; Onsager, 1949), whereas others argue it to be out of equilibrium in the final state (Pakter and Levin, 2018; Bouchet and Simonnet, 2009; Seyler et al., 1975; Verma and Chatterjee, 2022). Further, the initial condition of 2D Euler turbulence affects the asymptotic states considerably; in some cases, the asymptotic states are in equilibrium, but in many cases, they are out of equilibrium. (Bouchet and Simonnet, 2009; Bouchet and Venaille, 2012; Verma et al., 2020; Verma and Chatterjee, 2022).
Now, turning to particle simulations, many papers and books describe various properties of fluids in equilibrium or near-equilibrium using the Lattice Boltzmann Method (LBM), Direct Simulation Monte Carlo (DSMC), and Molecular Dynamics (MD). But, in this letter, we focus on the connections between turbulence and thermal particles. In past LBM and DSMC simulations, Gallis et al., McMullen et al., Martinez et al., Chen and Doolen, Bird (Martinez et al., 1994; Chen and Doolen, 1998; Bird, 1976; Gallis et al., 2017; McMullen et al., 2022) capture the physics at a granular resolution, in particular thermal fluctuations and hydrodynamic features, e.g., instabilities (Kadau et al., 2004; Ashwin and Sen, 2015a; Wani et al., 2022), waves (Tiwari et al., 2015). In particular, Gallis et al. (Gallis et al., 2017) and McMullen et al. (McMullen et al., 2022) simulated 3D decaying hydrodynamic turbulence using DSMC and showed that the resulting spectrum is a combination of $k^{-5/3}$ and exponential. Here, the energy of the structures cascades to small scales and finally to the thermal part, akin to those for 3D Euler turbulence.
Several important works attempt to understand fluid’s macroscopic and microscopic features simultaneously by coupling the Navier-Stokes equation to fluctuating hydrodynamics. In this framework, Bell et al. (Bell et al., 2022) and Bandak et al.(Bandak et al., 2022) showed that such systems exhibit Kolmogorov’s $k^{-5/3}$ and $\exp(-k\eta)$ ($\eta$ is the Kolmogorov’s length scale) energy spectra in the inertial-dissipation range, but it transitions to $k^{2}$ energy spectrum after the exponential part. The $k^{2}$ spectrum is attributed to the thermalization processes at the small scales. Betchov (Betchov, 1957) reported a similar crossover to the $k^{2}$ spectrum in his experiment on the turbulent jet. We remark that the connections between hydrodynamics and kinetic theory shed light on thermalization, the arrow of time, measure of entropy for hydrodynamics, etc. Frisch (1995); Bandak et al. (2022); Gallis et al. (2017); Verma (2019a, 2020); Verma and Chatterjee (2022).
In this letter, we perform molecular dynamics (MD) simulations of a large number of weakly interacting particles in a 2D box. For the initial condition, we chose a large hydrodynamic vortex embedded in a background of randomly moving thermal particles. We isolate this system from the heat bath to conserve the total kinetic energy. Hence, our system is similar to a noisy Euler flow whose total energy is conserved. We observe that our system starts with a nonequilibrium energy transfer from large scales to small scales, and it thermalizes after tens of eddy turnover times. Our MD simulation captures the dynamics more accurately than the DSMC and represents real gas density and the intrinsic transport properties of the medium. Our work demonstrates energy flux and thermalization in an LJ gas that includes realistic interactions, gas density, and intrinsic transport properties of the medium.

In this letter, we report the results of an MD simulation with $2.5\times 10^{5}$ particles using the open-source code LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) (Plimpton, 1995). We employ MPI-based parallel LAMMPS with a maximum of 400 CPU cores. The particles, confined in a periodic 2D (x-y geometry) square box, interact with each other via pairwise LJ potential of the form $\phi_{ij}(r_{ij})=4\epsilon\left[\left(\sigma/r_{ij}\right)^{12}-\left(\sigma/r_{ij}\right)^{6}\right]$, where $r_{ij}$ is the separation between the particles, $\sigma$ is the distance at which $\phi_{ij}$ is zero, and $\epsilon$ is the depth of the potential. In this letter, we choose a small $\sigma$ and isolate the system from the heat bath. Also, the kinetic energy dominates the potential energy for our chosen temperature.

In the present simulation, we employ the standard LJ units, in which $\sigma$ is the unit of length, and $\epsilon$ is the unit of energy. For our simulation, the box size $L$ is 500 LJ units, and the mean free path is $0.9$ LJ units, significantly shorter than the box size. Hence, our simulation satisfies the continuum approximation and effectively includes the impact of thermal fluctuations at hydrodynamic scales.
We start with randomly moving particles (white noise) whose rms value is $u_{\mathrm{th}}=\sqrt{2T}$, where $T$ is the temperature in LJ units (here, the particle mass $m=1$, and the Boltzmann constant $k_{B}=1$). In addition, we provide each particle with a coherent velocity given by ${\bf u(x},t=0)=A\left[\sin(k_{0}x)\cos(k_{0}y)\hat{x}-\cos(k_{0}x)\sin(k_{0}y)\hat{y}\right]$ where $k_{0}=\pi/L$, and $A=u_{\mathrm{th}}$ is a constant amplitude chosen so that the initial coherent and incoherent kinetic energies have similar magnitudes. Figure 1 illustrates the initial condition and its evolution. The above initial condition mimics realistic terrestrial flows that have large-scale macroscopic velocity and small-scale thermal velocity. For present simulations, we take $T=100$ LJ units. Hence, rms velocity $U=\sqrt{200}\approx 14.1$ LJ units, and the eddy turnover time at $t=0$ is $L/U\approx 35.5$ LJ units. Our simulation can be related to inert Ar gas in a box of size 0.17 $\mu m$ whose average inter-particle separation is 0.19 nm and the associated mean free path is 0.304 nm. We also remark that our numerical results are independent of particle numbers if it is more than tens of thousands.

After applying the vortex velocity profile, we let the LJ gas evolve from $t=0$ to $2000$ LJ units, which is approximately 56 eddy turnover time. In Figure 1, we illustrate three snapshots of the flow using the velocity and vorticity fields at $t=0$, 1000, and 2000 LJ time units. These plots illustrate the evolution from an ordered state to a thermalized state, consistent with the numerical results of 3D Euler turbulence Cichowlas et al. (2005); Verma and Chatterjee (2022). We quantify the thermalization process using the kinetic energy spectrum and flux related to the hydrodynamic velocity computed on a 40x40 grid. For a given snapshot, we compute the coarse-grained flow velocity ${\bf u(x},t)$ by averaging over several hundred particles in each grid. We Fourier transform ${\bf u(x},t)$ and compute ${\bf u(k},t)$, and subsequently, the modal energy $E({\bf k},t)=|{\bf u(k},t)|^{2}/2$. After this, we compute the 1D shell spectrum $E(k,t)=\sum_{k-1\leq k^{\prime}\leq k}E({\bf k}^{\prime},t)$ (Frisch, 1995). Our spectral plots in Fig. 3(a,b) exclude $E(k,t)$ at $k={k_{0},2k_{0}}$. However, all the available wavevectors are considered for the computation of nonlinear transfer, energy fluxes, and enstrophy.
We also note that the temperature of the fluid gradually rises as time progresses, as shown in fig. 2. This is because the prominent flow patterns gradually decline over time, with the flow reaching thermal equilibrium approximately within 50 time units or one eddy turnover time. This behavior arises due to the low viscosity of the flow and a high Reynolds number, approximately 812. Consequently, the temperature of the flow undergoes a gradual increase.

Finally, we summarize our results. We study hydrodynamic features such as turbulent energy flux and thermalization in a many-particle system. We isolate our system from the heat bath so as to conserve the total energy and start it with a hydrodynamic vortex. In the early phase, the system exhibits $E(k)\propto k^{-3}$ at small wavenumbers and $E(k)\propto k$ at large wavenumbers. The system gets fully thermalized at $t\gtrsim 500$ LJ units, at which point the coherent hydrodynamic energy is fully transferred to the incoherent hydrodynamic component ($E(k)\propto k$) and thermal noise. We strengthen our arguments using the energy and enstrophy transfers. Our MD simulation clearly demonstrates how large-scale kinetic energy cascades to small scales in a realistic turbulent system.

Our results are consistent with the route to thermalization observed earlier in Euler turbulence Cichowlas et al. (2005), in Navier-Stokes equation with fluctuating hydrodynamics Bell et al. (2022); Bandak et al. (2022), and in DSMC simulations McMullen et al. (2022). However, there are some subtle differences. In our simulation, some of the hydrodynamic energy is transferred to the thermal particles at the microscopic level, which is absent in 3D Euler turbulence. We also remark that particle dynamics is more realistic in MD simulations than in DSMC.

In our system, the large-scale initial structure generates $k^{-3}$ energy spectrum, as well as forward energy and enstrophy fluxes. However, the vortex structure in our system breaks down and thermalizes due to a lack of external forcing. We also note that our system exhibits near-incompressibility, with maximum and average density fluctuations of $0.15$ and $0.05$, respectively. We also observe similar thermalization for other initial conditions, such as sheared flow, wave excitation, Rayleigh-Taylor instability, and four-vortex during the validation process. Interestingly, in earlier simulations of 2D Euler systems with a set of vortices as initial conditions are out of equilibrium Bouchet and Simonnet (2009); Bouchet and Venaille (2012); Verma and Chatterjee (2022). Resolution of this issue needs further investigation.

Thus, our MD simulations of hydrodynamic structures embedded in a noisy environment provide valuable insights into the thermalization process in a realistic fluid system. The work also sheds light on the connection between hydrodynamics and thermodynamics, which are active at macroscopic and microscopic scales, respectively.

The author thanks Jörg Schumacher, Rodion Stepanov, and Soumyadeep Chatterjee for valuable comments. The authors acknowledge the use of AGASTYA HPC for present studies. RW and ST also acknowledge the support for this work through SERB Grant No. CRG/2020/003653, and MKV, the support of SERB Grant Nos. SERB/PHY/20215225 and SERB/PHY/2021473.