The structure and dynamics of confined suspensions of particles of arbitrary shape are of interest in multiple disciplines from biology to engineering. Theoretical studies are often limited by the complexity of long-range particle–particle and particle–wall forces, including many-body fluctuating hydrodynamic interactions. Here, we report a computational study on the diffusion of spherical and cylindrical particles confined in a spherical cavity. We rely on an immersed-boundary general geometry Ewald-like method to capture lubrication and long-range hydrodynamics and include appropriate non-slip conditions at the confining walls. A Chebyshev polynomial approximation is used to satisfy the fluctuation–dissipation theorem for the Brownian suspension. We explore how lubrication, long-range hydrodynamics, particle volume fraction, and shape affect the equilibrium structure and the diffusion of the particles. It is found that once the particle volume fraction is greater than 10%, the particles start to form layered aggregates that greatly influence particle dynamics. Hydrodynamic interactions strongly influence the particle diffusion by inducing spatially dependent short-time diffusion coefficients, stronger wall effects on the particle diffusion toward the walls, and a sub-diffusive regime—caused by crowding—in the long-time particle mobility. The level of asymmetry of the cylindrical particles considered here is enough to induce an orientational order in the layered structure, decreasing the diffusion rate and facilitating a transition to the crowded mobility regime at low particle concentrations. Our results offer fundamental insights into the diffusion and distribution of globular and fibrillar proteins inside cells.

## I. INTRODUCTION

Diffusion under confinement is central to multiple physical, chemical, and biological systems, including colloidal and protein suspensions, devices for particle separation, and transport through membranes. A model system to study the diffusion and structure of highly concentrated particles under confinement could offer insights into the dynamics of crowded macromolecules, such as proteins, inside cells where they typically occupy 20%–40% of the cytoplasm volume.^{1}

Previous studies have shown that crowding between macromolecules affects reaction rates of equilibrium reactions^{2–6} and hinders the diffusion of intra-cellular particles.^{7} *In vivo* experiments, using fluorescence recovery after photo-bleaching (FRAP) techniques, have reported that the apparent diffusion coefficient for green fluorescent proteins (GFP) in *E. coli*’s cytoplasm is about 11 times lower than that in water.^{8–10} Even though a variety of intracellular activities, namely, metabolism, cellular homeostasis, signaling, transcription, translation, and locomotion, strongly depend on diffusion of the intracellular macromolecules, the mechanisms behind the hindered diffusion are not fully understood. A review by Skolnick discussed several factors that could lead to hindrance, including the viscosity of cytoplasm, steric effects, hydrodynamic interactions (HI), and other short-range interactions between particles.^{11} A significant effort that relied on Brownian simulations of 50 types of macromolecules modeled as spheres inside *E. coli*’s cytoplasm was able to reproduce the translational diffusion coefficient of GFP.^{12,13} The authors demonstrated that steric repulsions cannot explain the hindered diffusion and suggested that electrostatic and other short-range interactions are essential variables to consider. However, short-range interactions modeled as van der Waals potentials can always be tuned to match experimental results, and these authors were unable to provide conclusive remarks regarding the effects of short- and long-range HI. In other work, Brownian Dynamics (BD) simulations considering fluctuating HI have been performed in the bulk,^{14} confirming that HI plays a nontrivial role in the hindering of macromolecular diffusion. Unfortunately, the mobility variations induced by confinement were not considered in this work. Confinement, using a network of wall particles that are constrained by a predefined potential, was included in a subsequent study to represent a cell membrane.^{15} In that work, the non-slip boundary condition was not strictly satisfied, allowing the flow to penetrate the membrane and driving a tangential component of the particle average velocity. From the point of view of the role of hydrodynamic interactions on mobility, previous efforts were centered on the study of dynamics for single or many particles immersed in an unconfined viscous fluid^{16–24} and that of particles moving near a wall^{25–30} or confined in a slit,^{26,31–34} in cylindrical geometries,^{26} and in a spherical cavity.^{25,35–37} Recently, a theoretical study examined the behavior of a concentrated colloidal dispersion of spheres confined in a spherical cavity.^{37} That framework relied on a set of hydrodynamic tensors that capture far- and near-field (lubrication) hydrodynamics between particles and walls. The authors studied the structure and diffusion of the system using a Stokesian dynamics (SD) approach^{36} and reported that the confinement, crowding, and HI collectively lead to an anisotropic micro-structure, which then induces position-dependent and anisotropic short- and long-time dynamics. That study was, however, limited to spherical particles.

We have developed an efficient computational framework to perform BD simulations of arbitrarily shaped particles confined in any type of geometry. We use an Immersed-Boundary (IB) method to represent the suspended particles, a parallel finite element general geometry Ewald-like method (pFE-GgEm)^{38} to calculate the confined Green’s functions, and a Chebyshev polynomial approximation to satisfy the fluctuation–dissipation theorem. In this work, we use this methodology to study how steric repulsion, short- and long-range hydrodynamic interactions, confinement, particle volume fraction, and particle shape affect the structure and the diffusion of spherical and cylindrical finite-size particles confined in a spherical cavity. The cylinders are selected to break the three-dimensional symmetry of the particles, a feature that is common in protein structures.

The paper starts with a description of the particle model, the geometry, the methods, and the dimensionless variables to characterize the system. The BD and the IB-pFE-GgEm methods are briefly summarized. We then proceed to present and discuss the results, including the spatial and orientational ordering of the particles and the short- and long-time diffusion behavior. We finish the paper with a summary of the most important findings.

## II. MATHEMATICAL MODEL AND SYSTEM

Let us consider *N* mono-disperse and semi-rigid particles embedded in a viscous fluid of viscosity *η*, confined in an spherical cavity of radius *R*. Under a zero Reynolds number condition, the *N*-body force/torque balance on the particles is

where **F**^{H} is the 6*N* hydrodynamic force/torque vector, **F**^{B} is the Brownian force/torque vector, **F**^{C} is the force/torque vector containing configurational forces, **F**^{EV} represents the force/torque vector due to excluded volume interactions, and **F**^{ext} includes all external forces/torques.

The evolution of the suspended particles, from Eq. (1), is carried out through the grand mobility or resistance tensors that relate the hydrodynamic forces/torques to the translational and rotational velocities of the particles.^{39–41} SD^{42–44} and boundary integral methods (BIM)^{39,45} have been used extensively to solve this “mobility problem.” The regularized Stokeslets,^{46} the accelerated BIM,^{47} and the Immersed Boundary (IB)^{48,49} approaches all provide computational efficiency and simplicity, typically to improve (or avoid) the calculation of the single- and double-layer hydrodynamic potentials of the suspended particles. In particular, the IB method represents the surfaces of the suspended solids as a distribution of discrete force densities that, together with a surface force description and Stokes equations, generate the temporal evolution of the suspended particles. This is the approach that we use in this work.

The surface of each suspended particle is discretized into a set of *N*_{IB} nodes that constitute a mesh, similar to boundary element methods.^{45} On the surface nodes, we define structural spring potentials that maintain particle shape, volume, and surface. The force balance on each of the *N* particles is then translated into the *N*_{IB} surface nodes as follows:

for every node *ν* = 1, …, *N*_{IB}, where $f\nu H$ is the hydrodynamic force, $f\nu B$ is the Brownian force, $f\nu \u2009C$ is the constitutive force, and $f\nu EV$ is the force from all the excluded volume interactions: particle–particle and particle–wall.

Assuming that the probability density for the nodal positions is a continuous density for the Fokker–Planck equation,^{50} an equivalent stochastic differential equation for the motion of the nodes is written as follows:^{51}

where $R=(x1,x2,\u2026,xN\xd7NIB)$ denotes a 3*N* × *N*_{IB} vector containing the spatial coordinates of the nodes, **D** = *k*_{B}*T***M** is the 3*N* × 3*N* diffusion tensor, *k*_{B} is the Boltzmann constant, *T* is the absolute temperature, and **M** is the (3*N* × *N*_{IB}) × (3*N* × *N*_{IB}) mobility tensor. In addition, **U** = **M** · **F** contains the 3*N* fluctuating velocities from the hydrodynamic interactions and **F** is the 3*N* × *N*_{IB} vector that contains the non-HI and non-Brownian forces on the nodes. The divergence of the diffusion tensor $\u2202\u2202R\u22c5D$ is the drift resulting from the configuration-dependent mobility of the confined particles and **dW** is a random vector, the components of which are obtained from a real-valued Gaussian distribution with zero mean and variance *dt*; **dW** is coupled to the diffusion tensor through the fluctuation–dissipation theorem: **D** = **B** · **B**^{T}.

In IB methods, the force distributions at moving solids are discretized as distributions of regularized point-forces. The “smoothing” function for the delta function scales as the distance between the surface nodes that are used to represent the moving particles. Consequently, the structural forces on the particles, $f\nu \u2009C$, define a force density as follows:

where *δ*_{IB}(*a*, **x**) is a smoothing function that is a modified Gaussian (details available in Ref. 52). The regularization parameter *ξ*_{IB} in *δ*_{IB} is related to the characteristic length *h* for the node spacing on the particle surface, i.e., *ξ*_{IB} ∼ *h*^{−1} ∼ *a*^{−1}. The rational behind *δ*_{IB}(*a*, **x**) is to ensure that the regularized force on each node is spread over the length scale of the associated surface elements, thereby preventing fluid from penetrating the particle surface.

In this work, we consider the particles as “semi-rigid,” where each node on the particle is linked to its neighboring nodes by elastic springs with a prescribed large stiffness constant. The nodes are also connected to the particle center-of-mass by an elastic spring to conserve the desired shape. For simplicity, we assumed that each link is a linear spring where the force acting on the point *ν* by the point *μ* is given by

Here, *k* is the spring elastic constant, *q*_{0} is the equilibrium spring size for each specific situation, **x**_{νμ} = **x**_{ν} − **x**_{μ}, and *r*_{νμ} = |**x**_{νμ}|. For each particle, a spring network is formed (see Fig. 1), resulting in an internal nodal force that resists the deformation of the particles. In addition, at each surface node of each particle, a purely repulsive Lennard-Jones (LJ) force is added to account for particle–particle and particle–wall excluded volume interactions. This LJ force comes from two contributions: (1) interactions with nodes from a different particle and (2) interactions with the wall. It is equivalent to the negative gradient of the LJ potential, which is defined by

for *r* ≤ 2^{1/6}*σ* and zero otherwise. Here, *r* is the Euclidean distance between nodes of two different particles or between the node and walls, and *σ* = 2.2*a* is chosen empirically to guarantee that each surface node has an excluded volume of radius *a*. The translational and rotational velocity of particles are calculated by integrating the velocity over the surface mesh of each particle, thereby satisfying the force and torque balance. For completeness, we include a validation of our IB approach in Appendixes A–C: (i) fulfillment of Stokes’ law of particles under confinement, (ii) validation of the fluctuation–dissipation theorem (diffusion of particles), and (iii) consistency of the particle shape as a function of the spring constants.

In what follows, we use dimensionless variables for all results and discussions. We use *a* as the characteristic length scale and the nodal spacing diffusion time, *a*^{2}*ζ*/*k*_{B}*T*, as the characteristic time scale. We set *k*_{B}*T* as the scale for energy and *k*_{B}*T*/*a* for the force. The friction coefficient *ζ* is related to the fluid viscosity *η* and *a* through Stokes’ law: *ζ* = 6*πηa*, and the nodal diffusion coefficient is *D*_{0} = *k*_{B}*T*/*ζ*.

Our simulation method, denoted by **IB-pFE-GgEm**, is an *O*(*N*) algorithm that includes hydrodynamic interactions for confined, large-scale suspensions of finite-size particles of arbitrary shape. Details can be found in our previous work.^{38} Briefly, the *O*(*N*) algorithm consists of three major components: (a) pFE-GgEm routines to calculate the Green’s function (Stokeslet) for any geometry, (b) Fixman’s mid-point algorithm for time integration, and (c) the Chebyshev polynomial approximation for the fluctuation–dissipation theorem. The parallel finite element GgEm (pFE-GgEm) routines are built using open source libraries, thereby facilitating usage of our software. The routines can be downloaded at http://miccomcodes.org as part of the Continuum-Particle Simulation Suite (COPSS) from the mid-west Integrated Center for Computational Materials (MICCoM).

## III. RESULTS

We consider spheres and cylinders, of equal volume, that are suspended in a spherical cavity of radius *R* = 15. The particles’ radius is *r*_{S} = 3 (volume $VHI=4/3\pi rS3$), while the cylinders’ size is determined by *r*_{C} = 2.62 and *h*_{C} = 2*r*_{C} (volume $VHI=\pi rC2hC$). Figure 1 shows the details of the system. According to our semi-rigid particle model, there are two ways to define the particle concentration in a cavity of volume *V*. One uses the hydrodynamic volume fraction, Φ_{HI} = *NV*_{HI}/*V*, and a second one is based on the excluded volume effective size, Φ_{EV} = *NV*_{EV}/*V*. Each surface node has an excluded volume of radius *a*, thus each spherical particle has an excluded volume of $VEV=4/3\pi (rS+a)3$, and each cylindrical particle has an excluded volume of $VEV=\pi (rC+a)2(hC+2a)$. We use the hydrodynamic volume fraction as the relevant scale for particle concentration. In this work, Φ_{HI} = [5%, 10%, 15%, 20%], which is equivalent to Φ_{EV} = [12%, 24%, 36%, 48%]. The lower volume concentration selected here is 5%, which represents a limiting volume fraction between the dilute and finite concentration regimes. We found that all results for dilute systems, starting from infinite dilution, are close to those for the limiting volume fraction of 5%. We start our analysis by exploring the structure of the particles using the particle number density as a function of radial position within the cavity. The number density (the probability that a particle is at a specific location) is calculated by discretizing the spherical cavity into *m* bins (spherical shells) with an even spacing in the radial direction. The shell radius of the *i*th bin is *b*_{i} = (*i* + 0.5)*R*/*m*. The particle number density *n*(*r*_{i}) = ⟨*N*(*r*_{i})/*V*_{i}⟩, where *N*(*r*_{i}) is the number of particles at shell that is located a distance *r*_{i} from the center of the cavity, *V*_{i} is the volume of the shell, and ⟨ ⟩ represents an ensemble average over time. Figure 2 shows the number density of the particles within the cavity as a function of particle concentration. Figure 2 (top) is for spherical particles, while Fig. 2 (bottom) is for the cylindrical ones. At low concentrations, Φ_{HI} = 5%, the number density is uniform throughout the cavity and goes to zero once the particles are in contact with the wall. Note that the maximum density is at *r* ∼ 10.5, which is smaller than wall contact. In our model, particles can never “touch” the wall because of the strong repulsive Lennard-Jones potential. This is common in particle-based and molecular dynamics simulations that consider infinite potentials at zero distances. As the concentration increases, the probability of finding particles near the wall increases, forming a layered structure. At moderate concentrations, Φ_{HI} = 10%, 15%, particles start to form the first layer next to the wall, and inner particles prefer to stay near the center of the cavity where the steric effects with the particles in the first layer are the weakest. At higher concentrations, Φ_{HI} = 20%, particles form a second layer, since there is not enough space at the center to accommodate them. Thus, the layered structure becomes more pronounced, and the layer separation is determined by the particle size. In the figure, we include the number density of a “bead” system where particles are represented as spheres with an excluded volume with radius *a*. The spherical cavity has a radius of *R* = 10*a*. In this system, particles interact only through far-field hydrodynamics and steric repulsions. However, the layered structure is also observed, becoming a characteristic of highly concentrated confined systems.

Cylindrical particles exhibit an orientational distribution within the cavity. Similar to liquid crystalline systems,^{53,54} we define an orientational order parameter $\lambda =12\u27e83\u2061cos2\u2061\theta \u22121\u27e9$, where $cos\u2061\theta =m\u22c5n\Vert m\Vert \u22c5\Vert n\Vert $, **m** is the vector parallel to the cylinder’s centerline, and **n** is the vector pointing from the cavity center to the cylinder’s center-of-mass. A random/disordered structure is characterized by *λ* = 0, whereas for ordered morphologies, *λ* = 1 when all cylinders are aligned parallel to the radial direction of the spherical cavity (radial phase), and *λ* = −1/2 when all cylinders are aligned perpendicular to the radial direction of the spherical cavity (concentric phase). Figure 3 shows the orientational order parameter of cylindrical particles within the cavity as a function of particle concentration. A major result is that for all concentrations, the cylinders are in a disordered state at the center of the cavity and oriented concentrically near the wall. For Φ_{HI} = 5%, the disordered morphology spans all locations in the cavity, whereas for Φ_{HI} = 15%, two concentric regions are observed, separated by a layer where the cylinders are oriented following a ⟨*θ*⟩ = 35*°* or 145*°*. Interestingly, the radial–centripetal ordering, *λ* = 1, is never observed; we believe that this is due to the small aspect ratio of the cylindrical particles.

We now proceed to analyze the short- and long-time diffusive behavior in an attempt to delineate the consequences of long- and short-range HI on the dynamics of the particles. The short-time diffusion coefficient and the mean squared displacement (MSD) of the particles are used to quantify these effects. The particles are suspended in a viscous fluid under zero Reynolds number conditions; they interact with other particles and walls through the HI. Recall that under these conditions, momentum transport is infinitely fast.^{55–58} Previous studies of confined suspensions have shown that there are multiple factors that originate from HI that should affect the diffusion (mobility) of the particles: (i) the reduction of the particle mobility due to confinement—particles in the bulk diffuse faster than in confined geometries, (ii) the space dependent mobility due to the non-slip conditions at the walls—particle diffusion is zero at the walls, and (iii) the decrease of particle mobility as the concentration increases—there is an interplay between lubrication and long-range HI that becomes important as the average inter-particle distance decreases.^{29,34,37,59} We seek to quantify and determine the consequences of such factors in a cavity enclosure using our model.

We start by measuring the particles’ radial and tangential short-time diffusivities within the cavity as a function of particle concentration. These transport coefficients are calculated from the relation between MSD and time from Stokes–Einstein following a directional decomposition,^{60} i.e.,

for *t* → 0 and where Δ**x** = **x( t** +

**−**

*dt*)**x(**, the radial displacement Δ

*t*)**x**

_{R}= Δ

**x**·

**x**/|

**x**|, the tangential displacement Δ

**x**

_{T}= Δ

**x**− Δ

**x**

_{R},

*D*

_{R}(

*r*

_{i}) and

*D*

_{T}(

*r*

_{i}) are the instantaneous radial and tangential short-time diffusivities at a distance

*r*

_{i}from the center of the cavity, and

*dt*is an infinitesimal time interval.

*D*

_{R}(

*r*

_{i}) and

*D*

_{T}(

*r*

_{i}) are blocked averaged during a typical simulation at each shell and then averaged over independent simulations. Figure 4 summarizes the diffusion coefficients within the cavity for spherical particles as a function of concentration. In the figure, the diffusivities are normalized by the bulk value, which is defined as single particle diffusivity at infinite dilution.

We find that the confinement hinders the particle diffusion in both directions. The highest value for the particle diffusivity at Φ_{HI} = 5% is around 60% of the bulk value. The lower particle mobility is directly related to the long-range character of the HI and the non-slip conditions at the walls. As the particle concentration increases, lubrication forces begin to dominate and particle diffusion decreases monotonically. In addition, the short-time coefficients are not constant within the cavity, showing sudden decrease as the particles approach the wall. For Φ_{HI} = 10%, we calculated the averaged “inner” coefficients, which are represented by the orange dashed lines in Fig. 4. Importantly, the decrease in particle mobility at the walls has a stronger effect on the radial diffusion coefficient, indicated by (i) the radial particle mobility decreases at *r* = 7, when compared with the tangential diffusivity that “feels” the presence of the walls at *r* = 9, and (ii) the radial diffusivity decreases by 0.2 from *r*(*a*) = 7 to the wall contact, whereas the tangential diffusivity only decreases by about 0.12. For a spherical particle near the walls, it is easier to diffuse concentrically than radially. Finally, note that for Φ_{HI} = 20%, the coefficients adopt a non monotonic character, which is correlated with the layered structure of particle density in the cavity. The diffusivity data for *r*(*a*) < 2.3 are missing for the case of *ϕ*_{HI} = 20% because the short-time diffusion coefficient at a given location is measured on particles that appear at that location; in other words, the diffusion coefficients are only measured in regions where there is a finite concentration of particles [see Fig. 2 (top)]. Similar observations have also been reported by Aponte-Rivera *et al.* using a Stokesian dynamics (SD) approach.^{37}

In Fig. 5, we include the short-time diffusion coefficients of spheres with an excluded volume confined to a spherical cavity, with and without long-range HI (no lubrication forces). Our intention is to change the level of the particle description to isolate the HI contributions. Free draining point-particles do not undergo a space dependent diffusion and the concentration perturbs the diffusion rate only weakly. For HI point particles, the mobility decreases monotonically as the particles approach the walls, and the diffusion rate has a stronger dependence on concentration. Consequently, long-range HI and the zero mobility at the walls are responsible for the non-uniform particle diffusion inside the cavity. Lubrication, on the other hand, imposes a directionality on the short-time mobility—radial diffusion is different than concentric diffusion near the walls—and correlates the particle diffusion with the layered density at high volume fractions. The observed behavior is consistent with the experimental findings.^{61,62}

It is of interest to compare the short-time diffusion coefficient between spherical and cylindrical particles and to validate the effects of the level of confinement. Figure 6 includes the short-time diffusion coefficients for cylindrical particles (left) in a cavity with *R* = 15 and for spherical particles (right) in a cavity with *R* = 30. In the figure, Φ_{HI} = 5% and the results for spherical particles in a cavity with *R* = 15 are included for reference. Similar to the spherical particles, the cylinders exhibit position-dependent diffusion coefficients, and radial diffusion is affected strongly by the presence of the walls when compared to tangential diffusion. Interestingly, the shape of the cylindrical particles has an important effect on the rate of diffusion. Recall that the cylinders and spheres have the same volume and that the diffusion coefficients are normalized by the bulk diffusivity of spherical particles. Consequently, the short-time diffusion of particles of equal volume is decreased when the geometrical symmetry is broken. Finally, decreasing the level of confinement does not change the qualitative behavior of the short-time mobility but, as the confinement decreases, the inner diffusion coefficients approach the bulk value.

Finally, we examine the long time diffusive behavior. We use a generalized Stokes–Einstein relation where the MSD is linearly correlated with a mobility coefficient following a power law by

where **R** is the 3*N* particle coordinate vector, *M* is the particle mobility coefficient, and *α* is the power law exponent that characterizes the type of particle transport. For isotropic diffusion, if *α* = 1, the particles are diffusive and *M* = 6*D*, where *D* is the particle long-time diffusion coefficient. If *α* ≠ 1, the transport of Brownian particles is said to be in the anomalous diffusion regime: *α* < 1 is sub-diffusive, while *α* > 1 is super-diffusive. For confined systems, as *t* → *∞*, the walls impose long-time restrictions on the motion along the confining direction. Therefore, for particles confined in a spherical cavity, the MSD should reach a plateau on a time scale corresponding to the particle diffusion time over the cavity size.

Figure 7 shows the MSD of spherical (top) and cylindrical (bottom) particles in a cavity of *R* = 15 as a function of particle concentration. The time scale in Fig. 7 is in *a* diffusion time unit, *a*^{2}*ζ*/*k*_{B}*T*, where *ζ* = 6*πηa*. In the figures, we include the results for Φ_{HI} = 15% of the other particle shape for comparison purposes. The MSDs are collected from ten independent simulations with different random seeds; the particles were able to diffuse more than 300 *a* diffusion times. As expected, the MSD exhibits diffusive behavior as *t* → 0 and a plateau when MSD ∼ 200 ∼ *R*^{2} at *τ*_{R}, where *τ*_{R} is characteristic particle diffusion time over the cavity. At low concentrations, the diffusive behavior, for spheres and cylinders, spans from *t* = 0 to *t* = *τ*_{R} (for 5% spheres, *τ*_{R} ∼ 300 in *a*^{2}*ζ*/*k*_{B}*T* units). As the particle concentration is increased, the diffusion rate decreases; it is correlated with the short-time behavior, as indicated by the monotonic shift of the MSDs in Fig. 7. For spherical particles at Φ_{HI} ≥ 15%, there is a clear sub-diffusive regime over more than two decades that increases as the concentration is increased. Interestingly, the sub-diffusive regime for cylindrical particles starts at Φ_{HI} = 10%. Note that the MSD for cylindrical particles at Φ_{HI} = 15% is almost equal to the MSD for spherical particles at Φ_{HI} = 20%. In the inset, we have included the MSD of spheres with excluded volume but no lubrication to show how the sub-diffusive regime is never observed when only long-range HI is included, thereby suggesting that the sub-diffusive regime is a lubrication effect. It is then natural to attribute this anomalous diffusion to crowding driven by short-range HI. This regime is then characterized by the diffusion-to-crowding transition *τ*_{DtC} time and by the power law exponent *α*. In this regime, there are two major features: *τ*_{DtC} and *α* decrease as the particle concentration increases. In particular, *τ*_{DtC} = 20 and *α* = 0.85 for cylinders at Φ_{HI} = 10%, while *τ*_{DtC} = 5 and *α* = 0.65 for cylinders at Φ_{HI} = 15%. After the highly concentrated systems enter the crowding regime and given the fact that there are three dimensional restrictions in their motion as *t* approaches *τ*_{R}, the particles transition from the slow rate sub-diffusive behavior to the plateau. In Appendix C, we show that our mean square displacements for spheres agree with the isotropic mean square displacement analysis reported in the literature^{36,37} both qualitatively and quantitatively. Past work^{37} reported on the anisotropic behavior of the long-time MSD where both sub-diffusive and subsequent super-diffusive regimes appear at intermediate time scales in the radial direction for highly concentrated systems, while the super-diffusive regime is absent in the tangential direction. In view of that past contribution, our work focuses on the combined effects of different factors on the isotropic long-time MSD instead and provides a first step toward understanding the role of the particle shape on the underlying dynamics.

## IV. CONCLUSIONS

We have used an immersed boundary approach to study the structure and dynamics of suspended spherical and cylindrical particles confined in a spherical cavity.

At low concentrations, the particle number density distribution is uniform in the interior of the cavity. As the concentration increases, a layered structure appears. Cylindrical particles exhibit a random orientation at low concentrations and at the center of the cavity for all concentrations. Excluded volume interactions at the wall force the cylinders to orient concentrically. Interestingly, at high concentrations, the layered morphology of the cylinders is correlated with concentric layers that are separated by a layer of cylinders oriented at an average angle of 35*°* with respect to the radial direction. Cylinders are never found forming radial morphologies within the spherical cavity.

We used multiple hydrodynamical descriptions of the suspended particles to determine the specific influence of hydrodynamic interactions during the dynamic of the particles (diffusion/mobility): (i) free draining point-particles, (ii) long-range hydrodynamic “beads,” and (iii) finite size particles considering lubrication and long-range HI. We found that long-range HI leads to a position-dependent diffusion of the particles inside the cavity; the particles diffuse faster near the center of the cavity and slower near the walls. The HI also decrease the mobility of the suspended particles when compared with their diffusion in the bulk. The increase in particle concentration also results in a decrease in the particles’ diffusion coefficients; this effect is observed for free-draining and HI particles. However, the concentration decrease in the diffusion rate is stronger when HI are considered. Lubrication forces, or short-range HI, influence the dynamics of highly concentrated suspensions; they generate a direction dependent diffusion, where particles diffuse at a lower rate when moving toward the walls than when moving parallel to the walls. The non-slip conditions at the walls, i.e., zero particle mobility, work synergistically with lubrication forces, resulting in an stronger wall dependence of the diffusion coefficients in the radial direction.

Regarding the long–time dynamics, lubrication gives rise to a sub-diffusive regime at high particle concentrations. The sub-diffusive regime, characterized by the diffusion-to-crowding transition time and the mobility power law exponent, becomes more prominent as the concentration increases.

Introducing cylindrical particles has two major consequences: (i) cylindrical particles have lower short-time diffusion coefficients and (ii) the crowding regime is observed at lower concentrations compared with spheres of equal volume. These observations suggest that the shape of bio-molecules, particles, and polymers could determine their mobility and diffusion inside cells and tissues.

## AUTHOR’S CONTRIBUTIONS

J.L. and X.J. contributed equally to this work.

## ACKNOWLEDGMENTS

The development of fast computational codes for simulations of nanoparticles interacting through hydrodynamic and polarization interactions was supported by the Department of Energy, Basic Energy Sciences, Division of Materials Research, through the MICCoM center. The calculations on particle segregation and transport presented in this work were supported by the Department of Energy, Basic Energy Sciences, Division of Materials Research, through the AMEWS center.

### APPENDIX A: SUSPENDED BROWNIAN SPHERES WITH THE IB-pFE-GgEm

There are three important validations that are in order to verify the immersed boundary approach that we used in this work. We start by verifying the fulfillment of Stokes’ law by measuring the sedimentation velocity of a spherical particle that is confined between two parallel walls. Analytical values for this velocity are extracted from previous studies.^{22,32,33,44,63–65} For the specific IB-pFE-GgEm calculation, the particle radius *r*_{S} = 5 and the wall distance is *H* = 15. The periodic boundary condition (PBC) is enforced in the unconfined directions, which are set to a length of 200 (≫*H*) to avoid the influence of the particle periodic images. The surface of the spherical particle is discretized using 119 surface nodes, resulting in a nodal separation between *h*_{min} = 1.334 08 and *h*_{max} = 1.905 39 and a smoothing parameter *ξ*_{IB} = 1/0.76*h*_{min}. For the calculation, we used a GgEm parameter *α*_{GgEm} = 0.2 and a mesh resolution with a spacing of $1/2\alpha $. The resulting mesh was 60 × 60 × 6 HEX20 elements with 324 886 degrees of freedom. The particle is initially located at (0, 0, *d*) between the two parallel walls, and it moves parallel to the walls under a sedimenting force with a body force density of (0, 0, 1). The particle’s sedimenting velocity (*U*_{∥}) is calculated by averaging particle’s velocities over 100 time steps. Figure 8 (top) shows the sedimenting velocity, normalized by Stokes’ velocity *U*_{∥}/*U*_{0}, as a function of normalized location [(*d* − *R*)/*H*]. According to the results, the IB-pFE-GgEm provides an excellent agreement with analytical values.

In addition to Stokes’ law, it is important to verify that our combined approach, between the IB-pFE-GgEm, Fixman’s mid-point algorithm, and the Chebyshev polynomial approximation, is satisfying the fluctuation–dissipation theorem. A sphere that is confined in a spherical cavity and between parallel walls offers a scenario to validate the connection between the diffusion and the fluctuating tensors and the proper calculation of the mobility gradients inside the confined geometries. First, we used three different particle discretizations, using 20, 40, and 60 nodes, for the sphere of size *r*_{S} = 3 that is confined in a spherical cavity of size *R* = 15. Figure 8 (center) shows the comparison between the short-time diffusion coefficient of the sphere computed through our IB-pFE-GgEm algorithm and SD algorithm previously reported in the literature.^{37} These results suggest that even with poor surface descriptions, as long as the IB parameter is appropriately chosen, the IB-pFE-GgEm follows the correct fluctuating short-time behavior. Finally, the long-time diffusion for a sphere confined in a slit as a function of the separation between the walls is shown in Fig. 8 (bottom). Analytical and numerical solutions for this coefficient had been well assessed in the literature.^{68–71} For the IB-pFE-GgEm, the particle is initially located at the mid-plane between two parallel walls. During the particle’s Brownian movement, its motion is restricted to the plane of symmetry. The surface of the spherical particle is discretized using 20 surface nodes (*h*_{min} = 2.198 26, *h*_{max} = 2.526 08, and *ξ*_{IB} = 1/0.35*h*_{min}). The slit mesh resulted in 60 × 60 × 4 HEX20 elements with 228 872 degrees of freedom. To calculate the error bars of the MSD, five independent simulations for each confinement ratio are performed with a constant time step of 0.002.^{70} The diffusion coefficients from the IB-pFE-GgEm shows an excellent agreement with analytical and numerical results.

### APPENDIX B: RIGIDITY OF THE SUSPENDED PARTICLE

In our immersed boundary description of the finite-size particles, each surface node is linked to its neighboring nodes and the center-of-mass through linear springs with a prescribed stiffness constant *k*. This parameter controls the stiffness (shape) of particles. If the springs are too weak, particles are deformable and special care must be taken to forbid fluid penetration. On the other hand, if the springs are too stiff, the forces acting on surface nodes will be too large, requiring very small time steps to ensure numerical stability. We performed simulations of spherical particles that are confined in a spherical cavity at high concentrations varying the particle stiffness. We measured the particles’ moment of inertia for each particle *ν* as follows:

where *I*_{ν}(*t*) is the moment of inertia of particle *ν* at time *t* and *r*_{j}(*t*) is the distance between node *j* on surface and the particle center-of-mass. Rigid spheres will have an equal and constant moment of inertia. On the other hand, our semi-rigid particles will show a variation in *I*_{ν}(*t*). Figure 9 shows the averaged standard deviation of moment of inertia as a function of time for particles with different *k* for Φ_{HI} = 20%. The standard deviation is calculated by

where *I*_{ν}(0) is the moment of inertia of particle *ν* at time 0, which corresponds to a perfect sphere. As one can observe from the figure, in the case that *k* = 10, the moment of inertia jumps from 0.05 to 0.06 at *t* ≈ 10, indicating that *k* = 10 is not stiff enough to maintain particle shape. As the stiffness is increased, the shape variations decrease. In this paper, we used a *k* = 200 for all simulations.

### APPENDIX C: COMPARISON OF MSD WITH LITERATURE RESULTS

The equilibrium structure and diffusion in concentrated hydrodynamically interacting spherical particles confined in a spherical cavity were studied using SD by Aponte-Rivera *et al.*^{37} There are differences between methods in our and their work such as particles models and HI models. However, one can draw both quantitative and qualitative agreements between results from two methods. Here, we present the comparison of isotropic MSD at different volume fractions [Fig. 15(a) in Ref. 37 and the top panel in Fig. 7 in this work]. Qualitatively, a short-time diffusive regime and a long-time plateau for systems at all concentrations are observed in both works. At intermediate times, a sub-diffusive region emerges for highly concentrated systems only. To make a semi-quantitative comparison, we convert the MSD data in Ref. 37 to the units used in our work at the same level of excluded volume crowding. Specifically, the length and time reported by Aponte-Rivera *et al.* are normalized by the radius of a sphere with radius *r*_{s} and the characteristic diffusion time *t*_{0} of a sphere with radius *r*_{s}, respectively ($t0=rs2/D0$, where *D*_{0} = *k*_{B}*T*/6*πηr*_{s}). By multiplying the time, MSD, and volume fraction in the their work by $rs3$, $rs2$, and $(rs/rs+a)3$ (*r*_{s} = 3*a*), respectively, we can compare the isotropic MSD directly. Note that the hydrodynamic crowding, Φ_{HI}, and the excluded volume crowding, Φ_{EV}, are equivalent in the work of Aponte-Rivera *et al.* In contrast, in our work, $\Phi EV/\Phi HI=(rs+a)3/rs3$. Thus, the comparison is conducted at the same level of excluded volume crowding but not the same hydrodynamic crowding. For example, the highest concentration in Ref. 37 is 40%, which is comparable to the highest excluded volume crowding in our work, Φ_{EV} = 47% (Φ_{HI} = 20%). For this highest concentrated system, in the work of Aponte-Rivera *et al.* (after conversion to our units, the converted value is reported in parentheses in the following), the short-time diffusive regime lasts from t ≈ 10^{−3} (10^{−2}) to ≈10^{−1} (10^{0}), while the MSD increases from 10^{−3} (10^{−2}) to 10^{−1} (10^{0}); the intermediate sub-diffusive regime lasts from t ≈ 10^{−1} (10^{0}) to t ≈ 10^{1} (10^{2}), while the MSD increases from 10^{−1} (10^{0}) to 10^{0} (10^{1}); and the MSD reaches a plateau at 10^{1} (10^{2}) at t ≈ 10^{2} (10^{3}). Good agreement is found between our MSD data and those reported in Ref. 37.