Tagged-particle motion in quasi-confined colloidal hard-sphere liquids

In this paper we investigate the dynamics of a tracer (tagged) particle with diameter σ_s and bare diffusion coefficient ${D}_{0}^{\left(s\right)}$ suspended in a 3D quasi-confined colloidal liquid of identical hard spheres with diameter σ and bare diffusion coefficient D₀. All particles evolve according to Brownian dynamics, where hydrodynamic interactions are ignored and are restricted to a quasi-confined environment by applying periodic boundary conditions along the confining direction, which we refer to as the z-direction. Thus, the transversal coordinates are folded back to −L/2 ⩽ z ⩽ L/2. The in-plane coordinates are denoted by $\overrightarrow {r}=\left(x,y\right)$ and are infinitely extended in the thermodynamic limit, N → ∞, A → ∞, where N is the number of particles and A is the area in the dimensions perpendicular to the confining direction. We introduce the area density n₀ = N/A and the volume density n = n₀/L related to the packing fraction φ = nπσ³/6 of the system. All coordinates and quantities of the tracer particle are labeled with a superscript s.

In this section, we convert the relevant MCT equations from the more general class of layered colloidal liquids [35] to quasi-confined fluids. For this purpose we basically follow reference [33] in which the corresponding collective equations have been derived. Since most of the steps of the derivation are identical if collective quantities are replaced by the corresponding single-particle quantities, we just provide a short summary of the main equations, focusing on the relevant differences to the collective case.

The most relevant correlation function is the SISF, also called the incoherent intermediate scattering function,

$$\begin{equation}{S}_{\mu }^{\left(s\right)}\left(q,t\right)=\langle {\rho }_{\mu }^{\left(s\right)}{\left( \overrightarrow {q},t\right)}^{{\ast}}{\rho }_{\mu }^{\left(s\right)}\left( \overrightarrow {q}\right)\rangle ,\end{equation} \tag{ 1 }$$

given in terms of the fluctuating density modes

$$\begin{equation}{\rho }_{\mu }^{\left(s\right)}\left( \overrightarrow {q}\right)=\mathrm{exp}\left[\mathrm{i}{Q}_{\mu }{z}_{s}\right]{\mathrm{e}}^{\mathrm{i} \overrightarrow {q}\cdot { \overrightarrow {r}}_{s}},\end{equation} \tag{ 2 }$$

of the tagged particle. The initial value ${S}_{\mu }^{\left(s\right)}\left(q,t=0\right)=1$ is independent of q and μ. The two-dimensional wavevector $\overrightarrow {q}=\left({q}_{x},{q}_{y}\right)$ in the lateral direction becomes continuous in the thermodynamic limit, whereas the wavenumber ${Q}_{\mu }=2\pi \mu /L,\mu \in \mathbb{Z}$ is discrete. For L → ∞ rotational invariance with the 3D wavevector $\overrightarrow {k}=\left( \overrightarrow {q},{Q}_{\mu }\right)$ is restored and the bulk limit is reached as it has been demonstrated for the general class of layered fluids [39]. Due to periodic boundary conditions, the accessible z-coordinate positions do not differ between the tracer and the host-liquid particles, even if the tracer particle varies in size. This is in contrast to the case of a tracer particle confined between two parallel walls, where the effective confinement length L_s perpendicular to the walls for the tracer particle does not have to coincide with the accessible slit width L for the host-liquid particles [34]. Therefore, for the tagged particle within the slit geometry an additional plane-wave basis $\mathrm{exp}\left(-\mathrm{i}{Q}_{\mu }^{\left(s\right)}z\right)$, with ${Q}_{\mu }^{\left(s\right)}=2\pi \mu /{L}_{s},\mu \in \mathbb{Z}$ is required. Quasi-confinement, however, simplifies the situation and no additional plane-wave basis is necessary.

The exact equations of motion (e.o.m.) for ${S}_{\mu }^{\left(s\right)}\left(q,t\right)$ can be derived using the Zwanzig projection operator formalism [11] introducing an effective self-memory kernel ${M}_{\mu }^{\left(s\right)}\left(q,t\right)$ analogous to the collective case [33]

$$\begin{equation}{\dot {S}}_{\mu }^{\left(s\right)}\left(q,t\right)+{D}_{\mu }^{\left(s\right)}\left(q\right){S}_{\mu }^{\left(s\right)}\left(q,t\right)+{D}_{\mu }^{\left(s\right)}\left(q\right){\int }_{0}^{t}{M}_{\mu }^{\left(s\right)}\left(q,t-{t}^{\prime }\right){\dot {S}}_{\mu }^{\left(s\right)}\left(q,{t}^{\prime }\right)\mathrm{d}{t}^{\prime }=0,\end{equation} \tag{ 3 }$$

with initial decay ${D}_{\mu }^{\left(s\right)}\left(q\right)={D}_{0}^{\left(s\right)}\left({q}^{2}+{Q}_{\mu }^{2}\right)$ related to the tagged-particle bare diffusion coefficient ${D}_{0}^{\left(s\right)}$. Due to the quasi-confined geometry the relaxation naturally splits into channels parallel and perpendicular to the confining direction. Following the same procedure as for the collective case, the effective kernel is related to the matrix-valued irreducible memory kernel ${\boldsymbol{\mathcal{M}}}_{\mu }^{\left(s\right)}\left(q,t\right)$ with diagonal components ${\mathcal{M}}_{\mu }^{{\Vert}{\Vert},\left(s\right)}\left(q,t\right)$ and ${\mathcal{M}}_{\mu }^{\perp \perp ,\left(s\right)}\left(q,t\right)$ and the symmetric off-diagonal element ${\mathcal{M}}_{\mu }^{{\Vert}\perp ,\left(s\right)}\left(q,t\right)$ by the relation

$$\begin{align}\hfill {D}_{\mu }^{\left(s\right)}\left(q\right){M}_{\mu }^{\left(s\right)}\left(q,t\right)& +{D}_{0}^{\left(s\right)2}{\int }_{0}^{t}{M}_{\mu }^{\left(s\right)}\left(q,t-{t}^{\prime }\right){\alpha }_{\mu }^{\left(s\right)}\left(q,{t}^{\prime }\right)\mathrm{d}{t}^{\prime }={D}_{0}^{\left(s\right)}{\beta }_{\mu }^{\left(s\right)}\left(q,t\right)\hfill \\ \hfill & +{D}_{0}^{\left(s\right)2}{\int }_{0}^{t}{\mathcal{M}}_{\mu }^{{\Vert}{\Vert},\left(s\right)}\left(q,t-{t}^{\prime }\right){\mathcal{M}}_{\mu }^{\perp \perp ,\left(s\right)}\left(q,{t}^{\prime }\right)\mathrm{d}{t}^{\prime }\hfill \\ \hfill & -{D}_{0}^{\left(s\right)2}{\int }_{0}^{t}{\mathcal{M}}_{\mu }^{{\Vert}\perp ,\left(s\right)}\left(q,t-{t}^{\prime }\right){\mathcal{M}}_{\mu }^{{\Vert}\perp ,\left(s\right)}\left(q,{t}^{\prime }\right)\mathrm{d}{t}^{\prime }.\hfill \end{align} \tag{ 4 }$$

Therein, the coefficients

$$\begin{equation}{\alpha }_{\mu }^{\left(s\right)}\left(q,t\right)={Q}_{\mu }^{2}{\mathcal{M}}_{\mu }^{{\Vert}{\Vert},\left(s\right)}\left(q,t\right)+{q}^{2}{\mathcal{M}}_{\mu }^{\perp \perp ,\left(s\right)}\left(q,t\right)-2q{Q}_{\mu }{\mathcal{M}}_{\mu }^{{\Vert}\perp ,\left(s\right)}\left(q,t\right),\end{equation} \tag{ 5 }$$

and

$$\begin{equation}{\beta }_{\mu }^{\left(s\right)}\left(q,t\right)=\frac{{q}^{2}}{{q}^{2}+{Q}_{\mu }^{2}}{\mathcal{M}}_{\mu }^{{\Vert}{\Vert},\left(s\right)}\left(q,t\right)+\frac{{Q}_{\mu }^{2}}{{q}^{2}+{Q}_{\mu }^{2}}{\mathcal{M}}_{\mu }^{\perp \perp ,\left(s\right)}\left(q,t\right)+\frac{2q{Q}_{\mu }}{{q}^{2}+{Q}_{\mu }^{2}}{\mathcal{M}}_{\mu }^{{\Vert}\perp ,\left(s\right)}\left(q,t\right),\end{equation} \tag{ 6 }$$

are completely determined by the matrix entries of the irreducible memory kernel ${\boldsymbol{\mathcal{M}}}_{\mu }^{\left(s\right)}\left(q,t\right)$.

Finally, MCT approximates the irreducible memory kernel as a bilinear functional of the ISF and the SISF, where the coupling between different modes and wavenumbers is described by vertices. These are entirely determined by the direct correlation function ${c}_{\mu }^{\left(s\right)}\left(q\right)$ of the tagged particle, which is connected to the static structure factor ${S}_{\mu }^{\left(s\right)}\left(q\right)$ by a generalized Ornstein–Zernike relation [32]. Numerical efficiency can be remarkably enhanced by introducing bipolar coordinates [33]

$$\begin{align}\hfill {\mathcal{M}}_{\mu }^{\alpha \beta ,\left(s\right)}\left(q,t\right)& =\frac{n}{4{L}^{3}{\pi }^{2}}{\int }_{0}^{\infty }{\int }_{\vert q-{q}_{1}\vert }^{q+{q}_{1}}\sum\limits _{\begin{subarray}{c}{\mu }_{1}=-\infty \\ {\mu }_{2}=\mu -{\mu }_{1}\end{subarray}}^{\infty }{S}_{{\mu }_{1}}\left({q}_{1},t\right){S}_{{\mu }_{2}}^{\left(s\right)}\left({q}_{2},t\right)\hfill \\ \hfill & \quad {\times}\frac{4{q}_{1}{q}_{2}}{\sqrt{4{q}^{2}{q}_{1}^{2}-{\left({q}^{2}+{q}_{1}^{2}-{q}_{2}^{2}\right)}^{2}}}\left[{b}^{\alpha }\left(\frac{{q}^{2}+{q}_{1}^{2}-{q}_{2}^{2}}{2q},{Q}_{{\mu }_{1}}\right){c}_{{\mu }_{1}}^{\left(s\right)}\left({q}_{1}\right)\right]\hfill \\ \hfill & \quad {\times}\left[{b}^{\beta }\left(\frac{{q}^{2}+{q}_{1}^{2}-{q}_{2}^{2}}{2q},{Q}_{{\mu }_{1}}\right){c}_{{\mu }_{1}}^{\left(s\right)}\left({q}_{1}\right)\right]\mathrm{d}{q}_{1}\enspace \mathrm{d}{q}_{2},\hfill \end{align} \tag{ 7 }$$

where b^α (x, z) = xδ_α∥ + zδ_α⊥ is the selector with respect to the channel indices α, β ∈ {∥, ⊥}. The indices μ₁ and μ₂ run over all integer values satisfying the selection rule μ = μ₁ + μ₂, which is valid due to translational invariance along the confining direction. Since not only the SISF, but also the ISF enters the equations and all tagged-particle quantities within MCT are derived quantities, the collective dynamics have to be solved first. The algorithm for the tagged-particle dynamics is equivalent to the collective case. Solving the coupled equations for the effective (self-)memory kernel and the (S)ISF is not straightforward. A detailed description of the algorithm using a decimation scheme can be found in [33, 40].

To distinguish ergodic liquid states from arrested localized states it is useful to introduce the nonergodicity parameter or glass-form factor as the long-time limit of the SISF,

$$\begin{equation}{F}_{\mu }^{\left(s\right)}\left(q\right){:=}\underset{t\to \infty }{\mathrm{lim}}\enspace {S}_{\mu }^{\left(s\right)}\left(q,t\right).\end{equation} \tag{ 8 }$$

Liquid states then correspond to vanishing nonergodicity parameters, whereas ${F}_{\mu }^{\left(s\right)}\left(q\right)$ remains finite for glassy states. The nonergodicity parameters are directly accessible in simulations or experiments, and within the framework of MCT ${F}_{\mu }^{\left(s\right)}\left(q\right)$ can be determined without solving the full time-dependent equations explicitly [33, 34].

We now elaborate only the lateral dynamics of the quasi-confined liquid from a 'top-view' perspective and investigate the MSD and the VACF. Even though no confinement is present in the parallel direction, the dynamics of this two-dimensional projection are still sensitive to the confining length. To derive explicit equations for the MSD and the VACF we investigate the tagged-particle dynamics in the long-wavelength limit. The small q series of the SISF

$$\begin{equation}{S}_{\mu }^{\left(s\right)}\left(q,t\right)={P}_{\mu }\left(t\right)-{q}^{2}\delta {r}_{\mu }^{2}\left(t\right)/4+\mathcal{O}\left({q}^{4}\right),\end{equation} \tag{ 9 }$$

together with the out-of-plane SISF

$$\begin{equation}{P}_{\mu }\left(t\right)=\left\langle \mathrm{exp}\left[-\mathrm{i}{Q}_{\mu }\left(z\left(t\right)-z\left(0\right)\right)\right]\right\rangle ,\end{equation} \tag{ 10 }$$

can be written in terms of a generalized MSD

$$\begin{equation}\delta {r}_{\mu }^{2}\left(t\right){:=}\left\langle \mathrm{exp}\left[-\mathrm{i}{Q}_{\mu }\left(z\left(t\right)-z\left(0\right)\right)\right]{\left[ \overrightarrow {r}\left(t\right)- \overrightarrow {r}\left(0\right)\right]}^{2}\right\rangle .\end{equation} \tag{ 11 }$$

P_μ (t) does not depend on the in-plane wavevector $\overrightarrow {q}$, and probes the dynamics of the tracer only along the confinement direction. Inserting the expansion for ${S}_{0}^{\left(s\right)}\left(q,t\right)$ with P₀(t) = 1 in the e.o.m. for the tagged-particle correlator, equation (3), we obtain for the MSD

$$\begin{equation}\delta {r}_{0}^{2}\left(t\right)+{D}_{0}^{\left(s\right)}{\beta }^{2}{\int }_{0}^{t}{\zeta }_{0}^{\left(s\right)}\left(t-{t}^{\prime }\right)\delta {r}_{0}^{2}\left({t}^{\prime }\right)\mathrm{d}{t}^{\prime }=4{D}_{0}^{\left(s\right)}t,\end{equation} \tag{ 12 }$$

reminiscent of the corresponding e.o.m. in 2D systems [41]. However, in the case of quasi-confinement, a second equation for the long-wavelength limit of the effective self-memory kernel ${\beta }^{2}{\zeta }_{0}^{\left(s\right)}\left(t\right)={\mathrm{lim}}_{q\to 0}\enspace {q}^{2}{M}_{0}^{\left(s\right)}\left(q,t\right)$ is necessary. The irreducible memory kernel ${\boldsymbol{\mathcal{M}}}_{\mu }^{\left(s\right)}\left(q,t\right)$ for arbitrary mode index μ decouples in the long-wavelength limit, ${\mathrm{lim}}_{q\to 0}\enspace {\mathcal{M}}_{\mu }^{\alpha \beta ,\left(s\right)}\left(q,t\right)={\mathcal{M}}_{\mu }^{\alpha ,\left(s\right)}\left(t\right){\delta }_{\alpha \beta }$. The non-vanishing diagonal components are then given by

$$\begin{equation}{\mathcal{M}}_{\mu }^{{\Vert},\left(s\right)}\left(t\right)=\frac{n}{4\pi {L}^{3}}\int {k}^{3}\enspace \mathrm{d}k\enspace \sum\limits _{{\mu }_{1}}{\left({c}_{{\mu }_{1}}^{\left(s\right)}\left(k\right)\right)}^{2}{S}_{{\mu }_{1}}\left(k,t\right){S}_{\mu -{\mu }_{1}}^{\left(s\right)}\left(k,t\right),\end{equation} \tag{ 13 }$$

and

$$\begin{equation}{\mathcal{M}}_{\mu }^{\perp ,\left(s\right)}\left(t\right)=\frac{n}{2\pi {L}^{3}}\int k\enspace \mathrm{d}k\enspace \sum\limits _{{\mu }_{1}}{\left({Q}_{{\mu }_{1}}{c}_{{\mu }_{1}}^{\left(s\right)}\left(k\right)\right)}^{2}{S}_{{\mu }_{1}}\left(k,t\right){S}_{\mu -{\mu }_{1}}^{\left(s\right)}\left(k,t\right),\end{equation} \tag{ 14 }$$

enabling us to determine the equation for the effective memory kernel ${\beta }^{2}{\zeta }_{0}^{\left(s\right)}\left(t\right)$ in the long-wavelength limit

$$\begin{align}\hfill {D}_{0}^{\left(s\right)}{\beta }^{2}{\zeta }_{0}^{\left(s\right)}\left(t\right)& +{D}_{0}^{\left(s\right)2}{\beta }^{2}{\int }_{0}^{t}{\zeta }_{0}^{\left(s\right)}\left(t-{t}^{\prime }\right){\mathcal{M}}_{0}^{\perp ,\left(s\right)}\left({t}^{\prime }\right)\mathrm{d}{t}^{\prime }={D}_{0}^{\left(s\right)}{\mathcal{M}}_{0}^{{\Vert},\left(s\right)}\left(t\right)\hfill \\ \hfill & +{D}_{0}^{\left(s\right)2}{\int }_{0}^{t}{\mathcal{M}}_{0}^{{\Vert},\left(s\right)}\left(t-{t}^{\prime }\right){\mathcal{M}}_{0}^{\perp ,\left(s\right)}\left({t}^{\prime }\right)\mathrm{d}{t}^{\prime }.\hfill \end{align} \tag{ 15 }$$

The equation for the effective memory kernel is not restricted to μ = 0 relevant for the MSD, but can also be extended to arbitrary mode indices

$$\begin{align}\hfill {D}_{0}^{\left(s\right)}{Q}_{\mu }^{2}{M}_{\mu }^{\left(s\right)}\left(t\right)& +{D}_{0}^{\left(s\right)2}{\int }_{0}^{t}{Q}_{\mu }^{2}{M}_{\mu }^{\left(s\right)}\left(t-{t}^{\prime }\right){\mathcal{M}}_{\mu }^{{\Vert},\left(s\right)}\left({t}^{\prime }\right)\mathrm{d}{t}^{\prime }={D}_{0}^{\left(s\right)}{\mathcal{M}}_{\mu }^{\perp ,\left(s\right)}\left(t\right)\hfill \\ \hfill & +{D}_{0}^{\left(s\right)2}{\int }_{0}^{t}{\mathcal{M}}_{\mu }^{{\Vert},\left(s\right)}\left(t-{t}^{\prime }\right){\mathcal{M}}_{\mu }^{\perp ,\left(s\right)}\left({t}^{\prime }\right)\mathrm{d}{t}^{\prime },\quad \mu \ne 0,\hfill \end{align} \tag{ 16 }$$

with the long-wavelength limit ${M}_{\mu }^{\left(s\right)}\left(t\right)={\mathrm{lim}}_{q\to 0}\enspace {M}_{\mu }^{\left(s\right)}\left(q,t\right)$ for higher order effective self-memory kernels.

From the expansion for ${S}_{\mu }^{\left(s\right)}\left(q,t\right)$, equation (9), we get

$$\begin{equation}{\dot {P}}_{\mu }\left(t\right)+{D}_{0}^{\left(s\right)}{Q}_{\mu }^{2}{P}_{\mu }\left(t\right)+{D}_{0}^{\left(s\right)}{\int }_{0}^{t}{Q}_{\mu }^{2}{M}_{\mu }\left(t-{t}^{\prime }\right){\dot {P}}_{\mu }\left({t}^{\prime }\right)\mathrm{d}{t}^{\prime }=0,\end{equation} \tag{ 17 }$$

as the e.o.m. for P_μ (t). In contrast to the MSD, which only considers the dynamics in the parallel direction, it takes into account the dynamics perpendicular to the confining direction. The mathematical structure of this integro-differential equation is reminiscent of the e.o.m. for the tagged-particle correlator, equation (3), and therefore it is straightforward to solve it numerically.

The VACF Z₀(t) of a quasi-confined colloidal liquid can formally be connected to the MSD by the relation

$$\begin{equation}{Z}_{0}\left(t\right){:=}\frac{1}{4}\frac{{\mathrm{d}}^{2}}{\mathrm{d}{t}^{2}}\delta {r}_{0}^{2}\left(t\right).\end{equation} \tag{ 18 }$$

Therefore, the in-plane VACF can in principle be determined by taking derivatives of the MSD. Yet, due to discretization of the wavenumber grids in MCT, long-wavelength phenomena cannot be properly resolved. An exact e.o.m. can be readily derived using the Zwanzig projection operator formalism with the same irreducible effective memory kernel ${\beta }^{2}{\zeta }_{0}^{\left(s\right)}\left(t\right)$ as already introduced above in the e.o.m. for the MSD

$$\begin{equation}{Z}_{0}\left(t\right)+{D}_{0}^{\left(s\right)}{\beta }^{2}{\int }_{0}^{t}{\zeta }_{0}^{\left(s\right)}\left(t-{t}^{\prime }\right){Z}_{0}\left({t}^{\prime }\right)\mathrm{d}{t}^{\prime }=-{D}_{0}^{\left(s\right)2}{\beta }^{2}{\zeta }_{0}^{\left(s\right)}\left(t\right).\end{equation} \tag{ 19 }$$

We elaborate the long-time behavior of the VACF within the framework of MCT explicitly. For long times, the integrals for the irreducible memory kernel, equations (13) and (14), are dominated by small wavenumbers. In the long-wavelength limit, μ = 0 is the only mode which becomes slow and contributes to the sum. Thus, the perpendicular component vanishes, ${\mathcal{M}}_{0}^{\perp ,\left(s\right)}\left(t\right)\simeq 0$. Correspondingly, for long times the equation for the effective memory kernel, equation (15), reduces to ${D}_{0}^{\left(s\right)}{\beta }^{2}{\zeta }_{0}^{\left(s\right)}\left(t\right)={D}_{0}^{\left(s\right)}{\mathcal{M}}_{0}^{{\Vert},\left(s\right)}\left(t\right)$. Approximating the ISF and SISF for μ = 0 by collective diffusive motion, ${S}_{0}\left(q,t\right)\approx {S}_{0}\left(0\right)\mathrm{exp}\left(-{D}_{0}{q}^{2}t/{S}_{0}\left(0\right)\right)$ and self-diffusion ${S}_{0}^{\left(s\right)}\left(q,t\right)\approx \enspace \mathrm{exp}\left(-{D}_{0}^{L,\left(s\right)}{q}^{2}t\right)$, the long-time behavior of the VACF follows from the equation of motion, equation (19), using Tauber's theorem $Z\left(t\right)\simeq -{\left({D}_{0}^{L,\left(s\right)}\right)}^{2}{\beta }^{2}{\zeta }_{0}^{\left(s\right)}\left(t\right)$. The explicit representation of the long-time asymptote reads

$$\begin{equation}{Z}_{0}\left(t\right)\simeq -{\left({D}_{0}^{L,\left(s\right)}\right)}^{2}\frac{n}{4\pi {L}^{3}}\frac{{\left({c}_{0}^{\left(s\right)}\left(0\right)\right)}^{2}{S}_{0}\left(0\right)}{2{\left({D}_{0}^{L,\left(s\right)}+{D}_{0}/{S}_{0}\left(0\right)\right)}^{2}}{t}^{-2}.\end{equation} \tag{ 20 }$$

This power-law tail ∼−t⁻² is similar to the persistent anti-correlations found in bulk systems [38]. Due to the confinement, long-wavelength phenomena are only important along the lateral direction and we identify for long times the expected power-law for a two-dimensional system with a different prefactor.

In principle, autocorrelation functions similar to the VACF can also be defined for higher order modes, μ ⩾ 1, characterizing the out-of-plane dynamics with the same approach as for the MSD. However, in this work we focus on the long-time behavior of the VACF, since no long-time tails are expected for μ ≠ 0.

We perform event-driven molecular dynamics simulations for a 3D hard-sphere system, with periodic boundary conditions in all directions. The quasi-confinement is realized by setting the length of the z-direction of the simulation box much smaller than the other two directions. The extension of the unconfined directions is 25σ, whereas the confining length, L, has only a few particle diameters. Although the simulated particles are driven by Newtonian dynamics, we expect that the qualitative behavior, in particular for long times, will be comparable to MCT results for colloidal liquids undergoing Brownian motion [42–45]. To prevent crystallization polydisperse particles are used in the simulations. The particle diameters are given by an inverse-occupied volume distribution [8] with a standard deviation of 0.117σ, setting the polydispersity to 11.7%. It must be ensured that particles are initialized without overlap, which we do here by initializing particles in a larger volume, and compressing the system with standard algorithms [46, 47] to reach the desired packing fraction and system size [33]. The timescale of the simulation is determined by the thermal energy k_B T and the particle mass m, where we define ${t}_{0}=\sqrt{m{\sigma }^{2}/{k}_{\mathrm{B}}T}$. The SISF can directly be calculated from the density-mode representation, equation (1). In general, statistics for incoherent quantities can be significantly enhanced compared to collective ones by averaging over all particles.

As we have discussed previously [33], it is not practical to reach a fully equilibrated glassy fluid of polydisperse hard spheres in quasi-confinement, because as in the slit geometry [48], the particles will eventually demix and then crystallize. However, this process occurs on long timescales, outside what is experimentally relevant, and only accessible in simulations when relying on advanced, unphysical dynamical algorithms [8, 49]. Therefore, to study glassy dynamics here, we use Newtonian dynamics only, restrict ourselves to relatively short simulation times, and carefully check the results for any effects of aging. Our data is taken over a time 10⁵ t₀, after a sample preparation time 10⁵ t₀. We have also repeated all simulations using a longer sample preparation time of 10⁶ t₀, and found this does not change the qualitative behavior we will discuss in this section, and therefore the conclusions drawn are not affected by aging.

The simulation data for the SISF can be approximately described by the phenomenological Kohlrausch–William–Watts (KWW) stretched exponential [50]

$$\begin{equation}\text{KWW}\;\text{fit:}\;{S}_{\mu }^{\left(s\right)}\left(q,t\right)={F}_{\mu }^{\left(s\right)}\left(q\right)\mathrm{exp}\left\{-{\left[t/{\tau }_{\mu }^{\left(s\right)}\left(q\right)\right]}^{{\beta }_{\mu }^{\left(s\right)}\left(q\right)}\right\}.\end{equation} \tag{ 21 }$$

The data is fit in the time range t ∈ (10t₀, 10⁵ t₀) and the incoherent quantities, namely the Kohlrausch exponent ${\beta }_{\mu }^{\left(s\right)}\left(q\right)$, the relaxation time ${\tau }_{\mu }^{\left(s\right)}\left(q\right)$ and the self-nonergodicity parameter ${F}_{\mu }^{\left(s\right)}\left(q\right)$, can be extracted.

Unless otherwise stated, the numerical results are obtained within MCT working on an equidistant wavenumber grid $q=\hat{q}{\Delta}q+{q}_{0}$ with grid points $\hat{q}=0,\dots ,{N}_{q}-1$ parallel to the confinement. As input parameters q₀ σ = 0.1212, Δqσ = 0.404 and N_q = 100 are used. The discrete mode index associated with the confining direction is limited to |μ| ⩽ 15. We assume that the tagged particle is of the same species as the bath particles, therefore the identical direct correlation functions as in the collective case enter the MCT equations.

We start by analyzing the non-monotonic behavior of the nonergodicity parameters depending on the confinement length. We choose three different lengths L = 2.0σ, 2.5σ, 2.8σ where the effects are most pronounced. The non-monotonic behavior coincides qualitatively with the collective case [33], therefore we believe that discussing more different values does not contribute to a deeper understanding. Since it is known that MCT underestimates the critical packing fraction not only in the bulk, but also for confined systems [27, 51], MCT results for φ = 0.53 (figure 1 (dashed lines)) are compared to simulations with a higher packing fraction φ = 0.59 (figure 1 (solid lines)), therefore the dynamics are glassy in both cases. In principle, this discrepancy in the packing fraction can be reduced with modified MCT approaches [52], but it does not affect the validity of our conclusions in this work. Referring to the nonequilibrium-state diagram of quasi-confined liquids [33], MCT results are presented for φ = 0.53 in contrast to the remainder of this manuscript to ensure that the behavior is glassy for all investigated confinement lengths.

Zoom In Zoom Out Reset image size

Within MCT, the bell-shaped curves manifest non-monotonic behavior for all wavenumbers. The incommensurate value L = 2.5σ shows a higher structural arrest compared to the more commensurate values L = 2.0σ and L = 2.8σ. The confinement length is referred to as commensurate if the ratio L/σ is close to an integer and as incommensurate if it is close to a half-integer. The non-monotonic behavior reflects the collective behavior in quasi-confinement [33] and is also quite similar to the incoherent nonergodicity parameter within the slit geometry [51], although the non-monotonic behavior is less noticeable due to missing layering.

In contrast, for the simulations non-monotonic behavior is not obvious. Indeed, the structural arrest for L = 2.5σ is higher compared to L = 2.0σ, but it is almost indistinguishable from L = 2.8σ, where we expect a significant drop in the structural relaxation. Using the von Schweidler law [11] instead of KWW fits, a marginal difference can be observed. Yet, already in the slit case the non-monotonic behavior of the nonergodicity parameters is less pronounced for the simulations compared to the MCT results [51]. Since effects are considerably smaller in quasi-confinement anyway, the overlap of the curves for L = 2.5σ and L = 2.8σ is not surprising. This behavior indicates that the height of the plateau value of the SISF does not show a pronounced non-monotonic behavior, which is most likely due to the inclusion of polydispersity, in contrast to the MCT predictions and simulation data for the collective motion [33]. It has been shown for bulk systems [53] that the incoherent nonergodicity parameter is less sensitive to effects of different particle sizes of the constituents compared to the collective case. Therefore, we conclude that also effects of local order between commensurate and incommensurate packing have a minor influence on the incoherent nonergodicity parameter.

For both MCT and simulations the lowest order mode ${F}_{0}^{\left(s\right)}\left(q\right)$ reaches unity in the long-wavelength limit reflecting particle conservation in the system. For higher order modes, ${F}_{\mu }^{\left(s\right)}\left(q\to 0\right)$ corresponds to P_μ (t → ∞) according to equation (9), with a value different from unity in agreement with equation (10). P_μ (t) is only sensitive along the confining direction, and for large confinement length the nonergodicity parameter should approach the bulk limit. Then, it only depends on the magnitude of the 3D wavenumber $k=\sqrt{{q}^{2}+{Q}_{\mu }^{2}}$ and rotational invariance is restored, ${F}_{\mu }^{\left(s\right)}\left(q\right)\to {F}^{\left(s\right)}\left(k\right)$. The drop between ${F}_{0}^{\left(s\right)}\left(q\to 0\right)$ and ${F}_{1}^{\left(s\right)}\left(q\to 0\right)$ can be observed in figure 1 for simulations, whereas for MCT it becomes only visible for higher order modes, μ ⩾ 2, not included in the figure. Apart from the long-wavelength behavior the appearance of ${F}_{1}^{\left(s\right)}\left(q\right)$ is very similar to ${F}_{0}^{\left(s\right)}\left(q\right)$ indicating that the zero-order static quantities c₀(q) and S₀(q) are also the relevant quantities for the higher order self-nonergodicity parameters, since the static quantities c_μ (q) and S_μ (q) are significantly different between μ = 0 and μ = 1. Due to the periodic boundary conditions and the absence of layering, differences between the separate modes are less significant compared to the slit case [51].

We proceed with the full dynamic solution of the SISF. The unit of time for the MCT results is set by σ²/D₀. Figure 2(a) shows the first two modes ${S}_{0}^{\left(s\right)}\left(q,t\right)$ and ${S}_{1}^{\left(s\right)}\left(q,t\right)$ (inset) for MCT at packing fraction φ = 0.515, where it is known from the nonequilibrium-state diagram that a reentrant behavior occurs [33]. The related data for simulations at the higher packing fraction φ = 0.59 are displayed in figure 2(b).

Zoom In Zoom Out Reset image size

The MCT results look quite similar to the corresponding collective ISF [33] with clearly observable non-monotonic behavior depending on the confinement length. The MCT solution either decays to zero for the more commensurate values (L = 2.0σ, L = 2.8σ), identifying liquid states, or it reaches a non-vanishing plateau value corresponding to an arrested glassy state for the incommensurate value (L = 2.5σ). The increased longitudinal diffusion in the case of commensurate packing prevents slowing down of the dynamics, and therefore, the arrested glassy state is shifted to higher packing fractions compared to incommensurate packing [33].

The behavior of the first higher-order mode ${S}_{1}^{\left(s\right)}\left(q,t\right)$ almost duplicates the one of ${S}_{0}^{\left(s\right)}\left(q,t\right)$. Whereas ${S}_{0}^{\left(s\right)}\left(q,t\right)$ only considers the dynamics parallel to the confining direction, higher order modes of the SISF also include the perpendicular direction. The similarity between the modes of the SISF suggests that not only the long-time limits, but the whole dynamics are mainly driven by the zero-order mode of the static quantities, S₀(q) and c₀(q), in particular by the first sharp diffraction peak, like in the collective case [33]. Notably, ${S}_{1}^{\left(s\right)}\left(q,t\right)$ decays faster than ${S}_{0}^{\left(s\right)}\left(q,t\right)$ for all studied lengths.

In close resemblance to the collective ISF, the simulation data only show an intermediate plateau for all studied lengths, accompanied by a stretched relaxation. The general shape of the SISF coincides with the collective ISF with non-monotonic dependence on the confinement length. Yet, there are minor differences compared to the ISF, in particular the precise plateau values (incoherent nonergodicity parameters) differ as it has been discussed in the preceding section. Even though the non-monotonic dependence is not obvious in the plateau value, it still can be observed in the corresponding wavenumber dependent relaxation times ${\tau }_{0}^{\left(s\right)}\left(q\right)$ (figure 3(a)) and the Kohlrausch exponents ${\beta }_{0}^{\left(s\right)}\left(q\right)$ (figure 3(b)). It would be interesting to analyze if the large-wavenumber limit of the Kohlrausch exponent coincides with the von Schweidler exponent similar to other recent studies comparing the von Schweidler and the KWW equation [54]. Unfortunately, we cannot accurately determine the von Schweidler exponent form our data. Even for higher order modes, where the SISF almost overlap (inset of figure 2(b)), the related relaxation times ${\tau }_{1}^{\left(s\right)}\left(q\right)$ and Kohlrausch exponents ${\beta }_{1}^{\left(s\right)}\left(q\right)$ reveal considerable non-monotonic behavior (figure 3 (insets)). Therefore, we conclude that within the simulations the non-monotonic dependence is absent to a large extent in the plateau heights, but it is still present in the dynamics manifested by the relaxation times and the Kohlrausch exponents. Simulation data also confirms that higher modes of the SISF decay faster than the zero order quantity.

Zoom In Zoom Out Reset image size

The MSD is calculated for different confining lengths L at the same packing fractions φ = 0.515 for MCT (figure 4(a)) and φ = 0.59 for simulations (figure 4(b)) as above for the SISF.

Zoom In Zoom Out Reset image size

Quantitative differences in the limit of small t arise due to different microscopic dynamics. Whereas the MCT equations are derived for Brownian microscopic dynamics characterizing colloidal liquids, the simulations rely on Newtonian molecular dynamics. Therefore, for MCT, diffusive motion $\delta {r}_{0}^{2}\left(t\right)\simeq 4{D}_{0}t$, t → 0 can be observed for all investigated confinement lengths in the limit of small t, in contrast to the simulations, which show ballistic motion. With increasing time the dynamics are slowed down in the crossover regime as has already been discussed in connection with the SISF. For even larger times, particles are trapped due to the cage effect in the glassy state (L = 2.5σ), which becomes apparent in the MCT solution by a plateau value of the MSD. In contrast, the data for the other confinement lengths show liquid-like behavior with long-time diffusion $\delta {r}_{0}^{2}\left(t\right)=4{D}_{0}^{L,\left(s\right)}t$, t → ∞ after the crossover regime with a two-step relaxation process distinctive for MCT dynamics. This second rise is the α-process of the MSD, indicating the tagged-particle's exit of the cage. For simulations, particles move out of the cage eventually for all confinement lengths, and diffusive motion can be observed at long times. This behavior is consistent with the intermediate nature of the plateau in the SISF.

The inset of figure 4(a) shows the non-monotonic behavior of P₁(t) within MCT. The general shape of the relaxation is quite similar compared to the SISF for small q. The curves are indistinguishable for small times tD₀/σ² ≲ 10⁻¹ before the crossover regime sets in. All liquid curves decay to zero and the finite value for glassy states corresponds to ${F}_{1}^{\left(s\right)}\left(q\to 0\right)$ as already mentioned above.

A characteristic localization length ${r}_{0}^{\left(s\right)}$ can be defined by the small-q behavior of the nonergodicity parameter, ${F}_{0}^{\left(s\right)}\left(q\right)=1-{\left(q{r}_{0}^{\left(s\right)}\right)}^{2}+\mathcal{O}\left({q}^{4}\right)$ similar to the bulk case [55], which is related to the long-time limit of $\delta {r}_{0}^{2}\left(t\right)$ by

$$\begin{equation}\underset{t\to \infty }{\mathrm{lim}}\enspace \delta {r}_{0}^{2}\left(t\right)=4{r}_{0}^{\left(s\right)2}.\end{equation} \tag{ 22 }$$

In the glassy state (L = 2.5σ in figure 4) the finite value of the localization length ${r}_{0}^{\left(s\right)2}={\mathrm{lim}}_{t\to \infty }\enspace \delta {r}_{0}^{2}\left(t\right)/4$ related to the MSD corresponds to the size of cages in which the colloidal particles are trapped and reflects the behavior of colloidal suspensions in experiments for the slit geometry [22, 56], yet no physical boundaries are present in our quasi-confined system.

Table 1 indicates that the product of the long-time self-diffusion coefficient ${D}_{0}^{L,\left(s\right)}$ and the self-relaxation time ${\tau }_{0}^{\left(s\right)}$ at the structure factor peak, qσ = 7.0, is approximately constant, albeit the individual contributions vary significantly (by a factor of 50 for MCT and 3 for simulations, respectively) with the confinement length. Within MCT, we identify the relaxation time ${\tau }_{0}^{\left(s\right)}$ with the time where the SISF has reduced to 1/e. Since the dynamics are glassy within MCT for L = 2.5σ no self-relaxation time and self-diffusion coefficient can be extracted. MCT predicts a single α-relaxation scale, which implies the validity of the Stokes–Einstein (SE) relation for straight paths approaching the glass-transition line. However, since a reentrant phenomenon occurs between L = 2.0σ and L = 2.8σ, the SE relation does not have to hold here. Nevertheless, it is approximately fulfilled, albeit a slightly non-montonic behavior of ${\tau }_{0}^{\left(s\right)}{D}_{0}^{L,\left(s\right)}$ within the simulations. It is known that the violation of the SE relation sets in close to the critical packing fraction [57–60].

Table 1. Self-relaxation time ${\tau }_{0}^{\left(s\right)}$ for qσ = 7.0, long-time self-diffusion coefficient ${D}_{0}^{L,\left(s\right)}$ and the product ${\tau }_{0}^{\left(s\right)}{D}_{0}^{L,\left(s\right)}$ for different values of L at (a) packing fraction φ = 0.515 for MCT and (b) φ = 0.59 for simulations, respectively. For L = 2.5 the system is in the glassy state within MCT.

	L/σ	${\tau }_{0}^{\left(s\right)}/{\sigma }^{2}{D}_{0}^{-1}$	${D}_{0}^{L,\left(s\right)}/{D}_{0}$	${\tau }_{0}^{\left(s\right)}{D}_{0}^{L,\left(s\right)}/{\sigma }^{2}$
(a) MCT	2.0	1.11 × 102	1.41 × 10−4	1.57 × 10−2
	2.5	—	—	—
	2.8	5.33 × 103	2.62 × 10−6	1.40 × 10−2

	L/σ	${\tau }_{0}^{\left(s\right)}/{t}_{0}$	${D}_{0}^{L,\left(s\right)}/{\sigma }^{2}{t}_{0}^{-1}$	${\tau }_{0}^{\left(s\right)}{D}_{0}^{L,\left(s\right)}/{\sigma }^{2}$
(b) Simulations	2.0	1.36 × 104	1.81 × 10−5	2.46 × 10−1
	2.5	4.14 × 104	6.38 × 10−6	2.64 × 10−1
	2.8	1.76 × 104	1.43 × 10−5	2.45 × 10−1

Due to the dominating linear growth of diffusion some subtleties are usually hidden in the MSD, therefore, it is instructive to also investigate the VACF. Figure 5 shows the VACF Z₀(t) extracted from simulations for different confinement lengths. The prominent minimum in the VACF is related to the strong caging of particles for intermediate timescales. Quite similar curve shapes are present in the glass-like behavior of confined supercritical argon [61, 62] with inhomogeneous density profiles.

Zoom In Zoom Out Reset image size

As can be seen, non-monotonic behavior for short and intermediate times is only marginally pronounced, but still present (lower left inset). Since the crossover regime within the MSD for liquid states (figure 4(a)) extends over several orders of magnitude in time and just appears as an interpolation between the short-time diffusive motion $\delta {r}_{0}^{2}\left(t\right)=4{D}_{0}^{\left(s\right)}t,t\to 0$ and the long-time diffusion $\delta {r}_{0}^{2}\left(t\right)=4{D}_{0}^{L,\left(s\right)}t,t\to \infty$, we particularly want to investigate correlations of the VACF in the long-time regime. Unfortunately, due to insufficient statistics for the simulation data it is not possible to study long-time effects for the simulations, and we limit our discussion to MCT results. In principle, advanced simulation algorithms for colloidal bulk liquids have been developed recently to investigate the long-time behavior also for simulations [38].

The numerical evaluation within MCT is not straightforward and we have to rely on a logarithmic wavenumber grid similar to the 3D bulk case [38] to extract the long-time behavior. The reduced wavenumbers are then given by [x^−N , ..., x, 1]q_max σ with high-q cutoff q_max σ = 40, base x = 1.02 and N = 400 grid points. This corresponds to a minimal wavenumber q_min σ = 0.0148, ensuring that the long-time behavior of the VACF is properly resolved. Additionally, only |μ| ⩽ 10 modes are used to save computational resources.

In principle, there are more sophisticated grid types which have been discussed lately for two- and three-dimensional systems [63]. The main idea behind these nonuniform grids is to have a better representation with more grid points near the maxima and minima of the structure factor and the direct correlation function, but less wavenumbers in regions where the behavior is more linear. The implementation of this scheme for quasi-confined liquids would require a different wavenumber grid for each mode μ, since the extreme values of the static input sensitively depend on the wavenumber for different mode indices [32]. Still, we find for the long-time behavior of the VACF that it is sufficient to have a good long-wavelength resolution since the corresponding integrals are dominated by low wavenumbers instead of the structure factor peaks, and therefore, we simply rely on a logarithmic grid independent of the mode index.

At the arbitrarily chosen confinement length L = 2.0σ the numerical results for the reduced VACF −Z₀(t)/φ for MCT display non-trivial correlations beyond the crossover regime of the MSD with anti-correlated long-time tails Z₀(t) ∼ −t⁻² for all investigated packing fractions, see figure 6. Yet, in contrast to the bulk behavior, Z(t) ∼ −t^−5/2, the exponent indicates two-dimensional dynamics with the prefactor sensitively depending on the density, equation (20). Surprisingly, the non-monotonic dependence on the confinement length only becomes apparent for larger packing fractions φ ≳ 0.45 approaching the glass-transition (inset). For this comparatively high packing fraction non-monotonic behavior depending on the confinement length can be observed in two respects. Firstly, the crossover regime is shifted to larger times for L = 2.5σ compared to the other confining lengths and secondly, the tails are more suppressed for incommensurate values. This demonstrates that evidence of glassy-like behavior can be observed even at φ = 0.47, since these are the same values of L which correspond to glassy states at higher packing fraction. Indeed, the strong L³ dependence in the denominator of equation (20) is compensated for by the L-dependent structural quantities S₀(0), c₀(q) and the long-time self-diffusion coefficient ${D}_{0}^{L,\left(s\right)}$ if the density is not too high.

Zoom In Zoom Out Reset image size

For large confinement lengths (L ≳ 5.0σ) the dynamics of the particle are indistinguishable from bulk motion for short times with an intermediate window, where the VACF decays with a power-law, Z(t) ∼ −t^−5/2. Only for longer times, a crossover from the unconfined dynamics occurs and the expected power law for quasi-confinement, Z(t) ∼ −t⁻², is observed, similar to the tagged-particle within a slit [64]. Extreme confinement, however, cannot be considered as the two dimensional limit case since due to periodic boundary conditions the hard spheres are still able to move perpendicular to the confining direction.

Although, the long-time behavior in the simulations is covered by noise, indications of a similar intermediate tail Z(t) ∼ −t⁻² can be observed in the upper right inset of figure 4. Indeed, for much longer times the classical hydrodynamic Alder tails [37] are expected for a momentum-conserving system of Newtonian hard-spheres. A comparable intriguing crossover from positive correlations to anti-correlations has been reported for dense hard-sphere systems [65].

All MCT data show a divergent short-time behavior, Z(t) ∼ −t^−1/2 as t → 0, related to the hard-sphere potential. This divergence is not restricted to quasi-confined systems and can be observed in colloidal liquids with arbitrary dimension, e.g. in 3D bulk liquids [38]. The exponent does not depend on the confining geometry, contrary to the long-time tail.

MCT also gives a qualitative explanation for the strong variation of the long-time tail in terms of the static compressibility and slowing down of transport related to the 3D bulk case [38]. In accordance with equation (20) and [38] the tail is suppressed by orders of magnitude with increasing φ. In the dilute limit the self-diffusion coefficient in the denominator can be ignored and the prefactor of the long-time tail shows a direct dependence on the long-wavelength limit of the static structure factor, S₀(0). Approaching the glass transition the tails get even more suppressed due to the singularity of ${D}_{0}^{L,\left(s\right)}$. At intermediate times we then expect to observe two additional power laws, which are directly related to the critical beta decay [64].