1 Introduction

Hydrodynamic coupling and collective diffusivity can be significantly augmented by the proximity of microstructures to a surface. This is often relevant to transport phenomena in many applications, such as the electrophoresis of macromolecules and colloids in capillaries ^[1-6], processing of emulsions and polymers in microfluidic devices ^[7-11], or the behaviors of active suspensions such as swimming microorganisms near surfaces ^[12-13]. We develop stochastic computational methods to capture confinement effects of microstructures within channel geometries with no-slip walls. Our approach is based on the stochastic Eulerian-Lagrangian method (SELM) which provides tractable ways to incorporate thermal fluctuations into approximate descriptions of the hydrodynamic coupling between microstructures^[14]. In the SELM, fluctuating hydrodynamic equations similar to those introduced by Landau and Lifshitz^[15] are coupled and exchange momentum with microstructure conservation equations. The SELM framework provides criteria to ensure that the continuum stochastic description and related stochastic discretizations yield results consistent with statistical mechanics^[14]. Related computational methods for fluctuating hydrodynamics have also been introduced in Refs. [16]-[20].

We present the SELM approach for the channel geometry in Section 2. We develop stochastic numerical methods using a staggered finite volume discretization for the fluid velocity and pressure in Subsection 3.1 and Appendix A. To obtain a consistent discretization of the stochastic driving fields^{[14, 16]}, we impose an exact fluctuation-dissipation balance on our stochastic numerical methods taking into account the augmented dissipative properties of the discrete operators relative to their continuum differential counterparts in Subsection 3.1. We show that our numerical fluctuation-dissipation balance principle combined with an exact projection method for the incompressibility is sufficient to ensure results consistent with statistical mechanics in Subsection 3.2. In practice, a significant challenge is to generate efficiently the stochastic driving fields with the required covariance structure obtained from the fluctuation-dissipation balance condition. We show how a method with the cost O(Nlog(N)) can be developed for our discretization based on fast Fourier transforms (FFTs) to generate efficiently the driving fields in the presence of the no-slip walls in Section 3.3. We validate the computational methods in Section 4 by performing studies for the stochastic field generation in Subsection 4.1 and for the Brownian motion of a particle diffusing in a harmonic potential in Subsection 4.2. As a demonstration of our computational methods, we consider ellipsoidal particles within a channel and investigate the effects of confinement on both active hydrodynamic responses and passive diffusivity in Section 5. We first compare our computational model of ellipsoidal particles with those of analytic results for the translational and rotational mobilities in Subsection 5.1. We then investigate how the proximity of a particle to a wall within the channel affects the mobility of the ellipsoidal colloids in response to active forces or torques in Subsection 5.2. We then demonstrate through stochastic simulations that our fluctuating hydrodynamic methods capture the role of proximity to the no-slip walls as manifested in the empirical diffusivity of particles. We show that our stochastic methods yield diffusivities in close agreement in accordance with the Stokes-Einstein relations with the predictions from the active mobility studies in Subsection 5.4. We then make comparisons with conventional Langevin dynamics in Subsection 5.5. We study the collective dynamics of a large number of colloidal particles within the channel and their density relaxations. We find that our fluctuating hydrodynamics approach yields results having significant differences with simulations performed with conventional Langevin dynamics that neglects the hydrodynamic coupling. We conclude by investigating these differences by considering for a large number of interacting particles their collective drift-diffusion dynamics and associated correlation functions. Overall, we expect our introduced stochastic computational methods to be broadly applicable to problems involving confinement effects of microstructures subject to hydrodynamic coupling and thermal fluctuations.

2 Fluid-structure interactions subject to thermal fluctuations

Our description of the fluid-structure interactions subject to thermal fluctuations is based on the SELM^[14]. In the SELM, the microstructures exchange momentum with a fluctuating fluid to account simultaneously for hydrodynamic coupling and thermal fluctuations. The SELM provides a way to incorporate thermal fluctuations into the widely used approaches for approximating the fluid-structure interactions^[14]. We extend the SELM approach to account for no-slip walls when the microstructures and fluid are confined within a channel geometry. For the microstructure dynamics, we use

(1)

(2)

These are coupled to the incompressible fluctuating hydrodynamic equations,

(3)

(4)

We consider these on the domain Ω = [0, L_x]×[0, L_y]×[0, L_z] subject to the no-slip boundary condition on z=0 and z=L_z,

(5)

X denotes the collective vector of all of the degrees of freedom of the microstructure, v is the microstructure velocity, and m is the microstructure excess mass^{[14, 21]}. The fluid velocity is denoted by u, the fluid density is denoted by ρ, and the dynamic viscosity is denoted by μ. The pressure acts as a Lagrange multiplier to enforce the incompressibility constraint given in Eq. (4).

The operators Γ and Λ serve to couple the microstructure and fluid dynamics. The Γ operator serves to provide a local reference velocity from the fluid against which the microstructure velocity is compared. The term -Υ(v-Γu) acts as an effective drag force on the microstructure when its velocity differs from that of the surrounding fluid. Υ is assumed to be a positive-definite operator. The Λ operator is introduced to account for the drag of a moving particle's force equal-and-opposite effect on the fluid. The Λ operator converts a microstructure force into a corresponding force density in the fluid equations in Eq. (3). A particularly important property to obtain consistent results in the mechanics is that the coupling operators satisfy an adjoint condition Λ = Γ^{T[14, 21-22]}.

To account for thermal fluctuations, we let F_thm and f_thm be Gaussian random fields that are δ-correlated in time and have zero mean ^[23-24]. We determine the spatial covariance using the fluctuation-dissipation principle of statistical mechanics^{[14, 25]},

(6)

(7)

(8)

The adjoint condition Λ = Γ^T has been shown to play an important role when introducing the stochastic driving fields^[14]. Throughout our discussion, our stochastic differential equations should be given the Itô interpretation in stochastic calculus ^[23-24].

A widely used approximation to obtain a tractable description of the fluid-structure interactions is the immersed boundary method^[22]. This corresponds to the specific choice of coupling operators,

(9)

(10)

The kernel function η(z) is chosen to be the Peskin δ-function developed in Ref. [22]. The Peskin δ-function has important numerical properties that ensure a good approximation that there is the translational invariance of the coupling despite the breaking of this symmetry by the discretization lattice of the fluid^{[16, 22, 26]}. A smooth kernel that integrates to one and vanishes outside a finite volume (finite support) is chosen instead of a Dirac δ-function to ensure a model in which the mobility of even individual point particles has a finite effective hydrodynamic radius within the fluid^[16]. For the confined channel geometry and any particle near the boundary, the kernel function should be normalized over the part of the spatial domain that is within the channel. Since other choices for the fluid-structure coupling can also be useful, we only make use of the generic properties of Γ, Λ throughout the mathematical discussions. For coupling the fluid and microstructure, we use the Stokes drag,

(11)

where R is a length-scale which we take to be comparable to the support of the kernel function η(x). is the identity operator. Provided that Υ is sufficiently large, the precise value does not play a particularly central role in the dynamics^{[14, 21]}.

3 Semidiscretization: staggered finite volume method

We discretize the system using a finite volume approximation on a staggered grid for the velocity and pressure. The velocity components are represented at the cell faces and the pressure at the cell centers (see Fig. 1). A particular advantage of the staggered grid discretization is that the mass transport can be modeled naturally between the cells using fluxes based on the velocity at the cell faces. This is also useful in approximating the divergence-free incompressibility constraint since it arises in the continuum fluid mechanics from the continuity equations and the requirement that a uniform mass density remains constant under the fluid flow. Using this interpretation in the finite volume setting yields a constraint on the discrete velocity field expressed in terms of a discrete divergence operator which can be imposed exactly using the cell centered pressure and its discrete gradient. We approximate a finite volume formulation of the system to derive a finite difference method with a set of discrete operators playing roles very similar to their continuum differential counterparts. We shall show that, in the stochastic setting, these features have important implications for the statistical mechanics of our stochastic numerical methods. Related finite volume schemes for fluctuating hydrodynamics have been introduced ^{[16-17, 19]}.

We use for our discretization a computational domain decomposed into N_x× N_y × N_z cells indexed by m = (i, j, k). We decompose the fluid velocity u = (u⁽¹⁾, u⁽²⁾, u⁽³⁾) = (u, v, w) into its Cartesian components at each cell face. We index the fluid velocity in the x-direction by half indices u_{i±1/2, j, k} and similarly for the other directions. The pressure p is defined at the cell centers and indexed by (i, j, k) as p_{i, j, k}. A notable feature of the staggered grid discretization is that only one of the components of the velocity field is stored at each face center. This is the minimal information required to compute a well-defined divergence operator on the fluid velocity. This means that, for a mesh cell, the x-faces only store the x-component of the fluid velocity, the y-faces only store the y-component of the fluid velocity, and the z-faces only store the z-component of the fluid velocity. The no-slip boundary conditions are imposed by requiring for mesh cells bordering the wall that the z-faces have the zero fluid velocity for the component w_{i, j, -1/2}=w_{i, j, Nz-1/2}=0. The no-slip boundary conditions are also enforced through the definition of the discrete operators which are truncated to have zero weights for any values involving points outside the domain. This has the same effect as setting all the velocity components to be zero for any points outside the channel domain. The boundaries in the x-direction and y-direction are treated as periodic with the condition on velocity components u_{i, j, -1/2}=u_{i, j, Nz+1/2}, v_{i, j, -1/2}=v_{i, j, Ny+1/2}. To simplify the notation, we shall often use the vector-index notation u_I, p_I with I = (i, j, k). We let Ω_c denote the collection of cell indices within the interior of the domain and let ∂Ω_c denote the cell index on the boundary. Similarly, we let Ω_f denote the collection of face index within the interior of the domain and let ∂Ω_f denote the index on the boundary.

We use the finite volume interpretation to define the discrete divergence operator by

(12)

where e_d is the standard basis vector in the d-direction. We define the discrete gradient operator using the negative adjoint of the discrete divergence operator which yields

(13)

We define the face centered Laplacian as

(14)

We define a cell centered Laplacian using the discrete divergence and gradient operators,

(15)

We use our discrete operators to approximate Eqs. (1)-(4) to obtain the semidiscretization

(16)

(17)

(18)

(19)

These are subject to the boundary conditions

(20)

(21)

(22)

(23)

Using the finite volume interpretation of our discretized system, we associate the energy

(24)

where Δx_m³ denotes the volume of the cell with the index m. We integrate Eqs. (16)-(19) numerically in the time using the Euler-Maruyama method^[27].

An important consideration for the semidiscretized system is that the discrete operators have different dissipative properties from their continuum differential counterparts. It becomes important in the discrete stochastic setting to take this into account in the choice of the stochastic driving fields F_thm and f_thm to ensure the appropriate thermal fluctuations consistent with statistical mechanics^{[14, 16]}.

3.1 Stochastic driving fields for the discretization

To account for the thermal fluctuations in the discrete setting, we must approximate the stochastic driving fields in Eqs. (1)-(3). We must also take into account in the discrete setting the role of the incompressibility constraint. Another important consideration is that the dissipative properties of the discrete operators are often significantly different from their continuum differential counterparts. The properties of the specific discretization and how constraints are handled have important implications for how fluctuations propagate through the degrees of freedom of the discretized system^{[14, 16]}.

We treat the incompressibility constraint by imposing exactly the discrete divergence-free condition (17) and using the projection method^[28]. This approach allows for the semidiscretization in Eqs. (16)-(19) to be expressed as

(25)

(26)

(27)

The projection operator is

(28)

We approximate the stochastic driving fields by imposing on the discretization the following fluctuation-dissipation condition^{[14, 16, 25]}:

(29)

where F = [f_thm, F_thm] represents the stochastic terms in the system of Eqs. (16)-(19), and is the grand dissipative operator given by

(30)

An important part of our discretization method for the stochastic fields is the choice of covariance C we impose for the equilibrium fluctuations of the discrete system. We introduce an energy based on our finite volume interpretation in Eq. (24). Using the Gibbs-Boltzmann statistics associated with our discrete energy, we impose the covariance

(31)

Now that and C are specified, our discrete fluctuation-dissipation condition given in Eq. (29) provides our approximation of the stochastic driving fields F through the covariance,

(32)

By using the fluctuation-dissipation condition (29) to discretize the stochastic driving fields, we have taken into account the properties of the specific choice of discretization and how this influences the propagation of fluctuations throughout the system^{[14, 16]}. This approach drawing on insights from statistical mechanics also circumvents a number of potentially subtle issues in how directly to interpret and approximate the stochastic fluid equations which have solutions only in the generalized sense of distributions^{[14, 16]}.

We remark for the sake of our subsequent presentation and calculations that the precise form of the adjoint condition for Λ and Γ depends on the utilized representation of the operators and the inner-product. When representing the operators as standard matrices and using the standard matrix-vector inner product, the adjoint condition takes the form of Λ^T = Δx_n^-3Γ, and hence the apparent asymmetric factor of Δx_n^-3 appears in Eq. (32). Because of this interpretation, this in fact does not pose an issue in practice.

3.2 Statistical mechanics of the discretization subject to the incompressibility constraint

The stochastic driving fields are discretized without explicit consideration of the incompressibility constraint. We perform the analysis to explore the implications of this constraint on the statistical mechanics of the stochastic dynamics. In our approach, the incompressibility is imposed each time step using the Lagrange multiplier which is obtained from the discrete pressure p by solving

(33)

(34)

In practice to solve efficiently these equations, we use FFTs adapted to the channel geometry. We use standard FFTs in the periodic x-and y-directions ^[29-30]. To account for the Neumann boundary conditions in the z-direction, we use a fast cosine transform (FCT)^[30]. By substituting the solution into Eq. (16), we obtain the same result as applying the projection operator given in Eq. (25). In this manner, the fluctuating hydrodynamics satisfies exactly in the discrete setting the divergence-free incompressibility condition ·u = 0.

Under the projected stochastic dynamics, we show invariance of the Gibbs-Boltzmann distribution,

(35)

where Z denotes the partition function. The Gibbs-Boltzmann distribution under the projected stochastic dynamics satisfies

(36)

with the probability flux on the phase-space

(37)

The stochastic driving field when subjected to the projection operator yields the term . Hence, the covariance for the stochastic driving field when subjected to the projection operator is given by . The divergence of the flux can be expressed as

(38)

with

(39)

The energy of the discretized system we introduced in Eq. (24) has gradients

(40)

(41)

(42)

In our notation to simplify the expressions, we have suppressed explicitly denoting the fields on which the operators act with this information easily inferred from Eq. (25).

To show the invariance of Ψ_GB under the stochastic dynamics, we show that each of the terms A₁, A₂, A₃, A₄ is zero. The expression A₁ is closely related to the time derivative of the energy under the non-dissipative dynamics of the microstructures. We find after direct substitution of the energy gradients and cancellations that A₁ = 0. The term A₂ has important implications for the transport of probability mass on the phase-space of (u, v, X). It can be interpreted as the phase-space compressibility associated with transport by the vector field of the non-dissipative dynamics of the system (0, -m^-1∇_XΦ, v). A₂ is zero since each term depends on distinct degrees of freedom from those appearing in the gradient being taken. This is closely related to the Hamiltonian structure of our non-dissipative dynamics which have conjugate configuration and momentum degrees of freedom^[14]. The term A₃ arises from fluxes driven by multiplicative noise in the stochastic driving fields of the system^[23]. In the present setting, the multiplicative noise has a rather special form in which only the particle configuration X plays a role in modulating the noise covariance . This enters through the terms Γ and Λ which appear only in and (see Eq. (32)). Since each term depends on distinct degrees of freedom from those appearing in the gradient being taken, we have the term A₃ = 0.

This leaves A₄, which accounts for the balance between the dissipation and stochastic fluctuations of the system. An important issue is how the projection operator enforcing incompressibility impacts this fluctuation-dissipation balance. The first covariance term in A₄ in Eq. (39) can be expanded using properties of our projection operator and discrete energy to obtain

(43)

We have used the property that our discrete defined in Eq. (28) is an exact projection operator so that We have also used that the incompressibility constraint is imposed exactly so that the fluid velocity field at any given time satisfies By a similar argument, we have

(44)

(45)

We have used throughout that the projection operator when extended to all of the degrees of freedom of the system does not affect directly the microstructure configuration or velocity so in fact v = v and X = X. The components of associated with X are all zero, i.e., Using these results in Eq. (39) and canceling common terms yield that A₄ = 0. This shows that ∇·J = 0.

We have shown that the Gibbs-Boltzmann distribution is invariant under the projected stochastic dynamics given in Eqs. (25)-(27). These results show that it is sufficient to ensure thermal fluctuations consistent with statistical mechanics by using our combined approach of (ⅰ) imposing exactly a fluctuation-dissipation condition to obtain a discretization of the stochastic driving fields^{[14, 16, 31]} and (ⅱ) imposing the incompressibility constraint with an exact projection ^{[14, 16, 21]}. The fact that the first condition ensures the discretization properly balances the stochastic driving fields with the dissipative properties of the discrete operators. The second condition using an exact projection ensures that the incompressibility constraint does not introduce new discretization artifacts which adversely affect the propagation of thermal fluctuations through the discretized system. The general principles we have presented here, and in the prior works^{[14, 16, 31]}, can be used broadly to obtain discretizations for incompressible fluctuating hydrodynamics and other spatially extended stochastic systems.

An interesting point to remark is that the no-slip boundary conditions did not require any explicit consideration in our analysis above. However, they were tacitly included through the precise definition of the discrete dissipative operators and the particular form taken by the projection operator . The effective action of these operators on the degrees of freedom representing the state of the system in the interior of the domain is influenced by the boundary conditions. For instance, the stencils of the discrete operators change at locations adjacent to the no-slip boundary. An important feature of our discretization approach for the stochastic driving fields is that our fluctuation-dissipation condition naturally handles the augmented behaviors of the dissipative operators near the boundaries caused by the no-slip conditions. These features highlight the utility of our approach of using statistical mechanics principles to approximate the stochastic driving fields to obtain practical numerical methods for spatially extended stochastic equations^{[14, 16, 31]}.

3.3 Efficient generation of the stochastic driving fields

We have shown that imposing an exact fluctuation-dissipation balance on the semidiscretized stochastic dynamics provides a derivation of effective stochastic driving fields in the discrete setting that yields thermal fluctuations consistent with statistical mechanics. An important challenge in practice is to generate efficiently these discrete stochastic driving fields with the specific spatial covariance given in Eq. (32). A natural approach for generating multivariate Gaussians ξ with a specified covariance structure 〈ξξ^T〉 = G is to perform a factorization G = HH^T and use standard Gaussian variates η with 〈ηη^T〉 = I to generate ξ = Hη. This follows since 〈ξξ^T〉 = 〈Hηη^TH^T〉 = HH^T = G. For this to be effective, a key issue is to have a factor H whose action on a vector can be computed efficiently and to determine efficiently the factor H itself from G. Unfortunately, the most straightforward method of using Cholesky factorization to obtain H from G has a computational cost of O(N³) and typically yields a dense factor which would incur a cost of O(N²) each time we generate the stochastic driving fields. Here, N is the total number of mesh degrees of freedom which for problems in three spatial dimensions would be rather large.

We shall take another approach to generate the random variates with the computational cost O(Nlog(N)) by using special properties of the stochastic dynamics and the discrete operators. We first factor G into a form that decouples the microstructure and fluid degrees of freedom. This is accomplished by expressing the stochastic driving fields f_thm and F_thm in terms of g_thm=f_thm + ΛF_thm as in Ref. [14]. This has the convenient property that the spatial covariances become

(46)

(47)

(48)

A particularly important property is that F_thm and g_thm are now decorrelated Gaussian variates that can be generated independently. The original stochastic driving fields can be recovered by using that f_thm=g_thm -ΛF_thm. In the case that Υ is a scalar as it is assumed in our current presentation, the variate F_thm can be generated trivially with the computational cost O(N). This reduces the problem to that of generating efficiently the variate g_thm.

To generate g_thm, we shall use that its covariance structure G given in Eq. (46) can be diagonalized into the form D = PGP^T, where D is diagonal, and P provides a unitary change of basis. We generate the random field using , where η is a complex-valued Gaussian random variate with 〈ηη^T〉 = I. This requires that we determine the diagonal factor D and the necessary transforms to compute the action of P^T. Since η is complex-valued, to fully determine its statistics also requires that we determine the covariance between the real and imaginary components and also impose conditions that ensure a real-valued stochastic field g_thm.

For the specific spatial discretization introduced for the channel geometry in Section 3, we have that the discrete Laplacian subject to the no-slip boundary conditions is diagonalized by the transform P = S_zF_yF_x, where S_z denotes for the z-direction the fast sine transform (FST), and F_x and F_y denote for the x-and y-directions the standard FFTs ^[29-30]. The Fourier symbols yielding D for our discrete Laplacian can be computed with the computational cost O(Nlog(N)), which needs to be performed only once for the discretization mesh. We use in our algorithms the complex exponential form of the FFTs. The generation of the stochastic driving fields is accomplished by computing , where the complex-valued η satisfies 〈ηη^T〉 = I.

To make this generation procedure work in practice, there are some important additional considerations. Since the complex values are used in the Fourier-space, there are multiple ways to satisfy the condition 〈ηη^T〉 = I. In fact, to determine fully the statistics of the Gaussian variate η also requires specification of the other covariance components between the various real and imaginary components of η. A related issue is that the generated stochastic field g_thm must be real-valued. This latter requirement provides additional important conditions on η. We now discuss in detail for the discretization our procedure for generating such a complex-valued Gaussian.

Corresponding to the transforms F_x, F_y, and S_z, we index the Fourier modes using the wave-vector notation k = (, m, n). We define the conjugate mode to k as the mode with the index k = σk = (-, -m, n). We define the Hermitian transpose as η_k^H =η_k^T which we also refer to as the conjugate transpose. For the Fourier modes, under the inverse transform P^T = F_x^TF_y^TS_z^T, to yield a real-valued field requires the conjugacy condition η_k = η_σk. From these considerations, we have for Fourier modes corresponding to a real-valued field that η_k^T = η_k^T = η_σk^T. The condition 〈ηη^T〉 = I takes the form for the individual modes 〈η_k η_k^T〉 = I, where I for the modes denotes the three by three identity matrices. The real-valued conjugacy condition gives that η_k = η_σk. By letting k''= σk', we have the two conditions,

(49)

(50)

δ_{m, n} denotes the Kronecker δ-function. We now show for η that these two conditions fully determine the statistics of the real and imaginary components η = α + iβ. The conditions take the form of

(51)

(52)

Equating the real and imaginary parts in each expression gives the conditions

(53)

(54)

(55)

(56)

By linearly combining the equations, we yield the conditions,

(57)

(58)

(59)

The last condition shows that the real and imaginary components of the modes of the Gaussian η should be generated to be statistically independent. The condition that the generated stochastic field is real-valued has the important consequence that the random Gaussian variate for the real part α_k must always be exactly the same as the value α_σk, and the random Gaussian variate for the imaginary part β_k must always be exactly negative of the value β_σk. In the special case of self-conjugate modes k = σk, the conditions require that the imaginary part is zero, β_k = 0, and the contributions to the random field are only made by the real part 〈α_kα_k^T〉 = I. This provides the details of how to generate the modes η_k in the Fourier space to obtain the required complex-valued Gaussian random variate η. This can be accomplished in the Fourier space with a computational cost of O(N). By performing the inverse transforms, this provides an efficient method with the computational cost O(Nlog(N)) to generate the stochastic driving field g_thm. We remark that the effects of the no-slip boundary conditions are tacitly taken into account in our stochastic field generation through the specific form taken by the Fourier symbols in the diagonal operator D and the form of the inverse transforms used.

4 Validation studies 4.1 Covariance obtained from the stochastic field generator

To validate our methods for generating the stochastic driving field g_thm in Section 3, we consider the discretization of the channel domain with 20×20×20 mesh cells in each direction introduced in Section 3. We estimate empirically from generated samples of the stochastic driving fields the spatial covariance structure in the real-space . This is predicted to have the entries of the covariance matrix . These entries correspond proportionally to our face-centered approximation of the Laplacian by central differences. The correlation structure is sparse and involves a positive entry when m = m₀ and a negative entry for each m corresponding to the nearest neighbors in each of the Cartesian directions. We show the spatial correlation estimated empirically for a cross-section of the mesh for each of the velocity components in Fig. 2. We find good agreement between the predicted covariance structure and those obtained from our stochastic field generation method introduced in Subsection 3.3.

4.2 Brownian motion of a particle diffusing in a harmonic potential

We validate our computational methods for the microstructure dynamics coupled using fluctuating hydrodynamics by considering the Brownian motion of a particle diffusing in the harmonic potential . We use the spring stiffness K = k_BT, the particle mass m = 19 200 u, and the temperature T = 298.15 K (k_BT = 2.479 nm²· u·ps^-2). The computational domain is taken to be 30 nm×30 nm×30 nm resolved with a grid having cells with the mesh-width 1 nm. For this choice of parameters, the Gibbs-Boltzmann distribution ρ(X) = (1/Z)exp(-Φ(X)/(k_BT)) predicts a standard deviation in the position of =1 nm. The Maxwellian for the velocity fluctuations predicts a standard deviation = 11.36 nm/ns^[25]. We perform simulations with Δt = 0.15 ns for 50 000 time steps and estimate the position and velocity distributions. We find that our stochastic numerical methods yield an appropriate effective temperature for the microstructure and that the fluctuations in configuration agree well with the Gibbs-Boltzmann distribution (see Fig. 3).

5 Simulations of ellipsoidal colloids and confinement effects

To simulate the diffusion and hydrodynamic responses of ellipsoidal colloids subject to confinement effects within a channel, we develop a model based on control nodes connected by harmonic bonds with a specified rest-length as shown in Fig. 4. The control node distribution is obtained by starting with a scaled icosahedron and bisecting the edges recursively. We refine the edges of the icosahedron so that the bond length is about half of the grid size. Similar types of models have been developed ^[32-34]. We choose the strength of the bonds sufficiently strong that there is a minimal deformation of the particle shape when subject to hydrodynamic flow and applied forces ensuring a rigid-body-like response. We also choose a sufficiently strong coupling parameter for the fluid and control particles to ensure a no-slip-like response (for details, see Table 1). We consider how our model behaves in practice for ellipsoidal particles of different aspect ratios by making comparisons with the classical analytic results of Perrin and others ^{[1-2, 35-36]}.

5.1 Translational and rotational mobility

We consider prolate ellipsoidal colloids with the semi-major axis a and semi-minor axes b and c (b=c). For such colloidal particles in a bulk Newtonian fluid, the translational mobility has been analytically predicted by Perrin and others as ^[35-36]

(60)

where . In this expression, we take into account the correction to the typo that was found in the original paper by Perrin in the statement of Eq. (41) ^[35-36]. For the rotational mobility of prolate ellipsoids in response to an applied torque, Chwang and Wu^[37] analytically predicted

(61)

where

We compare the hydrodynamic responses of our colloidal model with these results. To use our computational methods but to minimize the effects of the wall in these studies, we consider computational domains where the channel has a width of about 10 times the particle radius. To probe the translational mobility of our model, we apply a unit force over the surface of our particle and measure from our simulations the resulting steady-state velocity of the ellipsoidal particle. To probe the rotational mobility of our model, we apply a unit torque as an averaged force distribution over the surface of our particle and measure the angular velocities by averaging moments over the particle surface. We find that the translational and rotational hydrodynamic responses of our particle model agree quite well with the analytic predictions when using an effective hydrodynamic radius for our particle of a=5.45 nm (see Fig. 5). We remark that this is consistent with how the SELM handles the control point coupling to the mesh through the averaging operator Γ which blurs the fluid-structure interface on the scale of the mesh width Δx≈ 1 nm (see Eq. (9)). The translational mobilities show very good agreement for all of the aspect ratios considered. The primary discrepancies occur for the rotational mobilities as the aspect ratio becomes small. The parameters for these studies are given in Table 1.

5.2 Effects of confinement on mobility and hydrodynamic interactions of particles 5.2.1 Mobility of spherical particle within a slit channel

To further validate the methods and our particle model, we perform simulations of the motion of a spherical particle when confined between two stationary parallel walls. We compute the particle in response to an applied force either in a direction parallel to the wall to obtain the mobility M_‖ or in a direction perpendicular to the wall to obtain the mobility M_⊥. Our simulation results are compared with the linear superposition approximation (LSA) for mobility of a spherical particle in a slit channel ^[38-39]. The LSA predicts a mobility of

(62)

(63)

where a is the particle radius, L is the channel width determined by the separation between the walls, and is the distance of the center of mass of the particle to the nearest channel wall. We set the reference mobilities M_‖⁰ and M_⊥⁰ based on the simulated mobility of the spherical particle at the center of the channel. We compare our simulations of the channel confined spherical particle model with the LSA theory in Fig. 7 for setup as in Fig. 6. We find the SELM approach yields results with close agreement with the predictions of the LSA. We find that M_‖⁰ agrees with the free-space mobility of a spherical particle of the effective hydrodynamic radius = 5.45 nm. We remark that this is the same value we used previously in our Rotne-Prager comparisons in Subsection 5.1. These results provide a validation of our numerical methods, indicating that they capture well the physics of particles confined within a slit channel. We remark that the LSA theory is expected to be most accurate away from the channel walls where lubrication effects are not dominant in the particle-wall interactions.

5.2.2 Mobility of ellipsoidal particles within a slit channel

We explore the role that confinement effects play a role in mobility and hydrodynamic interactions between the particles. We investigate responses for a spherical particle with the radius r=5 nm and an ellipsoidal particle with a=5 nm and b=c=2.5 nm. We consider a channel with walls separated by 30 nm in the z-direction and a span of 120 nm in the periodic x-and y-directions. We resolve this channel geometry using a discretization with the mesh-width 1 nm in each direction.

We investigate how the translational and rotational mobilities change as the colloidal particle occupies distances to one of the walls ranging from 10 nm to 2 nm. We consider for the ellipsoidal particle the two cases when (ⅰ) the semi-major axis is parallel to the channel wall and (ⅱ) the semi-major axis is perpendicular to the channel wall. As a consequence of the symmetry of the channel system under rotation around the z-axis, we only need to consider the responses in the x-and z-directions.

We use our colloid model and computational methods to study the hydrodynamic responses of particles as they approach one of the channel walls (see Fig. 8). We consider the translational and rotational mobilities obtained from the active responses of particles to an applied force or torque. We find that the translational mobility is impacted most significantly by the particles proximity to a channel wall. In contrast, the rotational mobility is only affected by the channel wall when the particle is rather close. The active responses as characterized by the mobility M have important implications also for the passive rotational and translational diffusivities D of particles within the channel through the Stokes-Einstein relation D = 2k_BTM^[25]. These results indicate that the stochastic methods will also capture how the diffusivity of a particle is augmented by its proximity to the wall.

5.3 Reduced model of ellipsoidal particles as trimers

In our studies, we expect that many of the hydrodynamic responses of the colloidal particles can be captured to a good approximation by a less detailed geometric description using only a few control nodes as in Fig. 9 ^[32-34]. To improve the computational efficiency in simulations, we use a reduced description that replaces the full geometric approach introduced in Section 5 by instead a trimer of three control nodes linked by harmonic springs of non-zero rest-length. While this reflects the approximate responses of an ellipsoid for one particular aspect ratio, this can also be augmented by adding control nodes arranged either linearly or in a cluster to obtain the approximate responses for other aspect ratios. When comparing the responses of our trimer model with three control nodes that span about 2 mesh-widths of the discretization mesh with the results of Perrin's equation (60), we find the trimer approximates a prolate ellipsoid with the aspect ratio of p = c/a = 0.6. We use this choice in our subsequent studies.

5.4 Hydrodynamic interactions and pair diffusivity: translational and rotational pair mobility

We probe the effects of confinement on both the active responses mediated by hydrodynamics between a pair of particles and the passive diffusive correlations in the joint Brownian motions of a pair of particles. We investigate the pair mobility M by studying for a force applied to one particle the hydrodynamic response of the other particle. By the Stokes-Einstein relations, the pair diffusivity of particle D is predicted to be

(64)

In the deterministic case of two particles in response to an applied load, the grand-mobility tensor M can be expressed in terms of the components of the particle configuration X = (X₁, X₂) as

(65)

X₁ gives the configuration of the first particle, and X₂ gives the configuration of the second particle. In the current setting, we are most interested in the block diagonal entry M₁₁ and the block off-diagonal entry M₁₂. M₁₁ characterizes how a force acting on a particle affects the motion of this same particle. M₁₂ characterizes how a force acting on one of the particles affects the motion of the other particle. We determine these responses in the deterministic setting (no thermal fluctuations) empirically by applying such a force and measuring the resulting velocity in our computational simulations (see Appendix B).

To estimate empirically the pair diffusivity D of our particles, we use

(66)

where ΔX = X(δt) -X(0). The components in which we are most interested are the diagonal terms of D₁₁ and D₁₂ which correspond to how the diffusive motion of one of the particles is correlated with the other in the same direction.

To investigate the role of confinement in particle interactions, we consider both the active force responses and the passive diffusivities and how they depend on both the proximity of the particles to one of the channel walls and the distance of separation between the particles. For the wall proximity, we consider the two specific cases when the particle is located at the positions z=2 nm and z = 10 nm.

We find that confinement plays a significant role both in the active hydrodynamic interactions between particles and in the passive diffusivity of their correlated Brownian motion. We find that our stochastic numerical methods yield results in close agreement with the Stokes-Einstein relations given in Eq. (64), as shown in Fig. 10. We consider the role of proximity to the channel wall both on the active hydrodynamic responses in the interactions between a pair of ellipsoidal particles (the aspect ratio p = c/a = 0.6) and on the passive diffusivities in the correlated Brownian motions of the two particles. We find confinement significantly impacts both the active responses and passive diffusivity. We also find that the stochastic computational methods agree to a good approximation with the predictions of the Stokes-Einstein relations D = 2k_BTM^[25]. We remark that we think the slight discrepancy appearing in our results comes from the inherent issues associated with how to attribute a mobility/diffusivity to the discrete particle models. In the deterministic setting, we estimate a mobility by measuring the translational and rotational velocities in response to an applied force/torque realized by a symmetric force distribution over the surface control points of the particle (for details, see Appendix B). The diffusivity is measured by averaging the translational and rotational motions of the collection of surface control points. As a consequence of these differences in the estimators, we expect some small systematic discrepancies.

5.5 Collective drift-diffusion of rods

As a further demonstration of the approaches, we consider the collective drift-diffusion dynamics of the rods when subject to a driving field. This can be thought of as the dynamics that would arise if the rods have a charge and are subject to an electric field. To model steric repulsions between the rod particles, we introduce a Weeks-Chandler-Andersen (WCA) interaction between the particles and between the particles and wall surface^[40]. We study these effects in a channel with a geometry of 60 nm× 30 nm× 20 nm resolved with a mesh-width 1 nm. The driving field acts in the x-direction. This corresponds to an electric field with the particles having screened net-positive charges.

We perform simulations in different density regimes to investigate the collective stochastic dynamics of the rods and role of hydrodynamic coupling. In particular, we perform simulations of (ⅰ) a dilute suspension of 100 particles and (ⅱ) a dense suspension of 1 050 particles (see Fig. 11). To investigate the role of the hydrodynamic coupling, we perform simulations (ⅰ) with hydrodynamic coupling using our stochastic methods for the SELM and (ⅱ) without hydrodynamic coupling using a standard Langevin dynamics. For a collection of particles within a channel, the hydrodynamic coupling is expected to drive a Poiseuille flow within the channel. This has the consequence that rods near the boundary are transported down the channel at a significantly slower rate. This also has as a consequence shearing of the rods when they are aligned perpendicular to the channel walls. For non-hydrodynamic simulations, these potentially important effects for some applications are neglected.

The profile of the velocity of the rods in the x-direction is shown when varying position along the cross-section of the channel in the z-direction (see Fig. 12). We consider the case when the particles are subject to a constant force field driving them down the channel in the x-direction. We show the velocity V_x(z) of the particles versus their location within the channel relative to the wall in the z-direction. We compare the case without hydrodynamic coupling obtained from Langevin simulations (left) with the case with hydrodynamic coupling obtained from our SELM approach (right). We find that the SELM approach produces velocity distributions in agreement with a Poiseuille flow. The profiles are normalized by the maximum velocity V_x0 across the channel.

We find the characteristic profile of a Poiseuille flow with some slip near the walls. Since the profile reflects the velocity of rods and not directly the fluid, this reflects slip that arises from the motions of the rods near the boundary. In fact, for our SELM formulation, there is allowed slip of the rods relative to the local fluid velocity which is set by the parameter Υ. In contrast, for Langevin simulations, no such spatial dependence in the velocity of the rods is found given the lack of hydrodynamic effects.

We next consider the role of thermal fluctuations and the collective diffusion of particles within the channel. We consider the equilibrium distribution of the translational and rotational configurations of rod. We show these results for the dense regime in Fig. 13. We consider the equilibrium distribution of the particle configurations in the dense regime. We find the particle center-of-mass distributes nearly uniformly across the channel diameter with the density decreasing in a layer near the boundary on a length-scale comparable to the particle size. We find that the angular distribution for pairs of ellipsoidal particles is well-approximated by . We define the angle as θ = arccos(a₁ · a₂/(|a₁||a₂ |)), where a₁ is the axis of the first trimer, and a₂ is the axis of the second trimer. These results are consistent with isotropic arrangements of the particles. The dashed curve for the predicted angle distribution in three dimensions is obtained by associating the trimer directors to points uniformly distributed over the surface of a sphere. We also find near the walls a small depletion in concentration of particles on a scale comparable to the finite size of the particles. This likely arises from the rotational entropy of the particle that would penalize configurations with the particle center-of-mass close to the channel wall.

Concerning the role of thermal fluctuations and hydrodynamic coupling in the collective diffusion of rods within the channel, we further explore the mean-squared-displacement (MSD) of rods within the channel. We consider both the translational and rotational MSD. Given the symmetries of the channel, we consider separately the translational MSD in the x-direction along the axis of the channel and in the z-direction perpendicular to the channel walls. We also consider the MSD when varying the concentration between our dense and dilute regimes (see Fig. 14). We consider the translational MSD and the rotational mean-squared angular (MSA) displacement. We find that the translational MSD is most affected by the increase in particle density and exhibits a noticeable decrease for the denser regime. We find that the MSA is rather insensitive to the density of particles.

We also find the translational MSD has a significant dependence on the direction of the diffusive motion and rod concentration. We find the rotational MSD of the rods does not depend significantly on the rod concentration in the regimes considered. The translational MSD is greater in the diluter regime than the dense regime considered. Interestingly, the translational MSD is significantly smaller when considering the z-direction with a greater difference depending on the rod concentration (see Fig. 14). This appears to arise from the stronger influence of the walls on hydrodynamic flows and trimer motions in the z-direction toward the channel walls which cause greater dissipation and steric restrictions relative to hydrodynamics flows and trimer motions aligned with the channel. We further investigate in the dense regime how the correlations of the fluid and trimers within the channel are influenced by hydrodynamic coupling and wall effects. We investigate how the autocorrelation of the fluid relaxes in response to spontaneous fluctuations in Fig. 15. We show how the velocity autocorrelation varies across the channel width. These results show that the computational methods capture the greater dissipation arising from proximity to the channel wall in the spontaneous hydrodynamics fluctuations. This is reflected in the significantly longer correlation time scale exhibited for velocity fluctuations in the z-slab near the channel center z₀ = 0 versus at distances closer to the wall. The results are normalized from 〈v(z, 0)v(z, 0)〉 = 2.29 nm²/ns² and shown for averages of 10⁶ time steps.

These results establish further that, for the spontenous fluctuations of the fluid velocity, the computational methods capture significant differences in the kinetics arising from local dissipation in the fluid from proximity to the channel walls. The velocity autocorrelation clearly exhibits decay over a significantly longer time scale in the near channel center z₀ = 0, where dissipation arising from shearing motions results in friction in reference to other fluid elements as opposed to near to the wall z₀ = -9.5, where shearing motions result in friction in reference to the stationary no-slip walls.

We characterize the contributions to the kinetics of the trimer motions related these spontaneous fluid velocity fluctuations in the spatial covariance function in Fig. 16. We show how the trimer autocorrelation varies across the channel width. We consider the motion of the center of mass of the trimer X as the location z₀ given by the displacement Δ_z0X = X(τ) -X(0). Here, the τ is 5 time steps and X(0) is in the z-slab corresponding to z₀. We compute the spatial correlation function Φ(z) = 〈Δ_z0XΔ_zX^T〉 and report the e₁^TΦ(z)e₁ component corresponding to the displacements in the direction of the channel^[41]. The results in Fig. 16 show that the computational methods capture the greater dissipation arising from proximity to the channel wall in the spontaneous hydrodynamic fluctuations. This is reflected for the trimers in a significantly stronger spatial correlation in the motions near the channel center z₀ = 0 relative to configurations near the channel wall z₀ = -9.5. This shows the spontaneous fluctuations with the stochastic computational methods capture an effective hydrodynamic screening from the prescence of the walls which augment collective diffusive motions of the trimers. The results are normalized from 〈v(z, 0)v(z, 0)〉 = 2.29 nm²/ns² and averaged from 10⁶ time steps.

The results of these studies show how our stochastic computational methods can capture in the collective dynamics of a large number of interacting particles the roles played by hydrodynamic coupling, thermal fluctuations, and wall effects. We expect our approaches to be useful in the further investigations of transport phenomena associated with particles confined within channel geometries.

6 Conclusions

We have developed stochastic computational methods to take into account confinement effects within channel geometries. Our approach is based on fluctuating hydrodynamics and the SELM for incorporating thermal fluctuations into approximate descriptions of the fluid-structure interactions. We have shown that our approach of imposing an exact fluctuation-dissipation condition to obtain a discretization of the stochastic driving fields combined with an exact projection method to enforce incompressibility is sufficient to ensure results consistent with statistical mechanics. We have shown that our stochastic discretization approach also handles naturally fluctuations in the context of the no-slip boundary conditions. We found that the stochastic computational methods work well to capture both active hydrodynamic responses and the passive diffusive responses when particles are affected by their proximity to the channel walls. We expect our introduced stochastic computational methods to be broadly applicable to applications in which confinement effects play an important role in the dynamics of microstructures subject to hydrodynamic coupling and thermal fluctuations.

Appendix A Staggered grid and discrete transforms

We discuss the Fourier/cosine transform used to compute the inverse Laplacian of the pressure and the discrete Fourier/sine transform used to diagonalize the covariance matrix in Section 3. Our computational domain is decomposed into a collection of N_x× N_y × N_z cells with periodic boundary conditions in the x-and y-directions and Dirichlet or Neumann boundary conditions in the z-direction (see Fig. 1). The locations of data at the staggered lattice sites in the z-direction are illustrated in Fig. A1.

For scalar fields, such as the fluid pressure, the values are stored at the cell center, and periodic boundary conditions are imposed in the x-and y-directions and the Neumann boundary condition in the z-direction. To determine the pressure and handle the incompressibility condition, the discrete Fourier transform (DFT) is used in the x-and y-directions, and the discrete cosine transform (DCT) is used in the z-direction. For the pressure p_k at the lattice site k, the DFTs in the x-and y-directions are computed by

(A1)

(A2)

For the pressure in the z-direction, the DCT is used with the pressure treated as even functions at the nodes and , where the pressure index is j=0, 1, 2, ..., N_z-1. This is computed by

(A3)

We get the full Fourier/cosine transform used for the pressure by combining the transforms above to obtain

(A4)

In practice, we use the FFTW library to compute the transforms^[42]. We use the standard DFT for x-and y-directions and DCT-Ⅱ in the z-direction.

For vector fields, such as the stochastic driving force, component values of F=(F_x, F_y, F_z) are stored on the cell faces of the staggered grid. The numbers of cell faces within the simulation domain and not on the boundary are different in each direction. This gives a slightly different sized lattice for each of the vector components. For F_x, the lattice is N_x× N_y × N_z, for F_y, the lattice is N_x× N_y × N_z, and for F_z, the lattice is N_x× N_y × (N_z-1). For the F_x, F_y, F_z components in the x-, y-, and z-directions, we do a standard DFT similar to the expressions in Eqs. (A1) and (A2). For the F_x and F_y components in the z-direction, we treat the components as odd functions at the index and , where the indices are j=0, 1, 2, ..., (N_z-1). We use the transform

(A5)

(A6)

We treat the F_z component in the z-direction as an odd function at j=-1 and j=N_z-1. We use the transform

(A7)

We get the full transform

(A8)

We summarize in each direction for the vector field components the index ranges a:b inclusive from a to b (see Table A1).

The staggered grid and indexing are illustrated in Figs. 1 and A1. In practice, we use the FFTW library to compute the transforms with the standard DFT in the periodic directions and the discrete sine transform (DST)-Ⅰ and DST-Ⅱ for the z-direction^[42].

Appendix B Protocol for estimating mobility from SELM particle model

We discuss the protocol used for estimating the particle mobility from the shell model and trimer model introduced in Section 5. Our approach estimates mobility by considering the velocity response of a particle to an applied force. This is done using the following steps:

(Ⅰ) A force F_j is applied to each of the control points indexed by j. To probe the response to an applied force or torque in the direction e_k (standard basis vector), we apply forces over the shell model such that either the total force or the torque . This ensures the force or torque sums to unity in the direction k.

(Ⅱ) In the simulation, the control points are allowed to relax for a short duration of time until the system reaches a steady state as indicated by observing the control point velocities. We use the specific criteria on successive time steps that the velocity v_j of each control point satisfies ‖v_jⁿ⁺¹-v_jⁿ‖/‖v_jⁿ‖ < 10^-6.

(Ⅲ) The translational mobility in the direction k is estimated by M_ek = v_avg/‖F‖, where v_avg = with the total number N of control points. The rotational mobility in the direction k is estimated by M_ek = v_ang/‖τ‖, where with the total number N of control points.

We remark that, for the translational mobility estimates, we can achieve the total force sum easily by applying the same uniform force F_j to each control point. However, for the rotational mobility estimates from a discrete collection of control points, the situation is more complicated. We use the convention of applying a symmetric force over the control points using , where , and α_j is a weight. The torque of the discrete collection of points can be expressed using the triple vector product formula . For a configuration of control points with symmetry under reflection about the plane with normal τ, a weight of α_j = 1 can be used. In this case, the first term vanishes, and the second term has . In practice, a discrete collection of control points may not have this symmetry. We setup the control points of our shell model so that the initial configuration has a load force only on a subset of points respecting the planar reflection symmetry. This allows us to use a uniform collection of non-zero weights. In practice, in our simulations, we have found that, to reach the velocity steady-state condition, the total displacements of the control points are negligible. They have typical displacement value less than 10^-2Δx. We take this as an indication, and the particle displacement does not play a significant role in our mobility estimates. As a further test of the rigidity of the particle shell model and our choice of bond stiffness, we have also done further tests by applying a symmetric load force only on a few control points of the surface of the shell model. We find that the particle shell model retains its shape and that the same mobilities are obtained as in the case when we apply load forces over all of the control points of the shell.