Diffusion under confinement: hydrodynamic finite-size effects in simulation

TL;DR

This paper analyzes how finite-size effects influence diffusion measurements in confined fluids through molecular dynamics simulations and hydrodynamic theory, providing analytical corrections for simulation artifacts in narrow and wide nanopores.

Contribution

It derives analytical expressions to correct finite-size effects in diffusion simulations of confined fluids, applicable across various simulation methods and pore geometries.

Findings

Finite-size effects significantly impact diffusion coefficients in confined fluids.

Analytical corrections match simulation results well except in very narrow pores.

Wide nanopore diffusion coefficients are affected by previously unrecognized finite-size artifacts.

Abstract

We investigate finite-size effects on diffusion in confined fluids using molecular dynamics simulations and hydrodynamic calculations. Specifically, we consider a Lennard-Jones fluid in slit pores without slip at the interface and show that the use of periodic boundary conditions in the directions along the surfaces results in dramatic finite-size effects, in addition to that of the physically relevant confining length. As in the simulation of bulk fluids, these effects arise from spurious hydrodynamic interactions between periodic images and from the constraint of total momentum conservation. We derive analytical expressions for the correction to the diffusion coefficient in the limits of both elongated and flat systems, which are in excellent agreement with the molecular simulation results except for the narrowest pores, where the discreteness of the fluid particles starts to play a role. The present work implies that the diffusion coefficients for wide nanopores computed using elongated boxes suffer from finite-size artifacts which had not been previously appreciated. In addition, our analytical expression provides the correction to be applied to the simulation results for finite (possibly small) systems. It applies not only to molecular but also to all mesoscopic hydrodynamic simulations, including Lattice-Boltzmann, Multi-Particle Collision Dynamics or Dissipative Particle Dynamics, which are often used to investigate confined soft matter involving colloidal particles and polymers.

Equations

$$D_{\parallel }(H,\infty )$$

$$\Delta D_{\parallel }=D_{\parallel }(H,L)-D_{\parallel }(H,\infty ),$$

$$\Delta D_{\parallel }$$

$$⟨\delta {\mathbb{T}}^i⟩=\frac{1}{H}{\int }_0^H\mathrm{d}z\sum _{(m,n)\ne (0,0)}{\mathbb{T}}^{\infty }(-mL,-nL,z,z).$$

$$⟨\delta {\mathbb{T}}^b⟩=\frac{1}{H^2}{\iint }_0^H\mathrm{d}z\mathrm{d}{z}_0{\iint }_{-\infty }^{\infty }\frac{\mathrm{d}{x}_0\mathrm{d}{y}_0}{L^2}{\mathbb{T}}^{\infty }(-x_0,-y_0,z,z_0).$$

$$\Delta D_{\parallel }(H>L)$$

$$\Delta D_{\parallel }(H<L)$$

$$D_{\parallel }(z)$$

$$\frac{⟨D_{\parallel }⟩}{D_{\infty }}$$

$$=1-\frac{9\sigma }{8H}\sqrt{\frac{7\sigma +16H}{25\sigma +16H}}\mathrm{Argth}\left(\frac{16H}{\sqrt{(7\sigma +16H)(25\sigma +16H)}}\right)$$

$$\eta {\nabla }^2\mathbf{v}(\mathbf{r})+\delta (x)\delta (y)\delta (z-z_0)F_{\parallel }{\mathbf{e}}_{\mathbf{x}}-\nabla p(\mathbf{r})$$

$$\nabla \cdot \mathbf{v}(\mathbf{r})$$

$$\mathbf{v}(\mathbf{r})$$

$${\nabla }^{2}p(\mathbf{r})={\delta }^{\prime }(x)\delta (y)\delta (z-z_0)F_{\parallel }.$$

$$h(\mathbf{q},z)={\iint }_{-\infty }^{\infty }\mathrm{d}x\mathrm{d}yh(x,y,z)e^{i{q}_{x}x+i{q}_{y}y}$$

$$-q^{2}p(\mathbf{q},z)+\frac{{\partial }^{2}p(\mathbf{q},z)}{\partial z^2}$$

$$\eta q^2{\mathbf{v}}_{\parallel }(\mathbf{q},z)-\eta \frac{{{\partial }^2\mathbf{v}}_{\parallel }(\mathbf{q},z)}{\partial z^2}$$

$$\eta q^2{v}_z(\mathbf{q},z)-\eta \frac{{\partial }^2{v}_z(\mathbf{q},z)}{\partial z^2}$$

$$i\mathbf{q}\cdot {\mathbf{v}}_{\parallel }(\mathbf{q},z)+\frac{\partial v_z(\mathbf{q},z)}{\partial z}$$

$$\mathbf{v}(\mathbf{q},z)$$

$$\{\begin{array}{ll}p(\mathbf{q},z)& =A\cosh (qz)+B\sinh (qz)+\frac{i{q}_x}{2q}\exp (q|z-z_0|)F_{\parallel }\\ -\eta v_x(\mathbf{q},z)& =A^x\cosh (qz)+B^x\sinh (qz)-\frac{i{q}_x}{2q}zA\sinh (qz)-\frac{i{q}_x}{2q}zB\cosh (qz)\\ & +\frac{q_x^2}{4q^2}|z-z_0|\exp (q|z-z_0|)F_{\parallel }+\frac{1}{2q}\exp (q|z-z_0|)F_{\parallel }-\frac{q_x^2}{4q^3}\exp (q|z-z_0|)F_{\parallel }\\ -\eta v_y(\mathbf{q},z)& =A^y\cosh (qz)+B^y\sinh (qz)-\frac{i{q}_y}{2q}zA\sinh (qz)-\frac{i{q}_y}{2q}zB\cosh (qz)\\ & +\frac{q_x{q}_y}{4q^2}|z-z_0|\exp (q|z-z_0|)F_{\parallel }-\frac{q_x{q}_y}{4q^3}\exp (q|z-z_0|)F_{\parallel }\\ -\eta v_z(\mathbf{q},z)& =A^z\sinh (qz)+B^z\cosh (qz)-\frac{1}{2}Az\cosh (qz)-\frac{1}{2}Bz\cosh (qz)\\ & -\frac{i{q}_x(z-z_0)}{4q}\exp (q|z-z_0|)F_{\parallel }\end{array}$$

$$-\eta \frac{\partial v_z(\mathbf{q},z)}{\partial z}$$

$$-\frac{1}{2}Aqz\sinh (qz)-\frac{1}{2}Bqz\cosh (qz)$$

$$-\frac{i{q}_x}{4q}\exp (q|z-z_0|)F_{\parallel }-\frac{i{q}_x|z-z_0|}{4}\exp (q|z-z_0|)F_{\parallel }$$

$$\{\begin{array}{ll}\cosh (qz)(i{q}_x{A}^x+i{q}_y{A}^y+q{A}^z-\frac{1}{2}A)& =0\\ \sinh (qz)(i{q}_x{B}^x+i{q}_y{B}^y+q{B}^z-\frac{1}{2}B)& =0\end{array}$$

$$\{\begin{array}{l}A=2i\mathbf{Q}\cdot \mathbf{A}\\ B=2i\mathbf{Q}\cdot \mathbf{B}\end{array}$$

$${\bf A}=\left(\begin{array}[]{c}A^{x}\\ A^{y}\\ A^{z}\end{array}\right)\;;\;{\bf B}=\left(\begin{array}[]{c}B^{x}\\ B^{y}\\ B^{z}\end{array}\right)\;;\;{\bf Q}=\left(\begin{array}[]{c}q_{x}\\ q_{y}\\ -iq\end{array}\right)$$

$$\{\begin{array}{l}\begin{array}{ccc}-\eta v_x(\mathbf{q},z)& =& A^x\cosh (qz)+B^x\sinh (qz)-\frac{i{q}_x}{2q}zA\sinh (qz)-\frac{i{q}_x}{2q}zB\cosh (qz)-\eta v_x^0(z)\\ -\eta v_y(\mathbf{q},z)& =& A^y\cosh (qz)+B^y\sinh (qz)-\frac{i{q}_y}{2q}zA\sinh (qz)-\frac{i{q}_y}{2q}zB\cosh (qz)-\eta v_y^0(z)\\ -\eta v_z(\mathbf{q},z)& =& A^z\sinh (qz)+B^z\cosh (qz)-\frac{1}{2}Az\cosh (qz)-\frac{1}{2}Bz\cosh (qz)-\eta v_z^0(z)\end{array}\end{array}$$

$$\eta {\mathbf{v}}^0(z)$$

$${\mathbf{V}}^{(0)+}=\frac{1}{2}\left[{\mathbf{v}}^0(H/2)+{\mathbf{v}}^0(-H/2)\right]$$

$${\mathbf{V}}^{(0)-}=\frac{1}{2}\left[{\mathbf{v}}^0(H/2)-{\mathbf{v}}^0(-H/2)\right]$$

$$\{\begin{array}{l}\begin{array}{ccc}0& =& A^x\cosh (qH/2)-\frac{i{q}_x}{2q}\frac{H}{2}A\sinh (qH/2)-\eta V_x^{(0)+}\\ 0& =& B^x\sinh (qH/2)-\frac{i{q}_x}{2q}\frac{H}{2}B\cosh (qH/2)-\eta V_x^{(0)-}\\ 0& =& A^y\cosh (qH/2)-\frac{i{q}_y}{2q}\frac{H}{2}A\sinh (qH/2)-\eta V_y^{(0)+}\\ 0& =& B^y\sinh (qH/2)-\frac{i{q}_y}{2q}\frac{H}{2}B\cosh (qH/2)-\eta V_y^{(0)-}\\ 0& =& A^z\sinh (qH/2)-\frac{A}{2}\frac{H}{2}\cosh (qH/2)-\eta V_z^{(0)-}\\ 0& =& B^z\cosh (qH/2)-\frac{B}{2}\frac{H}{2}\sinh (qH/2)-\eta V_z^{(0)+}\end{array}\end{array}$$

$$\mathbf{A}$$

$$\mathbf{B}$$

$$A=\frac{2i\eta }{1-\frac{qH}{\sinh (qH)}}\left[\frac{1}{\cosh (qH/2)}\{q_x{V}_x^{(0)+}+q_y{V}_y^{(0)+}\}-\frac{iq}{\sinh (qH/2)}{V}_z^{(0)-}\right]$$

$$B=\frac{2i\eta }{1+\frac{qH}{\sinh (qH)}}\left[\frac{1}{\sinh (qH/2)}\{q_x{V}_x^{(0)-}+q_y{V}_y^{(0)-}\}-\frac{iq}{\cosh (qH/2)}{V}_z^{(0)+}\right]$$

Full text

The dynamics of fluids can be dramatically modified under confinement down to the molecular scale in nanotubes Striolo (2006); Falk et al. (2010) or nanopores Falk et al. (2015), due to the discreteness of matter and to the interfacial fluid-solid interactions. Even in larger pores, wider than tens of molecular sizes which are typical of nanofluidic devices Mathwig et al. (2012); Siria et al. (2013) and for which continuum hydrodynamics hold Bocquet and Charlaix (2010), confining walls influence the dynamics of the fluids and solutes Hagen et al. (1997); Huang and Szlufarska (2015). In particular, the diffusion coefficient of particles along a wall is generally reduced due to the friction at the interface.

Quantitatively, the solution of the Stokes equation in a slit pore involves a series of contributions of hydrodynamic images situated inside the solid walls. It predicts a decrease in the diffusion coefficient along the surface with respect to the bulk value, of leading order $\sigma/d$, with $\sigma$ the diameter of the particle and $d$ the distance to the surface Faxén (1922); Blake and Chwang (1974); Swan and Brady (2010). After averaging over the hydrodynamic slab width ${H}$, this results in a decrease governed by Sup :

[TABLE]

where $\infty$ refers to infinite lateral dimensions of the slit on the left-hand side and to the bulk fluid on the right-hand side. Using mode-coupling theory, Bocquet and Barrat obtained a similar scaling, in good agreement with molecular dynamics (MD) simulations, and emphasized its origin as the suppression of long-wavelength modes due to confinement Bocquet and Barrat (1995). Slippage at the interface changes the dependence of the diffusion coefficient with distance to the surface Lauga and Squires (2005) and may result in some cases in an average diffusion coefficient larger than the bulk value Saugey et al. (2005).

MD simulations have generally confirmed the decrease in the diffusion coefficient near a variety of model and realistic boundaries Liu et al. (2004); Marry et al. (2008); Sendner et al. (2009); Botan et al. (2011); Wei et al. (2012); Hoang and Galliero (2012); Siboulet et al. (2013); von Hansen et al. (2013) as well as the importance of molecular details in the first adsorbed fluid layers. However, the use of periodic boundary conditions (PBC) in such simulations introduces finite-size effects on the diffusion coefficient, never appreciated under confinement to date, due to the distortion of hydrodynamic flows and the spurious hydrodynamic interaction between periodic images. For bulk fluids in a cubic simulation box, the correction to the diffusion coefficient reads Dünweg and Kremer (1993); Fushiki (2003); Yeh and Hummer (2004); Tazi et al. (2012) $D(L)=D_{\infty}-\xi k_{B}T/6\pi\eta L$, with $L$ the box size, $k_{B}$ Boltzmann’s constant, $T$ the temperature, $\eta$ the viscosity and $\xi\approx 2.837$. Recently, the effect of the box shape has also been considered for bulk fluids: The components of the diffusion tensor in anisotropic boxes may in some cases be larger than $D_{\infty}$ and diverge with system size in the limit of highly elongated boxes Rozmanov and Kusalik (2012); Kikugawa et al. (2015a). These features can also be explained from hydrodynamics Botan et al. (2015); Kikugawa et al. (2015b); Vögele and Hummer (2016).

Surprisingly, despite the ever growing importance of simulations to characterize the dynamics of confined fluids, such finite-size effects have not been investigated under confinement. Here we show that the diffusion of confined fluids is not only affected by the confining distance but is also influenced within simulations by finite-size effects due to the use of PBC in the directions parallel to the surfaces. We demonstrate this fact using MD simulations of a simple fluid confined inside a slit pore. We show that the diffusion coefficient is generally larger than the value for the unconfined fluid, by a factor which is in fact significant for typical box shapes, thereby illustrating the limitations of previous simulation results. Using continuum hydrodynamics, we further obtain the scaling with system size in the limits of elongated and flat systems. As in the case of bulk fluids, the present analysis opens the way to systematic extrapolation to the limit of infinite systems (for a fixed confining distance).

We simulate a Lennard-Jones (LJ) fluid under the same conditions as in previous work illustrating finite-size effects in bulk fluids Yeh and Hummer (2004); Botan et al. (2015), namely a reduced density $\rho^{*}=\rho\sigma^{3}=0.7$ and reduced temperature $T^{*}=k_{B}T/\epsilon=2.75$ with $\sigma$ and $\epsilon$ the LJ diameter and energy, respectively. The fluid is confined between rigid planar surfaces consisting of a square lattice with a spacing of 1$\sigma$. The fluid-surface interactions are characterized by a diameter $\sigma_{FS}=\sigma$ and an energy $\epsilon_{FS}$ equal to $\epsilon$ for 3/4 of the wall atoms and $3\epsilon$ for the remaining 1/4 (see the square lattice in Figure 1). This pattern ensures the absence of slippage at the wall, as demonstrated by the velocity profile under shear (with a wall velocity of $v_{wall}=\pm 0.5$ LJ units) in Figure 1. The distance between the LJ walls is ${H}+2\sigma$, with ${H}$ the distance between the shear planes on both surfaces, while the size of the simulation box ${L}$ is the same in the other two directions. In order to assess the finite-size effects due to both the PBC and the confinement, we consider systems with ${L}$ between 6 and 80 $\sigma$ and ${H}$ ranging from 10 to 160 $\sigma$, corresponding to particle numbers from 587 to 30720.

All MD simulations are performed using the LAMMPS simulation package LAMMPS , using a time step of $10^{-3}t^{*}$, with $t^{*}=\sigma\sqrt{m/\epsilon}$ and a cut-off distance $2.5\sigma$ to compute the LJ interactions. The systems are first equilibrated in the $NVT$ ensemble during $200t^{*}$, using a Nosé-Hoover thermostat with a time constant of $t^{*}$. After equilibration, trajectories in the $NVE$ ensemble are produced for $4\,10^{4}t^{*}$ up to $1.6\,10^{5}t^{*}$ depending on system size, from which diffusion coefficients parallel to the surfaces are computed from the slope of the mean-square displacement (MSD) in the time range $5\,10^{3}~{}-~{}1.5\,10^{4}t^{*}$ ($4-8\,10^{5}t^{*}$ for ${H}\geq 80\sigma$ with ${L}=6\sigma$). For each system, reported results correspond to averages and standard errors over 16 independent runs.

Figure 2 reports the diffusion coefficient parallel to the surfaces $D_{\parallel}$ as a function of $1/{L}$ (by analogy with the scaling for a cubic box of bulk fluid) for the various confining distances ${H}$. This figure also shows the results for an infinite bulk system, equal to $D_{\infty}=0.312\pm 0.005$ LJ units in the present case Botan et al. (2015), as well as the prediction of Eq. 1 for the limit ${L}\to\infty$ under confinement. We first note that the diffusion coefficient increases with ${H}$, as expected, and that the results for large ${L}$ are consistent with Eq. 1 (except for the narrowest pore). However, there is also a dramatic influence of the periodicity along the surface which increases with pore width ${H}$. For the smaller ${L}$, the diffusion coefficient is larger than $D_{\infty}$ for an infinite unconfined fluid. For wide thin pores (large ${H}$, small ${L}$), it can be several times larger than $D_{\infty}$.

For elongated boxes (${H}>{L}$), the results of Figure 2 seem to suggest a correction to the diffusion coefficient proportional to ${H}/{L}^{2}$. Such a scaling in this regime has been predicted on the basis of hydrodynamic arguments by Detcheverry and Bocquet, who computed the enhancement of diffusion in nanometric pipes due to the thermal fluctuations of the center of mass of the fluid Detcheverry and Bocquet (2012, 2013). This increase:

[TABLE]

with $D_{\parallel}({H},\infty)$ given by Eq. 1, is related to the total fluid-wall friction coefficient $\lambda$ via the Einstein relation $\Delta D_{\parallel}=k_{B}T/\lambda$. The friction coefficient is then derived from the total force on the walls exerted by the fluid with typical velocity $v$ as $F\sim\lambda v$, which is also proportional to the fluid-wall contact area $2{L}^{2}$ and the viscous stress at the boundary $\sim\eta v/{H}$. Therefore, the increase in the diffusion coefficient scales as $\Delta D_{\parallel}\sim\frac{k_{B}T}{\eta}\frac{{H}}{{L}^{2}}$. More precisely, they obtained the following result for the present slit geometry :

[TABLE]

with $\mathcal{D}_{c.o.m.}$ the diffusion coefficient for the center of mass of the fluid. Figure 3 reports the simulation results in a dimensionless form, namely $\Delta D_{\parallel}{H}\eta/k_{B}T$ as a function of ${H}/{L}$ (here we use the bulk value of the visosity, $\eta=1.28$ LJ units Botan et al. (2015)). The collapse of the simulation results on a master curve confirms the hydrodynamic origin of the finite-size effects. While Eq. 3 accounts qualitatively for the observed behavior it is not quantitative.

In order to go beyond this simple scaling argument, we now compute the correction to the diffusion coefficient due to both confinement and PBC along the surfaces by solving the full hydrodynamic problem. Following previous studies of bulk fluids Dünweg and Kremer (1993); Yeh and Hummer (2004); Botan et al. (2015), this is achieved by computing the mobility from the solution of the Stokes equation for an incompressible fluid, $\eta\nabla^{2}{\bf v}-\nabla p+{\bf f}=0$, with ${\bf v}$ and $p$ the velocity and pressure fields. In bulk fluids, the force density ${\bf f}={\bf F}\left[\delta({\bf r})-\frac{1}{V}\right]$ with ${\bf F}$ a force, $\delta$ the Dirac distribution and $V$ the volume of the system, includes both a perturbation at ${\bf r}=0$ and a compensating “background force” which enforces the constraint of vanishing total force on the system. This compensating force provides in fact the main contribution to the system-size dependence of the diffusion coefficient.

In the present case, the broken translational invariance in the direction of confinement renders the problem more difficult. The mobility tensor defined by ${\bf v}({\bf r})=\mathbb{T}({\bf r},{\bf r}_{0})\cdot{\bf F}({\bf r}_{0})$, provides the flow at ${\bf r}=(x,y,z)$ induced by a force ${\bf F}$ applied at ${\bf r}_{0}=(x_{0},y_{0},z_{0})$ and is related to the diffusion tensor as $\mathbb{T}=\mathbb{D}/k_{B}T$. The expression of the mobility tensor $\mathbb{T}^{\infty}$ for the case of an infinite slit pore (${L}\to\infty$) with non-slip boundary conditions at the walls (${\bf v}=0$ for $z=0$ and ${H}$), previously obtained by Liron and Mochon Liron and Mochon (1975) or Swan and Brady Swan and Brady (2010), is derived in a convenient form for the present work in the Supplementary Material Sup , which also provide the details of the following calculations. It depends on the relative position along the wall $(x-x_{0},y-y_{0})$ and on both $z$ and $z_{0}$. We note that contrary to the bulk case, the mobility can be obtained in the presence of solid walls even without compensating background. The effect of the latter, which remains necessary to enforce the constraint of total momentum conservation, can be introduced separately by subtracting the average mobility over the pore volume. The average correction to the diffusion tensor then reads $\Delta\mathbb{D}=k_{B}T\left[{\left\langle\delta\mathbb{T}^{i}\right\rangle}-{\left\langle\delta\mathbb{T}^{b}\right\rangle}\right]$ where the first term corresponds to the effect of periodic images and the second to that of the compensating background.

The effect of periodic images $(m{L},n{L})$ along the surfaces can then be expressed in real space as a correction to the mobility $\delta\mathbb{T}^{i}(x,y,z)=\sum_{(m,n)\neq(0,0)}\mathbb{T}^{\infty}(x-m{L},y-n{L},z,z_{0}=z)$, with $z_{0}=z$ since the images are located in the same plane. We consider here the mobility averaged over the whole pore, which from translational invariance along the surfaces simplifies to:

[TABLE]

The average contribution of the background force, summed over all periodic images, can be written as an integral over all space:

[TABLE]

Taking advantage of the symmetry of the system, the component along the surface is computed from ${\rm T}_{\parallel}=\frac{1}{2}\left(\mathbb{T}_{xx}+\mathbb{T}_{yy}\right)$. In addition, the integral over $x_{0}$ and $y_{0}$ in ${\left\langle\delta\mathbb{T}^{b}\right\rangle}$ is conveniently computed as the value for ${\bf q}=0$ of the 2D Fourier transform $\tilde{{\rm T}}^{\infty}_{\parallel}({\bf q},z,z_{0})$. We show in the Supplementary Material Sup that the contribution of the background is equal to ${\left\langle\delta{\rm T}_{\parallel}^{b}\right\rangle}={H}/24\eta{L}^{2}$.

For elongated systems (${H}>{L}$), the discrete sum over images in Eq. 52 can be estimated by the corresponding integral, which is equal to $7{H}/60\eta{L}^{2}$, after removing the term corresponding to $(m,n)=(0,0)$, which is given by $-3\ln(1+\sqrt{2})/4\pi\eta{L}$ (see Sup for both demonstrations). Subtracting the background, we finally obtain the correction to the diffusion coefficient:

[TABLE]

The first term is consistent with the scaling argument of Detcheverry and Bocquet (see Eq. 3) and the curvature is only 10% smaller. As shown in Figure 3, it is closer to the simulation results, but the second $\mathcal{O}(1/{L})$ term, which corrects for the spurious self-interaction introduced upon replacing the discrete sum by an integral, is necessary to describe the simulation results quantitatively. The agreement of Eq. 6 with the latter is excellent. This confirms the hydrodynamic origin of the observed finite-size effects due to the PBC, in addition to the effect of confinement correctly described by Eq. 1 for sufficiently wide pores.

In the opposite regime of flat simulation boxes (${H}<{L}$), the mobility ${\rm T}^{\infty}$ decays exponentially due to the screening of hydrodynamic interactions by the walls Sup . Therefore the sum (Eq. 52) becomes negligible compared to the effect of the background (Eq. 5), so that:

[TABLE]

This explains the decrease in the diffusion coefficient with $1/{L})$ observed for the narrower pores in Figure 2. Eq. 7 is in very good agreement with the results for ${H}=20\sigma$ (within 5%) and $40\sigma$ (within 1%), but still overestimates the diffusion coefficient for $H=10\sigma$ (by $\sim 25\%$). Since the fluid is perturbed by the surfaces over 2-3 layers on each wall (see the density profile in Figure 1), it is not surprising that the present continuum hydrodynamics calculations, which neglect molecular effects, are not quantitative in this case.

Vögele and Hummer recently obtained analytical expressions for the correction to the diffusion tensor for bulk fluids in anisotropic boxes Vögele and Hummer (2016), in good agreement with numerical and molecular simulation results Botan et al. (2015). For elongated systems, the correction for the component corresponding to $D_{\parallel}$ scales as $\frac{k_{B}T}{\eta}\left[{H}/12{L}^{2}-\xi_{1}/{L}\right]$, with $\xi_{1}\approx 0.23$, i.e. the same functional form as Eq. 6 but with different numerical factors. In contrast, for flat systems the scaling goes as $\frac{k_{B}T}{\eta}\left[\ln({L}/{H})/4\pi{H}-\xi_{2}/{H}\right]$, with $\xi_{2}\approx 0.15$ and therefore diverges as ${L}\to\infty$, which is of course not the case under confinement. Therefore, despite some similarities, the present case of confined fluids is fundamentally different due to the boundary conditions at the solid-liquid interface. We finally note that, as for bulk fluids, the role of the background force enforcing the constraint of total momentum conservation is essential. In the bulk, the correction corresponding to the minimum-image cell accounts for $\approx 85\%$ of the total correction Yeh and Hummer (2004). Under confinement, the total background contribution represents a large part of the $\mathcal{O}({H}/{L}^{2})$ term in Eq. 6 for elongated systems (the nearest-image cell corresponds to the $\mathcal{O}(1/{L})$ term) Sup and almost 100% of the effect for flat systems (Eq. 7).

Ideally, the limit ${L}\to\infty$ should be obtained from simulations with ${L}\gg{H}$ in order to minimize the finite-size effects due to PBC (see Figure 2). The only exception is that of confinement down to the molecular scale, where such finite-size effects are less dramatic (see e.g. Ref. 37) due to the discreteness of the fluid and the predominance of interfacial features. However, in practice typical boxes are rather elongated or cubic ($L\leq H$) than flat, because the latter regime is computationally more expensive. The present work suggests that it is possible to minimize the finite-size effects in the other limit of elongated boxes by choosing an aspect ratio of $H/L=10\ln(1+\sqrt{2})/\pi\approx 2.8$ for which $\Delta D_{\parallel}$ cancels (see Eq. 6). It further demonstrates the hydrodynamic origin of these effects on the diffusion coefficient, as for bulk fluids, and offers with Eq. 6 an estimate of the correction to be applied to the simulation results for finite (possibly small) systems.

The present work not only applies to molecular simulations, but also to all mesoscopic hydrodynamic simulations, including Lattice-Boltzmann, Multi-Particle Collision Dynamics or Dissipative Particle Dynamics. Therefore the effects of PBC should also be taken into account in the study of confined soft matter involving colloidal particles and polymers. The analysis could also be extended to other geometries such as nanotubes, as well as to slip boundary conditions at the interface, since slippage is known to have an influence on the diffusion coefficient Lauga and Squires (2005); Saugey et al. (2005). In both cases, the same strategy can be applied using the corresponding mobility tensor for the limit of an infinite system along the surfaces.

.1 Basic hypotheses and equations

We consider the flow of a Newtonian fluid at low Reynolds number. The motion of the fluid confined between two infinite plates $z=\pm H/2$ parallel to the $xy$ plane is governed by the Stokes equation $\eta\nabla^{2}{\bf v}({\bf r})+{\bf f}({\bf r})-\nabla p({\bf r})={\bf 0}$, its incompressibility (Eq. 11) and the non-slip boundary conditions at the wall (Eq. 12). Here $\eta$ denote the viscosity, while ${\bf f}({\bf r})$ and $p({\bf r})$ denote a body force field and the associated pressure. The mobility tensor $\mathbb{T}^{\infty}({\bf r},{\bf r}_{0})$ is obtained as the solution of the so-called Stokeslet problem:

[TABLE]

with a unit horizontal force density (${F_{\shortparallel}}=1$) acting without loss of generality at ${\bf r}_{0}=(0,0,z_{0})$ in the $x$ direction. This problem was fully solved by Liron and Mochon Liron and Mochon (1975). Their solution involves multiple reflexions at the walls, complex integral representations and series expansion techniques, which render the final form difficult to exploit for the present work. We present here a more direct route, along the lines of Swan and Brady Swan and Brady (2010).

.2 Solution of the Stokeslet problem under confinement

From the Stokes and incompressibility equations (Eqs. 10 and 11), the pressure obeys:

[TABLE]

It is convenient to introduce the partial Fourier transform with respect to $x$ and $y$ defined by:

[TABLE]

with ${\bf{q}}=(q_{x},q_{y})$. The local velocity may also be decomposed into parallel ${\bf v}_{\shortparallel}(x,y,z)$ and vertical directions ${\bf v}_{z}(x,y,z)$. Using these notations, the Fourier transforms $p({\bf q},z)$, ${\bf v}_{\shortparallel}({\bf q},z)$, and ${\bf v}_{z}({\bf q},z)$ are solutions of:

[TABLE]

These second order differential equations can be solved using standard methods, yielding:

[TABLE]

where $q=||{\bf q}||$ and where all the prefactors (to be determined in the following) depend on ${z_{0}}$. In the following this dependence is dropped for the sake of clarity, but it needs to be included when performing averages over the pore in the next sections. One must manipulate quantities such as $(z-z_{0})$ and $|z-z_{0}|$ with caution, especially when computing derivatives of $\exp(q|z-z_{0}|)$ with respect to $z$. The eight constants $A,B$, $A^{x},B^{x}$, $A^{y},B^{y}$ and $A^{z},B^{z}$ are determined from the boundary conditions and the divergence-free condition. Differentiating the last equation with respect to $z$, we obtain:

[TABLE]

Inserting this expression together with $v_{x}$ and $v_{y}$ into the incompressibility condition Eq. 18, one obtains the following simple conditions:

[TABLE]

Since this holds for every position $z$, the terms in parentheses must vanish. These conditions can be written in compact form as:

[TABLE]

where we have introduced the following vectors:

[TABLE]

Inserting these results in Eq 20 and using the non-slip boundary conditions on $z=\pm H/2$ provide the remaining 6 equations that permit to determine the unknowns $A$, ${\bf A}$, $B$ and ${\bf B}$. In order to proceed as simply as possible, we rewrite the velocity field under the following form:

[TABLE]

with

[TABLE]

By introducing:

[TABLE]

and exploiting the parity of hyperbolic functions, the six non-slip conditions may be rewritten as:

[TABLE]

From these relations we express ${\bf A}$ and ${\bf B}$ as a function of $A$ and $B$:

[TABLE]

Finally, inserting these relations in Eq. 23 we obtain the expression of $A$ and $B$ as a function of ${\bf V}^{(0)+}$ and ${\bf V}^{(0)-}$:

[TABLE]

Inserting Eq .2 in Eqs 39 and 43, and substituting into Eq 20 provides, after straightforward although tedious calculations, the desired solution for $v_{x}({\bf{q}},z,z_{0})$.

The $xx$ component of the mobility tensor $\mathbb{T}^{\infty}({\bf r},{\bf r}_{0})$ then follows as the inverse 2D Fourier transform in the particular case ${F_{\shortparallel}}=1$, and the $yy$ component can be obtained similarly by computing $v_{y}({\bf{q}},z,z_{0})$ under a perturbation along the $y$ axis. The parallel component is finally ${\rm T}_{\parallel}=\frac{1}{2}\left(\mathbb{T}^{\infty}_{xx}+\mathbb{T}^{\infty}_{yy}\right)$. The computation of the inverse 2D Fourier transform is difficult, but in fact unnecessary. Indeed, for our purpose we only need to compute averages of the mobility over the pore width. This involves integrals over $z$ and $z_{0}$, as explained in the next section, as well as over $x$ and $y$ which are computed directly from the value of the 2D Fourier transform for ${\bf q}={\bf 0}$.

.3 Computation of the associated averages along $z$

As explained in the main text, two different vertical averages must be computed. For the effect of the background ${\left\langle\delta\mathbb{T}^{b}\right\rangle}$ the average velocity is taken over $z$ and $z_{0}$ independently:

[TABLE]

Similarly the $yy$ component is obtained by applying the Stokeslet in the $y$ direction and computing $v_{y}$. This double integral can be performed using the full solution of section .2, with the result in tensorial form:

[TABLE]

with $\mathbbm{1}$ the 2D identity matrix and $\hat{\bf q}$ a unit vector in reciprocal space. The parallel component, given by the half-trace of this tensor, is particularly simple:

[TABLE]

Finally, the integral over $x-x_{0}$ and $y-y_{0}$ is obtained as the value of this 2D Fourier transform for ${\bf q}={\bf 0}$, namely: $\hat{f}_{z,z_{0}}(0)=H/24\eta$. Combined with the $1/L^{2}$ factor for the average in Eq. 5 of the main text, this leads to the final result:

[TABLE]

For the effect of periodic images ${\left\langle\delta\mathbb{T}^{i}\right\rangle}$, the average is taken over $z=z_{0}$, i.e.:

[TABLE]

The solution in that case is lengthier. The final result for the half-trace reads:

[TABLE]

Here again the value of the integral over $x-x_{0}$ and $y-y_{0}$ is obtained as the ${\bf q}={\bf 0}$ value of this 2D Fourier transform, namely:

[TABLE]

III. Asymptotic results

As explained in the main text, the correction to the diffusion coefficient due to periodic boundary conditions is given by $\Delta D_{\parallel}=k_{B}T\left[{\left\langle\delta{\rm T}^{i}_{\parallel}\right\rangle}-{\left\langle\delta{\rm T}^{b}_{\parallel}\right\rangle}\right]$, where we have already computed the background correction ${\left\langle\delta{\rm T}^{b}_{\parallel}\right\rangle}$ in Eq. 48. Here we derive asymptotic expressions in the regimes of elongated and flat systems for the effect of periodic images:

[TABLE]

.4 Elongated systems ($H\gg L$)

For elongated systems, we can approximate the discrete sum in Eq. 52 by an integral, provided that we remove the contribution corresponding to $(m,n)=(0,0)$ which is not included in the sum:

[TABLE]

The first integral is computed easily using the results of the previous section. From Eq. 51, one readily obtains for the average over $z$:

[TABLE]

The second integral can be rewritten as a convolution between ${\rm T}_{\parallel}^{\infty}$ and a rectangular function with value 1 if $(x_{0},y_{0})\in[-\frac{L}{2},\frac{L}{2}]\times[-\frac{L}{2},\frac{L}{2}]$ and zero otherwise. This convolution product is conveniently computed in Fourier space, using the result Eq. 50 and the well-known transform of the rectangular function. The average over $z$ then reads:

[TABLE]

with $\hat{f}_{z_{0},z_{0}}(q)$ given by Eq. 50 and $q=\sqrt{q_{x}^{2}+q_{y}^{2}}$. We first make a change of variables, $u_{x}=q_{x}L$ and $u_{y}=q_{y}L$:

[TABLE]

with $u=\sqrt{u_{x}^{2}+u_{y}^{2}}$. Now, the regime of elongated boxes correponds to $H/L\to\infty$, so that we can approximate $\hat{f}_{z_{0},z_{0}}(\frac{u}{L})$ by the asymptotic expansion of $\hat{f}_{z_{0},z_{0}}(q)$ for $q\to\infty$, namely:

[TABLE]

which can be derived from the full expression Eq. 50. Inserting this approximation into Eq. 58, we obtain:

[TABLE]

The integral $I$ defined by the second equality can be computed analytically by writing:

[TABLE]

We then rewrite:

[TABLE]

In the last line with have introduced the error function and computed the remaining integral analytically. Gathering the result with Eq. 58, we obtain ${\left\langle\delta{\rm T}^{\rm self}_{\parallel}\right\rangle}=3\ln(1+\sqrt{2})/4\pi\eta L$, which, together, with Eqs. 52, 53 and 54 provides:

[TABLE]

Finally, the complete solution for the correction to the diffusion coefficient is obtained by substracting the contribution of the background, Eq. 48:

[TABLE]

which is Eq. 6 of the main text.

.5 Flat systems ($L\gg H$)

Ê Here we show that the mobility tensor ${\rm T}_{\parallel}^{\infty}$ decays exponentially fast with distance in real space, so that the interaction between periodic images ${\left\langle\delta{\rm T}^{i}_{\parallel}\right\rangle}$ is negligible compared to the effect of the background $-{\left\langle\delta{\rm T}^{b}_{\parallel}\right\rangle}$. To that end, we need to compute the inverse Fourier transform of $\hat{f}_{z_{0},z_{0}}(q)$:

[TABLE]

where $\hat{f}_{z_{0},z_{0}}(q)$ is given by Eq. 50 and $J_{0}$ is the zeroth-order Bessel function of the first kind. This calculation is much more involved than the previous ones, which only required the $q\to 0$ and $q\to\infty$ limits of $\hat{f}_{z_{0},z_{0}}(q)$. Liron and Mochon Liron and Mochon (1975) evaluated such integrals using the Hankel contour illsutrated in Figure 4.

Ê The real-space function is then given by:

[TABLE]

where the function $\hat{f}_{z_{0},z_{0}}(q)$ given by Eq. 50 is now understood as a function of the complex variable $q$ and is itself complex valued. The function $H_{0}^{(1)}$ is a Hankel function given by $H_{0}^{(1)}(z)=J_{0}(z)+iY_{0}(z)$, with $Y_{0}$ the zeroth-order Bessel function of the second kind. The function $\hat{f}_{z_{0},z_{0}}(q)qH_{0}^{(1)}(qr)$ has no pole at the origin, but an infinite number of poles corresponding to the $\coth$ function (except at the origin) and additionnal poles denoted $s_{n}$ corresponding to the solutions of the transcendental equation $\sinh(s)^{2}=s^{2}$. For large $n$ the asymptotic behavior of these poles is given by:

[TABLE]

The net result appears as an absolutely convergent series of functions involving modified Bessel functions $K_{0}(n\pi r/H)\propto\displaystyle\sqrt{\frac{H}{r}}\exp(-n\pi r/H)$ coming from the residues arising from the $\coth$ function. An additionnal series of Bessel functions $|H_{0}^{(1)}(s_{n}r/H)|\propto\displaystyle\frac{\pi}{2}\sqrt{\frac{H}{r}}\exp(-y_{n}r/H)$ comes from the residues associated with the $s_{n}$ poles. As in both cases, the imaginary part of the poles behaves as an arithmetic sequence, the associated series can be approximated by their first term for large $r/H$. Overall, the asymptotic behavior for large $r/H$ is given by:

[TABLE]

i.e. an exponentially decreasing correction.Ê This result can in fact be recovered using equation (49) of Liron and Mochon Liron and Mochon (1975), by computing the required trace that suppreses the long-range contribution, and performing the associatedÊ$(z_{0},z_{0})$ average term by term in the resulting series. The exponential decay in real space implies that the sum over periodic images is dominated for flat systems ($L\gg H$) by the nearest images, so that this sum also decays exponentially fast with $L/H$. Overall, confinement screens the hydrodynamic interactions between the periodic images and the corresponding contribution ${\left\langle\delta{\rm T}^{i}_{\parallel}\right\rangle}$ is negligible compared to the effect of the background $-{\left\langle\delta{\rm T}^{b}_{\parallel}\right\rangle}$.