A lattice Boltzmann study on Brownian diffusion and friction of a particle in a confined multicomponent fluid

We study the diffusivity of a small particle immersed in a square box filled with a non-ideal multicomponent fluid in the presence of thermal fluctuations. Our approach is based on the numerical integration of fluctuating lattice Boltzmann models (LBM) for multicomponent mixtures. At changing the wettability on the particle's surface, we measure the mean square displacement (MSD) and compare with the prediction of the Stokes-Einstein theory. Two main set-ups are tested, involving periodic boundary conditions and wall boundary conditions realized on the computational box. We find that full periodic boundary conditions give rise to random advection after millions of lattice Boltzmann time steps, while this effect is mitigated in the presence of wall boundary conditions. The matching with the Stokes-Einstein relation is therefore guaranteed when we use the appropriate frictional properties measured in the presence of confinement. Our results will help the exploration of nanoscale applications with multicomponent fluids using LBM in the presence of thermal fluctuations. (This paper has been published in Journal of Computational Science)

Through the Chapman-Enskog expansion, the LBM can recover the hydrodynamic representation of the Navier-Stokes equations [18,19]. Recent studies have shown interests in integrating finite-size particles into the LBM framework, such as particle suspensions [20,21,22,23,24], pickering emulsions [25], self assembly particles [26,27,28]. However, there are relatively few LBM studies in the literature on particle-fluid interaction problems at the nanoscale, where the effects of the thermal fluctuations cannot be neglected. The deterministic hydrodynamics of Navier-Stokes equations is inadequate to describe the flow evolution at the nanoscales [29,30]. For this reason, there have been pioneering works to extend the LBM in the presence of thermal fluctuations, designing the so-called fluctuating lattice Boltzmann models (FLBM) [31,32,33,34,35,36]. Recently, FLBM have been used to study non-ideal multicomponent fluids [37,38] and also the effects of thermally excited capillary waves on the break-up properties of a thin liquid ligament [39]. Moreover, FLBM have been used to understand the influence of thermal fluctuations on the particle settling under confinement [40].
When the particle is coupled with the FLBM, it will naturally perform Brownian motion. In such conditions, the mean squared displacement (MSD) of the particle in an unconfined domain is expected to follow the Einstein's relation in time t [41,42] where the particle is located at x and the diffusion coefficient in unconfined domains is indicated with D unconf . Based on the fluctuation-dissipation relation [41], the diffusion coefficient D unconf can be represented by k B T/(6πηR p ), where k B T is the system thermal energy, η is the dynamic viscosity and R p is the particle radius. Notice that the denominator in D unconf is essentially the friction coefficient experienced by a spherical particle with radius R p in unconfined domains. Due to the large computational cost associated with the study of Brownian motion, only short-time Brownian motion studies of colloids have been proposed in the FLBM framework [2]. The long-time Brownian motion has rarely been studied. In this paper, we aim to quantitatively tackle the problem of the Brownian motion of a wetted particle immersed in a multicomponent fluid. In Section 2, we introduce our multicomponent FLBM algorithm which couples with finite-size wetted particles. In Section 3, we present our simulation setup. Results are discussed in Section 4, where we show that the particle experiences a random advection by using periodic boundary condition which does not allow a fair assessment of the Einstein's relation with the diffusion coefficient predicted by the fluctuation-dissipation relation [41]. We show how to remove this pathology when working with wall boundary conditions; the precise matching with the Einstein's relation, however, requires the knowledge of the friction in the presence of confinement. We draw our conclusions in Section 5

METHODOLOGY
In this section, we introduce the multicomponent FLBM for the binary mixture of two fluids [37]. Also, the wetted finite-size particle model [20,21,25], which is used to simulate the solid particle, is coupled with the fluctuating multicomponent fluids. Technical details of the FLBM have already been extensively presented in [37], and here we only briefly recall the most important aspects for the sake of completeness. We consider f li (x, t) as the discretized particle's probability distribution function on the i-th direction of a lattice cell with velocity teractions, and stochastic noise, can be written as: where Ω is a collision kernel [3,4], F li is forcing term for fluid-fluid interactions, and ξ li is a stochastic source. The collision kernel is chosen to be the multiple relaxation time (MRT) collision kernel [43,34,44]. The basic idea of the MRT kernel is to introduce a moment space to decompose the probability functions into "modes". The lower-order modes are related to hydrodynamic quantities (density, momentum and stress tensor). The higher-order modes are the "ghost modes", which do not contribute to the hydrodynamic behavior [43,34]. The collision kernel relaxes the distribution function towards the local equilibrium f eq li , which is given by in which ω i is a suitable weight needed to impose the isotropy in the interaction, ρ l (x, t) and u (x, t) are the hydrodynamic macroscopic density for each component and mixture velocity, respectively, which can be calculated from the distribution function: where ρtot(x, t) = l ρ l (x, t) is the total density of the two components. The term F li is the Shan-Chen forcing term which models the non-ideal interactions for the mixture [45,46,47,48,49]. The Shan-Chen forcing term is given by where G is a strength coefficient, ϕ l is known as the pseudo-potential function which is set to the fluid density ϕ l = ρ l . In all the simulations performed, the coupling coefficient is set to G = 1.5 lattice Boltzmann units (lbu). The term (2) is a stochastic force. The noise term does not influence on the conserved density modes, while higher modes receive the stochastic source following the fluctuation-dissipation relation [37]. Through the Chapman-Enskog expansion [34,44], we can recover the equations for the fluid densities and the hydrodynamical velocity [30] (superscript T means transposition): where the hydrodynamical velocity u (H) is given by u (H) = u+(F A +F B )/2ρtot, P b and µ are the bulk pressure and the chemical potential, respectively [49].
The mass diffusion coefficient M and η are related to the relaxation times of the fluid. In all the simulations performed, the mass diffusion coefficient is set to M = 1/6 lbu, and the dynamic viscosity is set to η = 0.383 lbu. We use the same kinematic viscosity in both fluids (i.e. the same relaxation time in the LBM framework). The stochastic stress (Σt) and the stochastic diffusion (Ψv) contribution to the equations of hydrodynamics are given by where k B T is the thermal energy, while and Wt and Wv are a Gaussian tensor and a Gaussian vector with independent and uncorrelated components and variance equal to unity. For the particle model, the technical details have been introduced in previous studies, readers can follow references [20,21,25]. Here we recall the basic ingredients. We integrate the wetted finite-size particle in the FLBM framework. The particle is built on lattices by declaring nodes belonging to the particle ("particle nodes"). The particle exchanges the momentum with surrounding fluids through the bounce back [20]. The particle moves due to the When w > 0 (hydrophobic case), the particle's contact angle is larger than 90 • ; when w < 0 (hydrophilic case), the particle's contact angle is smaller than 90 • . Black and red colors refer to two different resolutions for the particle radius (Rp = 10, 4 lbu). Small particle radii are enough to retrieve converged results.
action of the external forces on it; due to this motion, new lattice nodes that were originally belonging to the fluid will belong to the particle (cover-node behavior). Consistently, new fluid nodes are generated which were initially belonging to the particle (uncover-node behavior). To impose the total mass conservation, the mass correction algorithm described in Ref. [25] has been implemented in this paper. We introduce the virtual fluid in the outer-most layer of the particle for tuning the particle's wettability. The densities in the virtual fluid nodes for the two components are defined as ρ A,v and ρ B,v . The parameter w is a dimensionless number that regulates the affinity of the particle towards one of the two fluids [25]. The wettability properties described in the following (i.e. hydrophobic, neutral, hydrophilic) refer to the affinity of the particle towards the majority component in the bulk phase. A positive (negative) w will corre-spond to a hydrophobic (hydrophilic) case. We also placed the wetted finite-size particle at the fluid-fluid interface and measured the contact angle of the particle by tuning the parameter w. Fig. 1 shows the contact angle as function of w. Two different radii of the particle have been tested: R p = 4 lbu (red line) and R p = 10lbu (blue line). The results show that the particle wettability can be tuned from around 30 • contact angle to 150 • contact angle. The solvent solver is based on the implementation of a multicomponent fluctuating lattice Boltzmann model (FLBM) [37]. The wetted particle is coupled with the binary fluid [21,25].

NUMERICAL SETUP
The numerical set-up for the Brownian motion simulation is shown in Fig. 2.
We place a particle with radius R  Fig. 3 shows the concentration of the two components in proximity of the particle's surface and in the bulk. The density of the particle is set to ρ p = 2ρtot. We define the intrinsic diffusion timescale t diff as the time it takes for fluid perturbations to diffuse over a lengthscale comparable to a particle radius [2] t diff = R p 2 /ν.
To perform simulations and match the asymptotic Brownian motion behavior    with the Einstein's prediction, we need at least 10 4 time units. Thus, we cannot use very high resolution for the particle and we chose the radius of the particle equal to R p = 4 lbu. The domain size is set to 60 3 lbu, and the timescale t diff is t diff = 96 lbu. The wetting parameters of the particle are set to w = −0.88 (hydrophilic), w = 0 (neutral), and w = 0.32 (hydrophobic). The thermal energy is set to k B T = 1 · 10 −4 lbu. For each wetting parameter, we perform simulations with 10 7 lbu timesteps with 10 realizations of the noise to get enough statistics.

RESULTS
We start by analyzing a case with periodic boundary conditions [2]: the particle is initially placed in the center of the cubic box and is kicked by the fluctuating solvent. Fig. 4 presents the MSD as a function of the dimensionless time t/t diff at changing wetting parameters w. For t/t diff < 100, we observe a neat convergence towards the expected theoretical result. However, the convergent behavior is lost when t/t diff > 100, and the MSD eventually deviates from the Einstein's relation. The reason for the mismatch is traced back to the fact that the particle is modeled on the lattice, and its movement gives rise to a "cover-uncover" behavior on the lattice itself [2,25]. The cover-uncover behaviour, by definition, introduces spurious momentum contributions, because when nodes are covered, some others are uncovered, and there is inevitably an adaptation dynamics that is needed for the new uncovered nodes to "equilibrate" with the surrounding fluid. These spurious momentum contributions should be small for large particles; unfortunately, particles need to be sufficiently small to make the Brownian diffusion simulation affordable in terms of computational time. Consequently, after some time, the system will develop a spurious advection flow in a random direction. The random advection is removed if we design a new set-up, where we place walls with a neutral wetting boundary condition at the boundaries. Results of numerical simulations with this new set-up are reported in Fig. 6. Results are normalized by the expected asymptotic result for an unconfined domain. We observe, however, that the normalized MSD for the three wetting conditions does not approach a unitary value for large times.
To delve deeper into the observed discrepancy, we notice that the particle can experience confinement effects [40]. Hence, the theoretical estimate of the diffusion coefficient should be performed by including the friction in the presence of the bounding walls. The diffusion coefficients for the confined and unconfined domains, which are D conf and D unconf respectively, are defined as conf , where F p is the external forcing acting on the particle, and U (z) conf is the drift velocity which can be obtained from a dedicated simulation. We thus performed particle settling simulations under confinement without thermal fluctuations [50] for hydrophilic, neutral and hydrophobic cases. The domain size is the same as the Brownian motion simulation which is 60 3 lbu. We place the particle in the center of the domain, and apply the gravity acceleration g = 10 −6 lbu in the z direction. After a transient time, the particle's velocity reaches a stationary state, which is defining the drift velocity. We then define c m as the ratio between the particle's drift velocity under confinement and the Stokes' prediction for an unconfined domain driven by the same volume force: Figure 5 shows the ratio c m as function of the wetting parameter w. Three wetting parameters, which are w = −0.88, w = 0,w = 0.32 representing hydrophilic, neutral and hydrophobic cases respectively, have been investigated.
As we can see, the ratio c m is always smaller than 1, indicating that the friction under confinement is larger than the Stokes's friction. Also, we observe that a larger errorbar appears in the hydrophobic case, which indicates that the particle is more sensitive to external fluid kicks in this situation. According to Eq. (11) and Eq. (12), the diffusion coefficient D conf can be obtained as: Consequently, the theoretical expectation for the diffusion coefficient in the presence of confinement is smaller than the unconfined result. In Fig. 7 we report the normalized MSD as function of t/t diff at changing wetting conditions. At difference with respect to Fig. 6, we now normalize the MSD with the diffusion prediction based on D conf . Now, the simulation results approach a unitary value for the three wetting conditions. Summarizing, the presence of the bounding box with walls allows to remove spurious random advection; however, this introduces small -but measurable -confinement effects. These effects need to be taken into account for a quantitative matching with the diffusion prediction at large times.

CONCLUSIONS
In this paper, we simulated Brownian diffusion of a wetted particle with radius R p in a cubic box L × L × L filled with the fluctuating binary mixture.
The solvent is simulated with the multicomponent FLBM [37], and different wettabilities of the particle have been considered based on [25]. To quantitatively understand the Brownian motion of the finite-size particle, we have tested two boundary conditions on the cubic box: periodic and neutrally wetted walls.
For the periodic set-up, the particle's motion experienced a random advection flow which causes a deviation from the Einstein's relation, due to the "coveruncover" behavior triggered by the particle's motion [25]. The use of neutrally wetted walls removes such pathology; however, in the presence of confinement, the friction coefficient needs to be corrected to allow for a precise matching with the Einstein's relation. On a future perspective, one could consider the present contribution a precursor study for the more complicated situation where the wetted finite size particle rests at the interface separating 2 immiscible fluids [51]. Also, further investigations on the dynamics of Brownian particles in shear flow to study the Taylor dispersion could be interesting. Moreover, we would like to mention that the particle is still not very well resolved. Hence, it could be a challenging computational task to perform numerical simulations with larger resolutions to tackle these kind of problems.