Kinetic and hydrodynamic regimes in multi-particle-collision dynamics of a one-dimensional fluid with thermal walls

We study the non-equilibrium steady-states of a one-dimensional ($1D1V$) fluid in a finite space region of length $L$. Particles interact among themselves by multi-particle collisions and are in contact with two thermal-wall heat reservoirs, located at the boundaries of the region. After an initial ballistic regime, we find a crossover from a normal (kinetic) transport regime to an anomalous (hydrodynamic) one, above a characteristic size $L_*$. We argue that $L_*$ is proportional to the cube of the collision time among particles. Motivated by the physics of emissive divertors in fusion plasma, we consider the more general case of thermal walls injecting particles with given average (non-thermal) velocity. For fast and relatively cold particles, short systems fail to establish local equilibrium and display non-Maxwellian distributions of velocities.


I. INTRODUCTION
Statistical systems driven away from equilibrium by external agents, like concentration and/or thermal gradients, forced flows etc. are ubiquitous in nature and in technological applications. Especially in far-fromequilibrium regimes simulation of model systems is a basic tool to gain insight on the fundamental properties of the system under study. Molecular dynamics is the most natural approach, but methods based on effective stochastic processes may be a computationally convenient alternative. In particular, the Multi-Particle Collision (MPC) dynamics can be considered as a mesoscale simulation method where particles undergo stochastic collisions. 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 MPC dynamics is a useful tool to investigate concrete systems like polymers in solution, colloidal fluids etc. but also to address fundamental problems in statistical physics and in particular the effect of external sources.
In the present work, we investigate the transport property of a one-dimensional gas driven off-equilibrium by different reservoirs, modeled as thermal walls. This general setup has been considered very often in the literature [3,4] with several possible applications [5]. For the discussion that follows, it is worth recalling that the prediction of Nonlinear Fluctuating Hydrodynamics [6,7] for 1D systems, constrained by three global conservation laws, is that transport is generically anomalous and belongs to the Kardar-Parisi-Zhang (KPZ) dynamical universality class. This is confirmed by other approaches including hydroynamics [8,9], kinetic theory [10][11][12] and many numerical simulations [13][14][15]. However, for finite systems some deviations are found and it is thus of interest to clarify which are the different transport regimes (see e.g. Refs. [16,17] for an up-to-date account). In this framework, simple kinetic model are an unvaluable testbed.
Besides the general theoretical motivations, our interest in also in application of the technique in plasma physics [18,19]. We wish to investigate the MPC method as a tool to study nonequilibrium properties of fusion plasma which are subject to large temperature gradients. This is a relevant issue in the modeling of plasma transport in the direction parallel to the magnetic field in the edge of magnetic fusion devices where hot plasma regions are connected to a wall component. In this case a significant temperature gradient along the field line sets in between the hot region, that acts as a heat source, and the colder plasma region at the wall, that acts as a sink. Deviations from the classical Fourier law have been observed (see for example [20,21] and references therein) together with transitions from strongly collisional to collisionless regimes, exhibiting different heat conduction features. We want to point out that accounting for these kinetic effects by a fluid model of heat transport is of primary importance for implementing realistic and efficient hydrodynamic simulations. In this perspective, what contained in this work can be considered as a worthwhile evolution along the line of simplified 1D fluid models, already adopted for studying thermal transport in plasmas [22][23][24]. For the sake of brevity, we adopt for such model the label 1D1V used sometimes in the literature, meaning that we consider particles constrained on the line with one velocity component along it.
Besides this, in the edge plasma of magnetic fusion de-vices several phenomena, like for example neutral beam injection are related with the injection of hot and/or cold plasma particles. We refer for example to hot sources related to external heating of the plasma but also to particular regimes with emissive divertor targets on the wall due to the high temperatures associated with steady state loads, as expected for example on next step large fusion devices [25,26]. Here, we will thus attempt to model this situation, by generalizing the usual thermal wall method to include the possibility of injecting particles with an average non-thermal velocity. To our knowledge this case has not been treated in the statistical physics literature and it has thus an interest by itself. The outline is as follows. In Section II we recall the definition of the fluid dynamics and the thermal walls. Section III discusses the crossover from the so called kinetic (diffusive) regime to the anomalous, hydrodynamic one. The effect of general thermal walls is illustrated in Sec. IV with reference to the specific case of shifted Maxwell-Boltzmann distributions.

II. 1D1V MPC WITH THERMAL WALLS
For a detailed description of the MPC simulation scheme we address the reader to references [18,19]. Here, we just summarize the basic ingredients. We consider an ensemble of N point particles with equal masses m located in a finite space region, [0, L], partitioned into N c (fixed) cells of size a, L = N c a. The density of particles, d = N/L, is fixed, while all physical quantities, without prejudice of generality, are set to dimensionless units, namely a = 1, m = 1 and (for the isolated system) E c = 1, the latter denoting the total kinetic energy of the fluid of particles.
The MPC collision is performed at regular times steps, separated by a constant time interval of duration δt: the velocity of the j-th particle in the i-th cell is changed according to the update rule v j,old → v j,new = a i w j + b i , where w j is randomly sampled by a thermal distribution at the cell temperature T i , while a i and b i are celldependent parameters, determined by the condition of total momentum and total energy conservation in the cell [18]. After the collision all particles freely propagate, i.e.
(1) and they may move to a new cell. The non-equilibrium state is imposed by two thermal walls [27] acting at the edges of the space region: any particle crossing the boundaries of the interval [0, L] is reinjected in the opposite direction with a random velocity extracted form the Maxwell-Boltzmann distributions Re-injected particles propagate with their new velocities for the remaining time up to the next collision. Since arbitrarily large values of v are admitted, a particle may reach the opposite boundary during this time, although such a process, in practice, becomes extremely unlikely, i.e. negligible, for sufficiently large values of L, considering that in the adopted dimensionless units the typical value of v is O(1).
The heat fluxes at the two boundaries, Q 0 and Q L , respectively, are computed as the average kinetic energy exchanged by particles interacting with the walls [3, 4]. The relaxation of the fluid to a steady state, i.e. Q 0 ≈ −Q L = Q, is monitored by controlling the convergence of the running-time averages of the fluxes. The thermal conductivity is defined as K = Q/(∆T /L), where ∆T = T 0 − T L . Typically, in order to speed up the relaxation to a steady state, the initial velocities of the particles in cell i are randomly extracted from a Maxwell-Boltzmann distribution at temperature according to the linear temperature profile compatible with the Fourier's law.
The relevant kinetic parameters of the MPC protocol are the particle mean free path ℓ, which, in an uniform system, is proportional to the thermal velocity of the fluid v T = √ T , in formulae ℓ ∼ v T δt , and the thermal conductivity, that is expected to be proportional to the typical kinetic energy of the fluid, in formulae K = Cv T ℓ ∼ v 2 T δt, C being the specific heat at constant volume.
We conclude this section by making the reader aware of the main advantage of MPC simulations, with respect to those typically performed in nonlinear lattice models. The possibility of observing crossover between different transport regimes in the latter class of models depends on the fine tuning of various simulation parameters, e.g. energy density, integration time step, coupling with the thermal reservoirs, and all that are strongly model-dependent, because the crossover effects are ruled by the lifetime of typical nonlinear excitations, influencing the stationary transport process. As we show in what follows, in MPC simulations one can control different regimes essentially by a single natural parameter, the collision time δt, which is straightforwardly related to the mean-free path of the fluid particles, weak and strong collisional regimes corresponding to large and small values of δt, respectively.

III. FROM KINETIC TO HYDRODYNAMIC REGIMES
The kind of stationary heat transport phenomenon one is dealing with is characterized by the scaling of Q with the system size L. Diffusive transport, i.e. Fourier's law, yields Q ∝ L −1 . For the 1D1V MPC, based on general theoretical arguments [6, 7] we expect that transport would be anomalous and belonging to the KPZ universality class, as indicated by equilibrium correlation functions [18]. In the nonequilibrium setup used in the present work this would imply that However, this asymptotic anomalous regime may be eventually attained going through a crossover from a standard kinetic, i.e. diffusive, regime [17,28,29]. This scenario can be rationalized by assuming that the flux is made of two contributions where Q N is the "normal" part of the flux that can be written as (with q and r being suitable parameters) while Q A is the anomalous part that, for large enough values of L, should scale as and b ∼ O(1) is a characteristic length [17]. Formula (5) is suggested by kinetic theory to account for the crossover from initial ballistic regime to a diffusive one, see e.g.
[3]. It implies that up to system lengths of the order of the mean free path L ∼ ℓ transport is essentially ballistic, with a flux independent of L while for L ≫ ℓ, Q N ∼ ℓ L . According to the above formulae, for large enough L one can define a crossover length L * from the normal to the anomalous regime when Q N ≈ Q A , i.e.
Altogether, upon increasing L at fixed ℓ, one should see a first ballistic regime followed by a kinetic (diffusive) one until eventually the asymptotic hydrodynamic regime is attained. It should be however remarked, that the different regimes are observable only provided that the relevant length scale are widely separated, in such a way that the range of scales between the ballistic and anomalous regimes is large enough i.e. ℓ ≪ L * . The numerical results reported in Fig.1(a) provide a clear confirmation of this scenario. The energy flux Q(L) displays a crossover from diffusive to anomalous scaling upon decreasing the collision time δt. In Fig.1(b) we also show that, for low collisionality and small L, the thermal conductivity K, obtained for different values of δt, tend to collapse onto a single curve after rescaling K and L by δt, or, equivalently, by ℓ . This confirms that in the kinetic regime ℓ is the only relevant scale as expected.
A more refined analysis of the data, supporting Ansatz (4), can be obtained by estimating the two contributions in Eq.(4). This is accomplished by fitting rescaled conductivity data with the largest δt via the functional form AL B+L suggested by Eq. (5), with A, B being fitting parameters.
As seen in Fig.1(b) (dash-dotted line) the fitting is very accurate meaning that the anomalous contribution to the conductivity is negligible over the considered length range. From the fitting function, one can obtain straightforwardly the best estimate of Q N and thus the anomalous part as Q A = Q − Q N . As shown in Fig.2, a power law fit of Q A indicates an excellent agreement with the expected scaling Eq.
(3) on a remarkably large range of about three decades in L. To our knowledge, this one of the most solid numerical confirmations of the predictions of fluctuating hydrodynamics. We also note, that the characteristic length b is found to be independent of the collision parameter. The difference between the two regimes can be further appreciated also in the temperature and density profiles computed as averages of v 2 i and of the number of particles in each one of the N c cells. In order to deal properly with the large L limit, we plot the temperature and the density profiles, T (x) and n(x), respectively, as a function of rescaled variable x = i/L, where i is the cell index. This representation allows one to conclude also that local equilibrium is attained, since the product T (x)n(x) is x-independent, as expected for an ideal gas (data not shown).
In Fig.3 we report these profiles for two different values of the collision time and different sizes L. For L < L * the temperature profiles are in good agreement with the As shown in Fig.3(a) this functional form accounts for the measured shapes for large enough L, apart the temperature jumps at the edges, that are a typical manifestation of boundary impedance effects. A similar situation was found also for the HPG model [31,32]. Also, considerations of the same type apply to the density profiles in Fig.3(b), in view of the relation T (x)n(x) ≈ constant.
In the hydrodynamic regime, i.e. L > L * , the temperature profiles exhibit instead the typical signatures of anomalous transport, signaled by singularities at the boundaries (see Fig.3(c,d)). For instance, the inset in Fig.3(c) shows that T (x) ∼ (1 − x) µ for x → 1. This property is traced back to the fact that steady state profiles should satisfy a fractional heat equation [33]. The parameter µ has been termed meniscus exponent µ, as it described the characteristic curvature close to the edges. In the simulations here µ ≈ 3/4, a value already found for other models with reflecting boundary conditions [33,34].

IV. GENERAL THERMAL WALLS
Having in mind the possibility of simulating the injection of hot and/or cold particles into a fusion plasma by the MPC protocol, in this Section we describe how one can cope with this task by introducing suitable and more general boundary conditions. Within the adopted thermal wall scheme this can be achieved by changing the distributions of re-injected particles at one boundary of the space region. For instance, it can be assumed that the distribution p 0 on the left boundary is replaced by a function of the general form [35] p 0 = vH(v) .
By imposing the condition ∞ 0 vH(v)dv = 1, we ensure that the net flux of particles vanishes [35]. This distribution can be characterized by two parameters, namely the typical average velocity of the injected particles that measures the source spread [36]. As a specific example, here we consider a shifted Maxwellian distribution parametrized by the velocity V , where H 0 is a suitable normalization factor. Shifted Maxwellian sources are implemented in kinetic models to account for drift speed of plasmas [37]. Manifestly, both v 0 and ∆v are functions of V and T 0 . In the following simulations, we generate random velocities drawn from (9) and use V and T 0 as control parameters, reporting also the corresponding values of v 0 and ∆v for comparison [38].
In Fig.4 we present the results of MPC simulations for different values of δt and V . We first of all note that, even setting equal temperatures at the thermal walls, i.e. T 0 = T L , the system is out of equilibrium and there is a net flux of energy and particles. For what concerns the scaling with of the flux with the size L, the effect of these thermostats is pretty similar to the standard case discussed in the previous Section. As shown in Fig.4a, in the case of relatively large collisionality the anomalous scaling, Eq.(3) is attained. For weak collisionality, Fig.4b the diffusive scaling is found in the considered range. Although not observed in the data, we expect the same crossover described above beyond the hydrodynamic scale. In both cases, as expected, one observes that the magnitude of the current increases upon increasing the speed of the injected particles. The above results are not obvious, and indicate that regardless of the nature of the baths the overall scenario is mostly dictated by bulk properties.
To better appreciate the effect of this kind of thermal wall we measured the position-dependent steady state distributions f (x, v) in phase space for a moderate collisional parameter δt = 10, see Fig.5. We consider two cases: 1. the typical entry speed v 0 is of the same order of the thermal velocity v 0 ∼ ∆v = √ T 0 (left panels in Fig.5). In this case the situation is very similar to the standard thermal bath: temperature and density gradients are formed, while velocity distributions are Maxwellian even relatively close to the left boundary, thus indicating that local equilibrium has been achieved pretty fast.
2. fast, i.e. "cold", particles are injected, i.e. the typical entry speed v 0 is larger than ∆v (right panels in Fig.5). In this case one observes that the beam of particles entering from the left wall propagates, while experiencing the effect of the interactions with the cloud of almost-thermalized particles only in the center of the space region. As a result, the velocity distributions are strongly non-Maxwellian and remain double-humped in the first half of the space region.
Adopting a kinetic point of view, we expect that the typical size of the region in which the fluid is not in local equilibrium would be of the order of mean free path ℓ. For the simulations above, we checked that upon increasing L such a region remains indeed finite. On the other hand, in anomalously conducting systems lack of thermal equilibrium close to the sources is found [39] and we cannot exclude that the same occurs here for larger systems.
To conclude this Section, we briefly discuss the dependence of fluxes and profiles on the thermal wall parameters. In the top panel of Fig. 6 we plot the contour levels of the energy flux Q as a function of V and T 0 keeping the other parameters fixed. As remarked above, Q increases with V as intuitively expected. However there is an interesting feature occurring for small enough T 0 and V . For T 0 < T L and V small, injected particles are still relatively slow and thermalise quickly so energy mostly flows from right to left, Q < 0. The corresponding distributions and profiles T (x) and n(x) are shown in the leftmost panels of Fig. 6, where it is seen that T is decreasing. However, increasing the speed of the injected particles leads to a sign reversal of the current, Q > 0 and of the associated gradients (right panels of Fig. 6). A rough estimate of the conditions where this inversion takes place, can be obtained by equating the kinetic temperature T L of the rightmost reservoir with (twice) of the kinetic energies of the injected particle, which roughly gives T L ≈ T 0 + V 2 . As seen in Fig. 6(a) this is in agreement with the numerical data. The above estimate (dashed white line) indeed separate the two regions of the parameters plane having different signs of the energy current. On the basis of this observation, one may argue that the velocity V can be used as a further control parameter of the flux direction.

V. CONCLUSIONS
In this work we have considered properties of nonequilibrium steady-states and transport in a one- For Maxwellian walls we demonstrated a crossover from a normal (kinetic) transport regime to an anomalous (hydrodynamic) one, above a characteristic size L * . Since L * grows cubically with the collision time δt, the anomalous regime may hardly be detected with the sizes typically used in simulations. So we argue that a standard diffusive description will account for the observation for short enough systems. Altogether, the results nicely fit in the general framework observed for other models [17,28,29]. They call for a critical consideration of the usual assumptions made in kinetic theory, since memory effects and long-time tail may play a major role in predicting heat and particle transport in low-dimensions.
A natural question is whether the same scenario occurs in two dimensional fluids. Here, the expectation is that in the hydrodynamic regime anomalous transport would be associated with a logarithmic divergence of the conductivity with size. For the 2D MPC dynamics this is confirmed by equilibrium simulations [19]. Thus in principle a similar crossover scenario should occur, although it may be hard to be observed in view of such a weak logarithmic dependence.
In the second part, we considered the more general case of thermal walls injecting particles with given av-erage (non-thermal) velocity. We showed, that the kinetic and hydrodynamic regimes are observed confirming that transport is dominated by bulk correlations in both cases. For injection of fast and relatively cold particles, short systems fail to establish local equilibrium and display non-Maxwellian distributions of velocities. We also showed that the injection velocity may be used to control current reversal.
To conclude, let us comment on the perspective of applying the model to confined plasma. This requires considering the effect of the self-consistent electrostatic field, solution of the associated Poisson equation. As already noticed in Refs. [41,42], the case of equal masses and charges in one dimension can be mapped in the socalled ding-dong model [43]. The latter consists of a one dimensional array of identical harmonic oscillators (each one with frequency ω p , the plasma frequency) undergoing hard-core collisions [43]. Non-equilibrium and transport properties of such a model have been investigated in detail [41][42][43] and normal diffusive transport has been demonstrated, as expected since energy remains the only conserved quantity. One may thus argue that adding MPC dynamics will not affect significantly such a property and that the hydrodynamic regime will not be present. However, for finite systems the situation may be more subtle. Actually, the field introduces another relevant time scale in the dynamics, i.e. the inverse of the plasma frequency ω p , associated to the collective oscillations and Langmuir waves. This time scale should be compared with the collision frequency 1/δt. Generally speaking, if ω p δt ≪ 1 energy transport will be driven mostly by collisions, and it is not a priori excluded that hydrodynamic effects may play a role over a relevant range of scales. Those issues deserve investigation and will be considered in the future.

ACKNOWLEDGMENTS
This work has been carried out within the framework of the EUROfusion consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. This work is part of the Eurofusion Enabling Research project ENR-MFE19.CEA-06 Emissive divertor ; SL and RL acknowledge partial support from project MIUR-PRIN2017 Coarse-grained description for non-equilibrium systems and transport phenomena (CO-NEST) n. 201798CZL.