Multi-Particle-Collision Simulation of Heat Transfer in Low-Dimensional Fluids

1. Introduction

Simulation is often the only viable tool to study many-particle systems driven away from equilibrium, especially where external mechanical and thermal forces are strong. Molecular dynamics is the most natural approach, but may be computationally expensive. Since one is often interested in large-scale properties, details of microscopic interactions may not be essential, since only the basic conservation laws should matter. It is thus sensible to look for methods based on stochastic processes that may effectively account for molecular interactions. One popular approach is the Multi-Particle Collision (MPC) dynamics, a mesoscale description where individual particles undergo stochastic collisions, rather than genuine Newtonian forces. The implementation was originally proposed by Malevanets and Kapral [1,2] and consists of two distinct stages: a free streaming and a collision one. Collisions occur at fixed discrete time intervals, and space is discretized into cells that define the collision range. The method captures both thermal fluctuations and hydrodynamic interactions.

The MPC dynamics is a useful tool to investigate concrete systems and indeed has been used in the simulation of a variety of problems, like polymers in solution [1], colloidal fluids [3], plasmas and even dense stellar systems [4] etc. Besides its computational convenience it is also a useful approach to address fundamental problems in statistical physics [5,6,7,8] and in particular the effect of external sources.

In this work, we focus on the application of the MPC method to study heat transfer in a simple fluid both at equilibrium and in the presence of external reservoirs. In particular, we analyze the dimensionality effects on thermal transport in mesoscopic and confined fluids. Energy transport in low-dimensional systems has been thoroughly studied in recent decades and is crucial for achieving an understanding of macroscopic irreversible heat transfer on the nanoscale [9,10,11,12]. Also, it serves as a theoretical foundation for thermal energy control and management [13,14,15]. This is even more relevant at the nano- and microscale, where novel effects caused by reduced dimensionality, disorder, and nano-structuring affect natural and artificial materials [16]. One remarkable property is that in low-dimensional many-particle systems, energy propagates super-diffusively, implying a breakdown of classical Fourier’s law. This has been well studied, mostly in lattice systems but much less in fluids. We will show that the MPC dynamics is particularly effective and allows for a very accurate test of existing theories of anomalous transport.

The paper is organized as follows. In Section 2 we recall the basic definitions of the MPC dynamics and the thermal-wall method used to enforce the interaction with external heat baths. We then present results of both equilibrium and nonequilibrium for the one- (1D), two- (2D), and three-dimensional (3D) mesoscopic fluids, Section 3 to Section 5 respectively. In Section 6 we illustrate how, upon changing the aspect ratio of the simulation box, one can observe crossover behavior in size dependence of the transport coefficient. In Section 7, we briefly discuss the case 2D fluids with magnetic field to clarify that pseudomomentum conservation is not the necessary condition for the anomalous heat conduction behavior.

2. Mesoscopic Fluid and Thermal Walls

The 1D, 2D, and 3D mesoscopic fluid models we cosider in this study are shown in Figure 1. The fluid consists of N interacting point particles with equal mass m, confined in a system volume. As is shown: In 1D system, its volume has a length of L in the x direction; in 2D system, its volume also has a width of W in the y direction; in 3D system, its volume has an additional height of H in the z direction.

At the boundaries $x=0$ and $x=L$, the particles interact with a heat bath of temperature $T_h$ or $T_c$; the heat baths are modeled as thermal-walls [17,18]. When a particle crosses the $x=0$ (or $x=L$) boundary, it is reflected back with a new random velocity ($v_x$, $v_y$, and $v_z$ in the x, y, and z directions) assigned by sampling from a given distribution [17,18]:

$$\begin{array}{cc}\hfill f\left(v_x\right)=& {\displaystyle \frac{m|v_x|}{k_B{T}_{\alpha }}}\exp \left(-{\displaystyle \frac{m{v}_x^2}{2k_B{T}_{\alpha }}}\right),\hfill \\ \hfill f\left(v_{y,z}\right)=& \sqrt{{\displaystyle \frac{m}{2\pi k_B{T}_{\alpha }}}}\exp \left(-{\displaystyle \frac{m{v}_{y,z}^2}{2k_B{T}_{\alpha }}}\right),\hfill \end{array}$$

(1)

where $T_{\alpha }$ ($\alpha =h,c$) is the temperature of the heat bath in dimensionless units and $k_B$ is the Boltzmann constant. The particles are subject to periodic boundary conditions in the y and z directions. We point out that the numerical results also apply to fixed boundary conditions since, in both cases, $v_{y,z}>0$ and $v_{y,z}<0$ with equal probability $p=0.5$.

Interaction among particles as prescribed by the MPC method [1,3,19], amounts to first partition the simulation box into many cells of linear size a (as is shown in Figure 1 (b) for 2D case). The dynamics evolves in discrete time steps, each step consisting of a free propagation during a time interval $\tau$ followed by an instantaneous collision event. During propagation, the velocity ${\mathbf{v}}_i$ of a particle is unchanged, and its position is updated as

$${\mathbf{r}}_i\to {\mathbf{r}}_i+\tau {\mathbf{v}}_i.$$

(2)

For all particles in a given cell, their velocities are updated according to the following collision rules:

• In the 1D case, the velocity of the ith particle in the jth cell is changed according to the update rule

  $${\mathbf{v}}_i\to A_j{\omega }_i+B_j,$$

  (3)

  where ${\omega }_i$ is randomly sampled by a thermal distribution at the cell kinetic temperature $T_j$, while $A_j$ and $B_j$ are cell-dependent parameters, determined by the condition of total momentum and total energy conservation in the cell [20].
• In the 2D case, all particles found in the same cell are rotated around the z axis, with respect to their center of mass velocity ${\mathbf{V}}_{\mathrm{c}.\mathrm{m}.}$ by two angles, $\theta$ or $-\theta$, randomly chosen with equal probability. The velocity of the ith particle in a cell is thus updated as

  $${\mathbf{v}}_i\to {\mathbf{V}}_{\mathrm{c}.\mathrm{m}.}+{\widehat{\mathcal{R}}}^{\pm \theta }\left({\mathbf{v}}_i-{\mathbf{V}}_{\mathrm{c}.\mathrm{m}.}\right),$$

  (4)

  where ${\widehat{\mathcal{R}}}^{\pm \theta }$ are the rotation operators by the angle $\pm \theta$.
• In the 3D case, the velocity of the ith particle in a cell is updated as in 2D case, with the difference that the rotation axis is also randomly selected.

The resulting motion conserves the total momentum and energy of the fluid. Note that the angle $\theta =\pi /2$ for 2D and 3D cases corresponds to the most efficient mixing of the particle momenta. Note also that the probability of collision between particles increases as $\tau$ decreases. Thus the time interval $\tau$ can be changed to tune the strength of the interactions that, in turn, will affect the transport properties.

In the nonequilibrium setup, we set $T_{\alpha }$ to be slightly biased from the nominal temperature T, i.e., $T_{h,c}=T\pm \Delta T/2$, to investigate the dependence of the thermal conductivity $\kappa$ on the system length L. In our simulations, each particle is initially given by a random position uniformly distribution and a random velocity generated from the Maxwellian distribution at the temperature T. After the system reaches the steady state, we compute the thermal current J that crosses the system according to its definition (i.e., the average energy exchanged in the unit time and unit area between particles and heat bath) and taking it into Fourier’s law, $\kappa =JL/(T_h-T_c)$, to compute $\kappa$. We thus examine the dependence of $\kappa$ on L, to assess whether the heat conduction behavior of the system is anomalous or normal.

In the integrable case (i.e., when $\tau =\infty$, each particle maintains unchanged velocity as it crosses the system from one heat bath to the other), we can obtain an analytical expression for the thermal conductivity:

$$\kappa =L\left(d+1\right)\sqrt{{\displaystyle \frac{{\rho }^2{k}_B^3}{2\pi m}}}/\left({\displaystyle \frac{1}{\sqrt{T_h}}}+{\displaystyle \frac{1}{\sqrt{T_c}}}\right).$$

(5)

where d is the spatial dimension and $\rho$ is the particle density. It is shown that in the integrable case, transport is ballistic and the thermal conductivity is a linear function of the system length. This analytical result will be used to compare with our simulations as a numerical verification.

In the nonequilibrium setting, the difference between normal and abnormal heat conduction behaviors can be also appreciated upon examining the steady-state kinetic temperature profiles $T\left(x\right)$. In our simulations, $T\left(x\right)$ is measured as described in [21]. For systems with normal heat conduction, $T\left(x\right)$ is determined by solving the stationary heat equation assuming that the thermal conductivity is proportional to $\sqrt{T}$ as prescribed by standard kinetic theory, yielding [22]

$$T\left(x\right)={\left[T_h^{3/2}\left(1-{\displaystyle \frac{x}{L}}\right)+T_c^{3/2}{\displaystyle \frac{x}{L}}\right]}^{2/3}.$$

(6)

this prediction will be used to compare with our simulation results for normal heat conduction. On the other hand, for systems with abnormal heat conduction, $T\left(x\right)$ is expected to be qualitatively different, being solution of a fractional diffusion equation as demonstrated in several examples [23,24]. A typical feature is that the temperature profile is concave upwards in part of the system and concave downwards elsewhere, and this is true even for small temperature differences [25,26,27]. This will be used also to check our numerical simulations for abnormal heat conduction.

To check the results obtained in the nonequilibrium modeling, we will further turn to the comparison with linear-response results obtained in equilibrium modeling. Based on the celebrated Green-Kubo formula, which relates transport coefficients to the current time-correlation functions $C\left(t\right)$, the thermal conductivity can be expressed as [9,10,28]

$${\kappa }_{\mathrm{GK}}={\displaystyle \frac{\rho }{k_B{T}^2}}\underset{{\tau }_{\mathrm{tr}}\to \infty }{lim}\underset{N\to \infty }{lim}{\displaystyle \frac{1}{N}}{\int }_0^{{\tau }_{\mathrm{tr}}}C\left(t\right)dt.$$

(7)

In this formula, $C\left(t\right)\equiv 〈\mathcal{J}\left(0\right)\mathcal{J}\left(t\right)〉$ and $\mathcal{J}\equiv {\displaystyle \frac{1}{2}}{\Sigma }_i^N{\mathbf{v}}_i^2{v}_{x,i}$ represents the total heat current along the x coordinate in the equilibrium state. In the simulations, we consider an isolated fluid with periodic boundary conditions also in the x direction. The initial condition is randomly assigned with the constraints that the total momentum is zero and the total energy corresponds to T. The system is then evolved and after the equilibrium state is attained, we compute $C\left(t\right)$ and the integral in Eq. (7). Usually, the integral is truncated up to ${\tau }_{\mathrm{tr}}=L/v_s$ ($v_s$ is the sound speed) [9,10]. This results the superdiffusive heat transport ${\kappa }_{\mathrm{GK}}\sim L^{1-\lambda }$ as long as $C\left(t\right)$ decays as $\sim t^{-\lambda }$ with $\lambda <1$.

In order to compare $\kappa$ and ${\kappa }_{\mathrm{GK}}$ more accurately, we can resort to the spatiotemporal correlation function of local heat currents to compute the sound speed $v_s$ of the system in Eq. (7). The spatiotemporal correlation function of local heat currents is defined as [29,30,31]

$$C\left(x,t\right)\equiv 〈{\mathcal{J}}^{loc}(0,0){\mathcal{J}}^{loc}(x,t)〉.$$

(8)

Numerically, we compute $C\left(x,t\right)$ as performed in [32]: the system is divided into ${\displaystyle \frac{L}{b}}$ bins in space of equal width $b=0.2$; the local heat current in the kth bin and at time t is defined as ${\mathcal{J}}^{loc}(x,t)\equiv {\sum }_i{\displaystyle \frac{1}{2}}{m}_i{\mathbf{v}}_i^2{v}_{x,i}$, where $x\equiv kb$ and the summation is taken over all particles that reside in the kth bin. It is found that $C\left(x,t\right)$ features a pair of pulses moving oppositely away from $x=0$ at the sound speed [33,34], which are recognized to be the hydrodynamic sound modes. Their moving speeds of the two pulses allows to estimate the sound speed $v_s$ within the fluid.

As for the choice of parameters, in both nonequilibrium and equilibrium settings, we set $T=1$, $\Delta T=0.2$, $m=k_B=a=1$, $\theta =\pi /2$, and $\rho =5$ throughout the paper. In addition, for all data points shown in the figures, the errors are $\le 1\%$.

3. 1D fluid

3.1. Nonequilibrium Results

To start, let us discuss the dependence of the thermal conductivity $\kappa$ on the system length L with different interaction strengths. Here, the time interval $\tau$ between successive collisions will be used to tune the strength of the interactions. The particle mean free path ℓ, in an uniform system, is proportional to the thermal velocity of the fluid $v_T=\sqrt{T}$ and the time of the MPC move, namely $\ell \sim v_T\tau$.

To check the results and provide a numerical example, we first quantify the noninteracting system $\tau =\infty$ . In Figure 2, we report $\kappa$ of 1D case (Eq. (5)) with a red line is compared with our simulations (black cycles). It can be seen that they agree very well with each other. These simulations clearly strongly support our analysis.

We next turn to the interacting systems with $\tau <\infty$: It can be seen in Figure 2 that for weak interactions ($\tau =50,10$), $\kappa$ tends to saturate and becomes constant as L is increased, following the Fourier law. However, it can also be seen that for the strong interactions ($\tau =1,0.1$), $\kappa$ is no longer constant but diverges with L. In particular, for $\tau =0.1$, $\kappa$ eventually approaches the scaling $\kappa \sim L^{1/3}$ like that predicted in 1D momentum-conserving fluids [35]. This is the anomalous heat conduction behavior dominated by hydrodynamic effect, well known in nonequilibrium heat transport.

Altogether this can me understood as a crossover from ballistic, to diffusive (kinetic) and then to anomalous behavior controlled by corresponding timescales, see [36,37,38,39,40,41,42].

The difference between normal and abnormal heat conduction behaviors can be further appreciated also in the steady-state kinetic temperature profiles $T\left(x\right)$. For systems with normal heat conduction, $T\left(x\right)$ is predicted by Eq. (6). In Figure 3(a) this prediction is compared with our simulation results for $\tau =50$. It is seen that there is a good agreement between the results of our numerical simulations and Eq. (6). To better appreciate the deviations from the prediction, we also plot the differences $\delta T$ between the data and the black line. It is shown in the inset of Figure 3(a) that $\left|\delta T\right|$ decreases with increasing L, as expected since $|\nabla T|=\Delta T/L$ decreases when L increases, indicating that the linear response can correctly describe the transport properties of the system for large enough system length. As introduced in the Section 2, for systems with abnormal heat conduction, the typical feature of the temperature profile is that $T\left(x\right)$ is concave upwards in part of the system and concave downwards elsewhere, and this is true even for small temperature differences [25,26,27]. This is confirmed in Figure 3(b) by our numerical simulations for $\tau =0.1$. Note that the data for three different L overlap with each other, implying that the deviations from Fourier’s behavior are not finite-size effects. Altogether, those numerical results again support our findings based on the length-dependence of the thermal conductivity, that heat conduction in the weak interactions is normal while in the stronger interactions it is abnormal.

3.2. Equilibrium Results

To check the results obtained in the 1D nonequilibrium modeling, we now turn to the comparison with the results obtained by the Green-Kubo formula in the equilibrium modeling. The results for the current time-correlation functions $C\left(t\right)$ with different $\tau$ values are presented in Figure 4. It can be seen from Figure 4 (a) and (b) that for $\tau =0.1$ and $\tau =1$, the correlation function eventually attains a power-law decay $C\left(t\right)\sim t^{\gamma }$ with $\gamma =-2/3$, fully compatible with the theoretical prediction of the 1D case [43,44]. Substituting it in Eq. (7), and cutting off the integration as explained above, one obtains the superdiffusive scaling ${\kappa }_{\mathrm{GK}}\sim L^{1/3}$, in agreement with our nonequilibrium modeling.

However, it is clear in Figure 4 (c) and (d) that for $\tau =10$ and $\tau =50$, $C\left(t\right)$ undergoes a rapid decay at short times, and eventually, it begins to oscillate around zero (the negative values of $C\left(t\right)$ are not shown in this log-log scale.). Fitting function with the green solid line exhibits an exponential decay. To compare ${\kappa }_{\mathrm{GK}}$ and $\kappa$ more accurately, we compute the sound speed $v_s$ of the system with the help of the spatiotemporal correlation function $C\left(x,t\right)$ of local heat currents defined in Eq. (8). In Figure 5, we present $C\left(x,t\right)$ for the system size $L=32$. The two peaks representing the sound mode can be clearly identified in Figure 5 (a) and (b). Their moving speed $v_s$ is measured to be $v_s\simeq 2.40$ for $\tau =50$ and $v_s\simeq 1.875$ for $\tau =0.1$. In Figure 2, the horizontal line for $\tau =50$ is obtained truncating the integral in Eq. (7) upto $L/v_s$ with $v_s=2.40$. It can be seen that ${\kappa }_{\mathrm{GK}}\left(L\right)$ agrees with $\kappa$. Thus, the equilibrium simulations are fully consistent with the nonequilibrium simulations.

4. 2D fluid

4.1. Nonequilibrium Results

As for the 1D case, let us first reports the dependence of $\kappa$ on L for different $\tau$ values, see Figure 6. To provide a numerical check, we show that the noninteracting (integrable) expression Eq. (5) (red solid line) matches very accurately the simulation data (black circles). For the interacting systems it can be seen in Figure 6 that for the weak interactions ($\tau =10,5$), $\kappa$ tends to saturate and becomes constant as L is increased. This indicates that the normal heat conduction behavior for the nearly integrable 2D fluid system is also dominated by kinetic effect. However, as $\tau$ decreases, $\kappa$ is no longer constant but diverges with L. In particular, for $\tau =0.1$ and $\tau =0.01$, $\kappa$ eventually approaches the scaling $\kappa \sim ln\left(L\right)$ like that predicted in 2D momentum-conserving systems [11,45]. Therefore the scenario is similar to the 1D case and can be described as a crossover from kinetic to hydrodynamic regimes.

In order to demonstrate the crossover behavior, we assume that the flux can be decomposed as the sum into normal and anomalous contributions [39]

$$J=J_N+J_A,$$

(9)

From kinetic theory we expect that $J_N\left(L\right)=\mathfrak{a}/(\mathfrak{b}+L)$ so we can use the conductivity data for large $\tau$ to estimate $J_N$ and thus deduce $J_A=J-J_N$. In the inset of Figure 6 we show that for the strong interactions ($\tau <0.5$), indeed $J_A$ is proportional to $log\left(L\right)/L$. To our knowledge, this provides one of most convincing numerical evidence of the logarithmic divergence of the conductivity in 2D.

The difference between normal and abnormal heat conduction behaviors can be further verified by $T\left(x\right)$ as in the 1D case. For systems with normal heat conduction $T\left(x\right)$ is predicted by Eq. (6). In Figure 7(a) this prediction is compared with our simulation results for $\tau =10$. As expected, there is a good agreement between the results of our numerical simulations and Eq. (6). However, for $\tau =10$, $T\left(x\right)$ is concave downwards in the left part of the system and concave upwards in the right part of the system. This again conforms to the temperature distribution characteristics of abnormal heat conduction.

4.2. Equilibrium Results

To check the results obtained above, we now turn to the comparison with the results obtained by the Green-Kubo formula in the equilibrium modeling. The results for $C\left(t\right)$ with different $\tau$ values are presented in Figure 8. It can be seen from Figure 8 (a) and (b) that for $\tau =0.01$ and $\tau =0.1$, the correlation function eventually attains a power-law decay $C\left(t\right)\sim t^{\gamma }$ with $\gamma =-1$, fully compatible with the theoretical prediction of the 2D case [11,45]. Taking it in Eq. (7), one will obtain the superdiffusive heat transport ${\kappa }_{\mathrm{GK}}\sim ln\left(L\right)$, in agreement with nonequilibrium data.

However, it is clear in Figure 8 that as $\tau$ further increases from $\tau =0.1$ to $\tau =10$, $C\left(t\right)$ will change from power-law decay to exponential decay. This means that as $\tau$ increases, the kinetic effects will play a dominant role, and thus heat conduction will change from anomalous behavior to normal behavior, as observed in Figure 6.

To compare ${\kappa }_{\mathrm{GK}}$ and $\kappa$ more accurately, we compute the sound speed $v_s$ of the system as performed in 1D case. In Figure 9, we also present $C\left(x,t\right)$ for the system size $L=32$. The two peaks representing the sound mode can be clearly identified in Figure 9 (a) and (b). Their moving speed $v_s$ is measured to be $v_s\simeq 2.05$ for $\tau =10$ and $v_s\simeq 1.50$ for $\tau =0.1$. In Figure 6, the horizontal line for $\tau =10$ is obtained truncating the integral in Eq. (7) upto $L/v_s$ with $v_s=2.05$. It can be seen that ${\kappa }_{\mathrm{GK}}\left(L\right)$ agrees with $\kappa$. Again, the equilibrium simulations are fully consistent with the nonequilibrium simulations.

5. 3D fluid

5.1. Nonequilibrium Results

To complete the study of dimensionality effects, we also performed a series of simulations for the 3D case. The results in Figure 10 demonstrate that even here the formula Eq. (5) accounts very accurately for the non-interacting case (compare the red solid line with black circles). The data for the interacting case (symbols in Figure 10) confirms that, even for the smallest $\tau$ considered the conductivity converges to a finite value. As expected, Fourier’s law holds for 3D fluid system with the momentum conservation.

The normal heat conduction behavior can be further verified by $T\left(x\right)$. For systems with normal heat conduction $T\left(x\right)$ is predicted by Eq. (6). In Figure 11 this prediction is compared with our simulation results for $\tau =1$ and $\tau =0.1$. As expected, there is a good agreement between the results of our numerical simulations and Eq. (6), conforming to the temperature distribution characteristics of normal heat diffusion.

5.2. Equilibrium Results

To check what obtained in the 3D nonequilibrium modeling, we now turn to the comparison with the results obtained by the Green-Kubo formula in the equilibrium modeling. The results for $C\left(t\right)$ with different $\tau$ values are presented in Figure 12. It can be seen from Figure 12 (a) and (b) that for $\tau =0.01$ and $\tau =0.1$, the correlation function eventually attains a power-law decay $C\left(t\right)\sim t^{\gamma }$ with $\gamma =-1.5$, fully compatible with the theoretical prediction of the 3D case [11,45].

However, it is clear in Figure 12 that as $\tau$ further increases from $\tau =0.01$ to $\tau =1.0$, $C\left(t\right)$ will change from power-law decay to exponential decay. This means that as $\tau$ increases, the normal heat conduction behavior observed in Figure 10 will change from being dominated by the hydrodynamic effect to being dominated by the kinetic effect.

To compare ${\kappa }_{\mathrm{GK}}$ and $\kappa$ more accurately, we compute the sound speed $v_s$ of the 3D system. In Figure 13, we also present $C\left(x,t\right)$ for the system size $L=32$. The two peaks representing the sound mode can be clearly identified in Figure 13 (a) and (b). Their moving speed $v_s$ is measured to be $v_s\simeq 1.30$ for $\tau =1$ and $v_s\simeq 1.275$ for $\tau =0.1$. In Figure 10, the horizontal line for $\tau =1$ is obtained truncating the integral in Eq. (7) upto $L/v_s$ with $v_s=1.30$. It can be seen that ${\kappa }_{\mathrm{GK}}\left(L\right)$ agrees with $\kappa$. Once again, the equilibrium simulations are fully consistent with the nonequilibrium simulations.

6. Dimensional Crossovers

Dimensional-crossover is a relevant topic for thermal transport in low-dimensional materials [14]. Indeed, in 2014 it has been experimentally observed in the suspended single-layer graphene [46]. In this experimental setup, the width of the samples is kept fixed and the thermal conductivity changes upon increasing their length is measured. As the length increases, it is natural to expect that a dimensional-crossover behavior from two dimensions to quasi-one dimension will occur. These research results have greatly enriched our understanding of heat conduction in lattice systems.

We show that the MPC approach can be used successfully to investigate this issue, considering 2D and 3D mesoscopic fluid models with fixed transverse sizes ($W=1$ in 2D and $W=H=1$ in 3D) and study how $\kappa$ changes with L. It can be seen from Figure 14 (a) that in 2D fluid models, upon increasing the aspect ratio of the system, for both $\tau =10$ and $\tau =0.1$, $\kappa$ eventually follows the 1D divergence law $\kappa \sim L^{1/3}$ as L increases. This means that both in the case where dominated by the kinetic ($\tau =10$) or hydrodynamic effect ($\tau =0.1$), there exists dimensional-crossover above a given aspect ratio. As is shown in Figure 14 (b), in 3D fluid models, there is also a similar phenomenology. Altogether, these results confirm that the theories developed for the strictly 1D case effectively extend also quasi-1D, provided that the transverse extent of the sample is small enough.

To further support the dimensional-crossover behavior of heat conduction observed above, we perform the equilibrium simulations of $C\left(t\right)$ in 2D and 3D fluid models. The results for $C\left(t\right)$ with different $\tau$ values are presented in Figure 15. We can see from Figure 15 (a) and (b) that in 2D fluid models, under the condition of increasing the aspect ratio of the system, for heat conduction dominated by the hydrodynamic effect ($\tau =0.1$) or dominated by the kinetic effect ($\tau =10$), $C\left(t\right)$ will eventually change to a power-law decay $C\left(t\right)\sim t^{-2/3}$. As is shown in Figure 15 (c) and (d), in 3D fluid models, there is also a similar phenomenon that $C\left(t\right)$ will eventually change to a power-law divergence for the hydrodynamic effect ($\tau =0.1$) and the kinetic effect ($\tau =1$).

7. Heat Transfer with Magnetic Field

Another issue that can be studied through the MPC dynamics concerns the influence of a magnetic field on transport [47,48,49,50]. It is generally believed that heat conduction behavior is normal in low-dimensional systems where momentum is not conserved [10,11]. However, there is a counterexample in a low-dimensional system with magnetic field. Specifically, heat transport via the one-dimensional charged particle systems with transverse motions is studied in [47], where researchers studied two cases: case (I) with uniform charge and case (II) with alternate charge. An intriguing finding of this study is that in both cases involving non-zero magnetic fields, the heat conduction behaviors exhibit anomalies, similar to the case where momentum is conserved under the zero magnetic field condition. Remarkably, the abnormal behavior in case (I) is different from the case without magnetic field, suggesting a novel dynamical universality class. Due to the presence of the magnetic field, the standard momentum conservation in such a system is no longer satisfied but is replaced by the pseudomomentum conservation [51]. Thus, there are two relevant questions: (1) Does the pseudomomentum conservation of a system lead to abnormal heat conduction? (2) Can the abnormal behaviors in both cases also be observed in low-dimensional fluids under the same pseudomomentum conservation?

The above two questions have been well answered only recently in our research [52], where it is shown that under the same pseudomomentum conservation, the 2D fluid system with magnetic field can exhibit normal heat conduction behavior. Specifically, we consider a 2D system of charged particles as depicted in Figure 1 (b). In this system, a constant magnetic field perpendicular to the plane of motion, $\mathbf{B}=B\mathbf{k}$, is imposed. The particles interact via the modified MPC dynamics to maintain the pseudomomentum conservation of the system (see [52] for details). To compare with the results obtained in [47], we also consider two cases: case (I) with uniform charges $e_i=1$ and case (II) with opposite charges on each half of particles, say $e_i={(-1)}^i$. In Figure 16, we plot the relation of $\kappa$ vs L for various B obtained by nonequilibrium thermal-wall method. It is shown that for $B=0$, the system with momentum conservation exhibits the crossover from 2D to 1D behavior of the thermal conductivity under the condition of increasing the aspect ratio of the system. However, for $B\ne 0$, heat conduction behaviors in both cases with pseudomomentum conservation are normal because as L increases, $\kappa$ approaches a finite value.

Obviously, our above results are at variance with the findings in [47]. There, heat conduction in presence of pseudomomentum conservation in two cases are abnormal. This observation, together with our results, thus clarify that pseudomomentum conservation is not related to the normal and anomalous behaviors of heat conduction and provide an example of the difference in heat conduction between fluids and lattices in the presence of the magnetic field condition.

For completeness, we also mention that the MPC scheme can be extended to the case of charged particles that yields a self-consistent electric field. This situation is relevant for plasma physics and can be treated by coupling the MPC dynamics with a Poisson solver (see [53] and references therein for details). The effect of the electric field on heat transport can be studied by this method: simulations reveal that the field does not affect significantly the hydrodynamics of the model, at least for not too large amplitude fluctuations [53].

8. Conclusions

We have presented a series of numerical simulation demonstrating how the MPC method can be effectively employed to study the dimensionality effects on heat transfer in a simple confined fluid. Nonequilibrium dynamics can be simulated efficiently with the thermal-wall modeling of external reservoir and the results agree very well with Green-Kubo linear response. The data are statistically very accurate and span over a considerable range of system sizes. The overall theoretical scenario is confirmed by the data. It should be noticed that most of the publications in this context refer to lattice systems [11], so our results represent are a relevant extension to the case where particle are free to diffuse through the simulation box. This supports the general validity of low-dimensional hydrodynamic theories and of the tight connection with Kardar-Parisi-Zhang physics in transport problems [44].

Another relevant finding is that the crossover from diffusive to anomalous regimes, seen in quasi-integrable chains [39], extends to the somehow simpler case of fluids. In particular the decomposition of the current Eq.(9) is an effective and simple way to assess the divergence law when interaction is relatively weak and the accessible range of sizes too limited (as frequently occurs in practice).

We also illustrated how the important issue of dimensional crossovers and the effect of an applied magnetic field can be studied relatively easily via the MPC dynamics. A further extension would be to introduce the effect of chemical baths, namely to account for the exchange of particles with environment [6,54]. This would allow to study basic features of coupled transport process but also design and conceive novel possible applications.