Static properties of quasi-confined hard-sphere fluids

Charlotte F. Petersen [email protected] Lukas Schrack Thomas Franosch Institut für Theoretische Physik, Leopold-Franzens-Universität Innsbruck, Technikerstraße 21A, A-6020 Innsbruck, Austria.

Abstract

Confined fluids display complex behavior due to layering and local packing. Here, we disentangle these effects by confining a hard-sphere fluid to the surface of a cylinder, such the circumference extends only over a few particle diameters. We compare the static structure factor and the pressure measured in computer simulations to the Percus-Yevick closure in liquid state theory. A non-monotonic evolution of the static-structure-factor peak and the pressure is observed upon variation of the confining length, similar to a liquid confined between two plates. This indicates that the density profile and the particle correlations may not be intrinsically connected in real confined liquids.

Keywords: Structural correlations, Molecular dynamics

I Introduction

Confined liquids are prevalent in nature Beebe et al. (2002) and have many industrial applications Urbakh et al. (2004); Pinilla et al. (2005). Physical confinement is known to affect fundamental properties of a fluid, such as diffusion Mittal et al. (2008), the glass transition Varnik and Franosch (2016); Lang et al. (2010); Krakoviack (2005a, b), and the freezing transition Schmidt and Löwen (1996, 1997). Confinement induces ordering of the particles, which results in density modulation Hansen and McDonald (2013); Bunk et al. (2007). However, the local properties of the confined fluid are not simply related to bulk fluids with equivalent densities Goel et al. (2008); Mittal et al. (2008); Goel et al. (2009); Bollinger et al. (2015). A slit geometry, where a fluid is confined between two parallel walls, has been extensively used in investigations into the effects of confinement experimentally Nugent et al. (2007); Nygård et al. (2012, 2013, 2016, 2017); Nygård (2017), theoretically Deb et al. (2011); Mandal et al. (2014, 2017), and also computationally Mittal et al. (2008); Mandal and Franosch (2017); Geigenfeind et al. (2015). In this layered fluid, the density oscillates as a function of distance from the walls. For extreme confinement, diffusion is slowed down dramatically compared to the bulk Scheidler et al. (2004); Gokhale et al. (2016); Eral et al. (2009) and remarkably, both the diffusivity and the glass transition show a non-monotonic dependence on the plate-to-plate distance Mandal et al. (2014); Lang et al. (2010, 2012, 2013, 2014a). These non-monotonic effects are a result of a subtle interplay between the structure of the local cage created by neighbors, and the layering induced by the walls. There exists empirical evidence that these unusual long-range transport properties are intrinsically related to structure, through a scaling between diffusivities and excess entropies Ingebrigtsen et al. (2013). Similarly, the mode-coupling theory of the glass transition Götze (2009) predicts dynamic properties using only structural information as input. Thus, the relevance of static correlations to transport properties provides a strong motivation to further investigate the structure of fluids under confinement.

Apart from the inhomogeneous density profile, the structure of a confined liquid is also characterized by its pair-distribution function or equivalently its structure factor, which also show non-monotonic behavior with confinement length Mandal et al. (2014). The structure factor is a key measure of the correlations in a fluid, as it can be calculated in computer simulation, theoretically, and measured directly in scattering experiments. In liquid state theory a common approach is to use integral equations to formulate closure relations, the simplest being the Percus-Yevick (PY) approximation Hansen and McDonald (2013); Percus and Yevick (1958); Thiele (1963); Wertheim (1963). In bulk hard-sphere fluids an analytic solution of the structure factor exists in odd dimensions Freasier and Isbister (1981); Leutheusser (1984), and aside from 3D space, it has also been studied extensively in 1D Wertheim (1964), 5D Freasier and Isbister (1981); Leutheusser (1984); González et al. (1990), and 7D Robles et al. (2004). In 2D the PY equations have been solved numerically Lado (1968); Leutheusser (1986), and more recently using a convenient method which allows the results for all densities to be calculated at once Adda-Bedia et al. (2008). Unfortunately, the structure factors calculated with the PY approximation deviate somewhat from those calculated accurately with computer simulations Hansen and McDonald (2013). Yet, it still proves useful, since parameters can be readily changed without running new simulations and the solutions can be evaluated to arbitrary precision. The PY theory has been extended to confined systems, which are more complicated because they are inhomogeneous. The statics have been solved for hard spheres confined to a slit Sokołowski (1980); Kjellander and Sarman (1991); Boţan et al. (2009); Nygård et al. (2013), compared with experiment Nygård et al. (2012), and used to predict the dynamics Mandal et al. (2014); Lang et al. (2010).

In this confined system, the density profile changes as the degree of confinement is varied, and the correlations in the liquid display a non-trivial dependence on this confinement length. Thus a natural question arises; are the unusual correlations in the confined liquid are a mere reflection of the layering, or are they an independent manifestation of the confinement itself?

Refer to caption — Figure 1: Illustration of the 2D analog of our simulated system. (a) Periodic simulation cells. Space is tiled periodically with one simulation cell, such that particles at the top of the cell interact with particles at the bottom of the cell, and likewise for the left and right edges. The length of one side of the periodic box (the $z$ -direction) is much shorter than the other, almost bulk direction. The 3D simulations have 1 short dimension and 2 long dimensions, effectively a periodic slab. (b) The periodic cells in (a) can be represented as a single simulation cell, wrapped around the surface of a torus.

In this paper we investigate confinement in a liquid with uniform density. We do this by considering a 3D periodic liquid, that is, one which is confined to the surface of a torus or long cylinder. This quasi-confined fluid is illustrated in figure 1 for the two-dimensional analog. This model has the advantage of being translationally invariant, but confinement on the order of the particle diameter can still be imposed, allowing us to investigate strong confinement effects without the complications of walls and density modulation.

Similar systems have previously been used as models of confinement, for example in the study of non-equilibrium dynamics in a narrow channel Bénichou et al. (2013, 2016). In the turbulence community, a similar method of reducing the complexity of a system is used, referred to as Fourier-mode reduction or wavenumber decimation, to investigate the effects of local versus non-local dynamics Grossmann et al. (1996); Meneguzzi et al. (1996); De Lillo and Eckhardt (2007); Laval et al. (2001). A similar periodic geometry is used in finite size scaling analysis to probe the length scales relevant to the glass transition Karmakar et al. (2009); Berthier et al. (2012), but with a cubic periodic cell, where the edge length is much larger than a particle diameter. An anisotropic unit cell geometry has been used to calculate the critical Casimir force in a binary mixture close to criticality Puosi et al. (2016). In molecular dynamics simulations, finite size effects from the limited length of periodic simulation cells are known to cause anisotropy in calculated pair-distribution functions Denton and Egelstaff (1997); Pratt and Haan (1981a, b) and stress González-Melchor et al. (2005). These effects are particularly significant in non-equilibrium simulations where the periodic boundaries move Bernardi et al. (2015); Petravic and Evans (1998). Until now, these effects have usually been considered as an undesirably consequence of finite computing resources, rather than as a deliberate method of confinement. A more thorough understanding of the properties of quasi-confined liquids, as the degree of confinement approaches the particle diameter, will be important to gaining a more thorough understanding of the properties of the more complex case of confined liquids.

We study the static properties of quasi-confined hard-sphere fluids using theory and event-driven simulations. The quantities relevant to the study of quasi-confined fluids are presented in detail in section II, along with the integral equation theory and simulation methods. The results are presented in section III, were we find a good quantitative match between the PY theory and the simulations at low packing fractions, and semi-quantitative agreement for higher densities. For both methods we observe non-monotonic behavior of the static structure and the pressure with increasing confinement length, similar to the behavior of confined liquids in a slit geometry. In section IV we discuss the conclusions of the work.

II Methods

II.1 Theory of quasi-confined fluids

We consider $N$ hard spheres of diameter $\sigma$ in three dimensions using periodic boundary conditions, where one dimension is much smaller than the other two, shown schematically in figure 1(a). This can be thought of as spheres confined to the surface of a four dimensional torus, illustrated in 2D in figure 1(b). The small periodic dimension is designated as the $z$ -direction, such that $-L/2\leq z\leq L/2$ . We anticipate the thermodynamic limit, $N\to\infty$ , with lateral box size $L_{\text{box}}\to\infty$ such that the area density $n_{0}=N/L_{\text{box}}^{2}$ remains fixed. Thus, these fluids are a special case of the more general class of layered fluids, such as fluids next to a wall Roth (2010), in a slit geometry or fluids subject to a one-dimensional potential Bollinger et al. (2014); Saw and Dasgupta (2016); Herrera-Velarde and Castañeda-Priego (2007); Evers et al. (2013); Capellmann et al. (2018), which display translational symmetry along the $x$ – $y$ direction as well as rotational symmetry around the $z$ -axis.

The structural properties of layered fluids are encoded in averages of the multi-particle densities. Here we focus on the lowest order quantities, the simplest being the fluctuating density defined by

\rho(\mathbf{r},z)=\sum_{i=1}^{N}\delta(\mathbf{r}-\mathbf{r}_{i})\delta(z-z_{i}),

(1)

where $(\mathbf{r}_{i},z_{i})$ are the lateral and transverse coordinates of particle $i$ . For these layered fluids the equilibrium density $n(z):=\langle\rho(\mathbf{r},z)\rangle$ depends only on the transverse direction, where $\langle\cdot\rangle$ denotes ensemble averaging. For our quasi-confined liquid, this equilibrium density is even uniform $n(z)=n\equiv n_{0}/L$ .

The next simplest property is the two-particle density,

	$\displaystyle\rho^{(2)}(\mathbf{r},\mathbf{r\,}^{\prime},z,z^{\prime})$	$\displaystyle=\mathop{\sum\sum}_{i\neq j}^{N}\delta(\mathbf{r}-\mathbf{r}_{i})\delta(z-z_{i})\delta(\mathbf{r\,}^{\prime}-\mathbf{r}_{j})\delta(z^{\prime}-z_{j})$
		$\displaystyle=\rho(\mathbf{r},z)\rho(\mathbf{r\,}^{\prime},z^{\prime})-\delta(\mathbf{r}-\mathbf{r\,}^{\prime})\delta(z-z^{\prime})\rho(\mathbf{r},z).$		(2)

Its average for layered fluids depends only on the relative lateral distance $|\mathbf{r}-\mathbf{r\,}^{\prime}|$ , and this property in inherited by all two-point correlation functions. Then, the general definition of the pair-distribution function Hansen and McDonald (2013) reduces for layered fluids to

\displaystyle g(|\mathbf{r}-\mathbf{r\,}^{\prime}|,z,z^{\prime}):=\frac{\langle\rho^{(2)}(\mathbf{r},\mathbf{r\,}^{\prime},z,z^{\prime})\rangle}{n(z)n(z^{\prime})},

(3)

and provides readily interpretable information on the local packing.

In contrast, the density-density correlation function

\displaystyle G(|\mathbf{r}-\mathbf{r\,}^{\prime}|,z,z^{\prime}):=\frac{1}{n_{0}}\langle\delta\rho(\mathbf{r},z)\delta\rho(\mathbf{r\,}^{\prime},z^{\prime})\rangle,

(4)

where $\delta\rho(\mathbf{r},z)=\rho(\mathbf{r},z)-n(z)$ , is the fluctuating density, can be measured in scattering experiments and reflects layering as well as packing. Although both quantities can be related, for general layered fluids the density profile $n(z)$ explicitly enters. However, for quasi-confinement both quantities are trivially related, which can be made explicit in terms of the total correlation function $h:=g-1$ (which decays to zero for large distances). Then equations (3) and (4) yield

h(|\mathbf{r}-\mathbf{r\,}^{\prime}|,z,z^{\prime})=\frac{L^{2}}{n_{0}}G(|\mathbf{r}-\mathbf{r\,}^{\prime}|,z,z^{\prime})-\frac{1}{n}\delta(\mathbf{r}-\mathbf{r}^{\prime})\delta(z-z^{\prime}),

(5)

and the dependence on the transverse coordinates reduces to the relative distance $z-z^{\prime}$ . Furthermore, reflection symmetry reveals that it depends only on $|z-z^{\prime}|$ . Hence, the geometry of quasi-confinement, where confinement is introduced without walls or layering, is ideally suited to disentangle local packing and layering and greatly simplifies the analysis.

Here, we follow the strategy developed for the slit geometry Lang et al. (2010, 2012, 2013, 2014b) and expand the dependence on the transverse positions in all quantities in terms of discrete Fourier modes $\exp{\left(-iQ_{\mu}z\right)}$ in the $z$ -direction ( $Q_{\mu}=2\pi\mu/L,\mu\in\mathbb{Z}$ ). In contrast, the spatial dependence parallel to the $x$ – $y$ plane is decomposed into ordinary plane waves $\exp(-i\mathbf{q}\cdot\mathbf{r})$ with $\mathbf{q}=(q_{x},q_{y})$ . These wave vectors are taken as discrete initially $(q_{x},q_{y})\in(2\pi/L_{\text{box}})\mathbb{Z}^{2}$ , however in the thermodynamic limit they become continuous variables, such that sums are replaced by integrals $(1/A)\sum_{\mathbf{q}}\ldots\mapsto(2\pi)^{-2}\int\mathrm{d}^{2}q\ldots$ as usual, with $A=L_{\text{box}}^{2}$ . Then the following orthogonality and completeness relations hold

	$\displaystyle\frac{1}{A}\int_{A}e^{i(\mathbf{q}-\mathbf{q\,}^{\prime})\cdot\mathbf{r}}\mathrm{d}^{2}r=\delta_{\mathbf{q},\mathbf{q\,}^{\prime}},$		(6)
	$\displaystyle\frac{1}{A}\sum_{\mathbf{q}}e^{i\mathbf{q}\cdot(\mathbf{r}-\mathbf{r\,}^{\prime})}=\delta(\mathbf{r}-\mathbf{r\,}^{\prime}),$		(7)
	$\displaystyle\frac{1}{L}\int_{-L/2}^{L/2}\exp{\left[i(Q_{\mu}-Q_{\mu^{\prime}})z\right]}\mathrm{d}z=\delta_{\mu,\mu^{\prime}},$		(8)
	$\displaystyle\frac{1}{L}\sum_{\mu}\exp{\left[iQ_{\mu}(z-z^{\prime})\right]}=\delta(z-z^{\prime}).$		(9)

The microscopic density is now

\rho(\mathbf{r},z)=\frac{1}{A}\sum_{\mathbf{q}}\frac{1}{L}\sum_{\mu}\rho_{\mu}(\mathbf{q})\exp{\left(-iQ_{\mu}z\right)}e^{-i\mathbf{q}\cdot\mathbf{r}},

(10)

where the expansion coefficients

\rho_{\mu}(\mathbf{q})=\sum_{i=1}^{N}\exp{\left(iQ_{\mu}z_{i}\right)}e^{i\mathbf{q}\cdot\mathbf{r}_{i}}

(11)

take the role of the fundamental observables. By translational invariance in the unbounded directions

\langle\rho_{\mu}(\mathbf{q})\rangle=An_{\mu}\delta_{\mathbf{q},0}\,,

(12)

where $n_{\mu}=\int_{-L/2}^{L/2}n(z)\exp(iQ_{\mu}z)\mathrm{d}z$ is the Fourier coefficient of the density profile. For quasi-confinement only $n_{0}$ is non-vanishing, corresponding to the homogeneous density.

Translational symmetry in the $x$ – $y$ direction shows that the modes $\rho_{\mu}(\mathbf{q})$ are orthogonal for different wavenumbers. Thus the correlations between such modes are characterized by the generalized structure factor

S_{\mu\nu}(q)=\frac{1}{N}\langle\rho_{\mu}(\mathbf{q})^{*}\rho_{\nu}(\mathbf{q})\rangle.

(13)

By rotational invariance around the $z$ -axis the structure factors also depend only on the modulus $q=|\mathbf{q}|$ of the wave vector. For quasi-confinement, translational invariance along the z-direction implies that only the diagonal elements survive

S_{\mu}(q)\equiv S_{\mu\mu}(q)=\frac{1}{N}\langle|\rho_{\mu}(\mathbf{q})|^{2}\rangle.

(14)

Furthermore, in this case $S_{-\mu}(q)=S_{\mu}(q)$ as a consequence of reflection symmetry at the $x$ – $y$ plane.

The mode expansion of the density-density correlation function for layered fluids reads

\displaystyle G(|\mathbf{r}-\mathbf{r\,}^{\prime}|,z,z^{\prime})

\displaystyle=\frac{1}{A}\sum_{\mathbf{q}}\frac{1}{L^{2}}\sum_{\mu\nu}S_{\mu\nu}(q)e^{i\mathbf{q}\cdot(\mathbf{r}-\mathbf{r\,}^{\prime})}\exp{\left[i(Q_{\mu}z-Q_{\nu}z^{\prime})\right]},

(15)

where the symmetries carry over from real to Fourier space. Reversely, the generalized structure factor can be obtained evaluating the Fourier coefficients of $G(|\mathbf{r}-\mathbf{r\,}^{\prime}|,z,z^{\prime})$ . For quasi-confinement, the double sum over the mode indices collapses, drastically simplifying the situation.

II.2 Integral equation theory

The static structure factor can be obtained approximately by first introducing the direct correlation function via a suitably adapted Ornstein-Zernike (OZ) relation. Then a closure relation connects the direct correlation function back to the pair-distribution function. Here, we rely on the Percus-Yevick closure, which is particularly simple and also rather accurate for hard spheres Hansen and McDonald (2013); Henderson (1992).

The Ornstein-Zernike equation adapted for a quasi-confined liquid reads

\displaystyle h(|\mathbf{r}-\mathbf{r\,}^{\prime}|,z,z^{\prime})=

\displaystyle c(|\mathbf{r}-\mathbf{r\,}^{\prime}|,z,z^{\prime})+n\int\mathrm{d}^{2}r^{\prime\prime}\int_{-L/2}^{L/2}\mathrm{d}z^{\prime\prime}c(|\mathbf{r}-\mathbf{r\,}^{\prime\prime}|,z,z^{\prime\prime})\times h(|\mathbf{r\,}^{\prime\prime}-\mathbf{r\,}^{\prime}|,z^{\prime\prime},z^{\prime}),

(16)

and relates the total correlation function $h(|\mathbf{r}-\mathbf{r\,}^{\prime}|,z,z^{\prime})=g(|\mathbf{r}-\mathbf{r\,}^{\prime}|,z,z^{\prime})-1$ to the direct correlation function $c(|\mathbf{r}-\mathbf{r\,}^{\prime}|,z,z^{\prime})$ . Equation (16) is a special case of the more general Ornstein-Zernike relation for layered fluids.

By the convolution theorem the OZ relation simplifies in Fourier space

\displaystyle\hat{h}(q,z,z^{\prime})-\hat{c}(q,z,z^{\prime})

\displaystyle=n\int_{-L/2}^{L/2}\mathrm{d}z^{\prime\prime}~{}\hat{c}(q,z,z^{\prime\prime})\hat{h}(q,z^{\prime\prime},z^{\prime}),

(17)

where we introduced the Hankel transform

	$\displaystyle\hat{h}(q,z,z^{\prime})$	$\displaystyle=\int\mathrm{d}^{2}re^{i\mathbf{q}\cdot\mathbf{r}}h(\|\mathbf{r}\|,z,z^{\prime})$
		$\displaystyle=\int_{0}^{\infty}\mathrm{d}r{\mathrm{J}}_{0}(qr)rh(r,z,z^{\prime}),$		(18)

where ${\mathrm{J}}_{0}(\cdot)$ is the Bessel function of the first kind of order 0. Similarly, $\hat{c}(q,z,z^{\prime})$ is given in the same way.

Here, we supplement the OZ equation by the standard Percus-Yevick closure relation in real space

\displaystyle c=f\left(1+h-c\right),

(19)

with $f=\exp(-\beta U)-1$ denoting the Mayer function Hansen and McDonald (2013). For hard-core interactions the Mayer function evaluates to $-1$ for overlap and $0$ otherwise. Hence, equations (16) to (19) form a closed set of equations, which can be solved numerically.

The connection of the static structure factor, equation (14) to the total correlation function, equation (5) reads upon decomposition into Fourier modes

S_{\mu}(q)=1+\frac{n_{0}}{L^{2}}h_{\mu}(q),

(20)

where only the diagonal component

\displaystyle\hat{h}_{\mu}(q)\equiv\hat{h}_{\mu\mu}(q)=\int_{-L/2}^{L/2}\mathrm{d}z

\displaystyle\int_{-L/2}^{L/2}\mathrm{d}z^{\prime}\hat{h}(q,z,z^{\prime})\times\exp[-iQ_{\mu}(z-z^{\prime})]

(21)

enters.

We solve equations (17) to (19) iteratively. The Hankel transform and its corresponding inverse can be calculated by introducing logarithmic grids. Then equation (II.2) becomes a convolution integral, which can be computed efficiently using a fast Fourier transform Talman (1978); Hamilton (2000). We calculate the integrals for distances $r/\sigma$ and wave numbers $q\sigma$ lying within the range $[10^{-4},10^{4}]$ using 2048 grid points. To avoid an amplification of numerical errors during the iteration, cutoffs $r_{\text{min}}/\sigma=q_{\text{min}}\sigma=2.5\times 10^{-3}$ and $r_{\text{max}}/\sigma=q_{\text{max}}\sigma=2.5\times 10^{3}$ are introduced, where $h$ and $c$ are put to zero outside of the intervals $[r_{\text{min}}/\sigma,r_{\text{max}}/\sigma]$ and $[q_{\text{min}}\sigma,q_{\text{max}}\sigma]$ . The remaining $z$ integral in equation (17) is calculated using a simple trapezoidal rule on a uniformly discretized grid with spacing $0.01\sigma$ . The quasi-confinement is realized by using periodic boundary conditions.

II.3 Simulation

We perform event-driven simulations of hard spheres in three dimensions Rapaport (2004), in the periodic geometry illustrated in figure 1. We have simulated systems with confinement lengths $L=2.1\sigma$ up to $L=15\sigma$ . The size of the periodic box in the other two lateral dimensions is approximately $L_{\text{box}}=80\sigma$ , depending on the packing fraction and confinement length. This length was chosen to give suitable $q$ resolution in the structure factor calculations. We consider densities of $n=$ 0.35, 0.7 and 0.85, corresponding to packing fractions $\varphi:=N\pi\sigma^{3}/6AL$ of $\varphi=0.183,~{}0.37$ and 0.445. Dependent on $\varphi$ , $L_{\text{box}}$ and $L$ , the number of particles in each simulation ranges from 7,200-64,000. Due to the hard-sphere interaction, the thermal energy $k_{B}T$ does not affect the static properties, and enters only via the timescale of the simulation $t_{0}=\sqrt{m\sigma/k_{B}T}$ , where $m$ is the particle mass. The simulations are initialized with the particles in a periodic arrangement and velocities assigned randomly, but the system is allowed to equilibrate before any measurements are made. We checked that the results do not change if a longer equilibration period is used. The dynamics at the packing fractions under consideration are not glassy, and the correlations decay on the order of $10t_{0}$ . To ensure adequate statistics, 100 independent runs of each simulation are performed. In each of these repetitions, the simulations are run for a time $10^{3}t_{0}$ – $6\times 10^{3}t_{0}$ after the equilibration period, and the configuration is sampled every $50t_{0}$ . These configurations are used to calculate time-averaged static properties of the fluid.

III Results and Discussion

We first consider the real-space correlations, to ensure that our simulations have not crystallized at the packing fractions under consideration. The structure of the fluid can be measured in real space with the pair-distribution function, $g(|\mathbf{r}-\mathbf{r\,}^{\prime}|,z,z^{\prime})$ given by equation (3). Due to the symmetry of our system, we choose to sum only over particles which have approximately the same $z$ coordinate, $z=z^{\prime}$ . This correlation is also independent of the direction of $\mathbf{r}-\mathbf{r^{\prime}}$ , so we can average over all angles to give us a 2D radial distribution function, measuring the pair correlations perpendicular to the confinement direction. The pair-distribution function shows that the correlations between particles in the simulation die out at long distances, seen in figure 2, implying that the fluid has not crystallized. Already in this correlation measure we see non-monotonic effects of confinement, demonstrated by the contact value, plotted in the insert of figure 2. Of the three values of confinement length considered in the main plot, the contact value is lowest for the middle degree of confinement.

The main measure we use to characterize the structure of the fluid is the static structure factor, equation (14). It is a measure of the correlations in the system in Fourier space, and is a natural choice in the study of fluids because it is measurable directly experimentally by scattering experiments. It is also the starting point for dynamic theories of liquids and glasses, such as the mode-coupling theory Götze (2009). The lowest-mode structure factor $S_{0}(q)$ depends only on the lateral coordinates and correspondingly we refer to it as the in-plane structure factor. In contrast, the higher modes $S_{\mu}(q),\mu>0$ are sensitive to the arrangement along the confinement direction and reveal valuable information on the competition between local packing and confinement.

At low packing fraction, we find that the in-plane structure factor $S_{0}(q)$ is virtually independent of the confinement length $L$ , seen in figure 3(a). The effect of confinement can only be seen in higher modes, demonstrated in figure 3(d) and (g). In each of these cases the Percus-Yevick closure is a very close match to the simulations. At higher packing fractions we begin to observe quantitative differences between the structure factor calculated with the PY closure and simulations, seen in the right two columns of figure 3. Yet, the curve shapes are similar and the wave vector corresponding to each peak matches. Nevertheless, the height of the first peak is quite different, particularly at the highest packing fraction considered. However, it is at these higher packing fractions ( $\varphi=0.37$ and 0.445) where we also start to see interesting effects due to the confinement. In particular, the peak of the in-plane structure factor shows a non-monotonic dependence on the confinement length, seen in figure 4. The curves for the simulation and PY closure are qualitatively similar, with oscillations occurring at almost the same values of $L$ , but the oscillations are significantly less pronounced in the theoretical results. At large values of $L$ the structure factor peak appears to saturate at a constant value for both simulation and theory, as expected for approaching the bulk. Consistent with the literature, the PY approximation overestimates the height of the first peak of the static structure factor in the bulk Hansen and McDonald (2013). The height of the main peak of the structure factor is a measure of the near-ordering of the fluid. The oscillatory behavior of this peak suggests that the degree of ordering in the fluid changes from commensurate to incommensurate packing. This interpretation is corroborated by the fact that the period of the oscillates coincides with the hard-sphere diameter. In the simulation, the maximum appears at $L=2.3\sigma$ . This is the same behavior as liquids confined between hard walls Mandal et al. (2014).

The constant value of the peak of the in-plane static structure factor at large $L$ implies that the system is approaching the bulk limit. Then we expect that 3D rotational invariance is restored such that only the magnitude of the 3D wave vector $\mathbf{k}=(\mathbf{q},Q_{\mu})$ determines the structure factor

\displaystyle S_{\mu}(q)\to S(k)\qquad\text{with }k=\sqrt{q^{2}+Q_{\mu}^{2}},

(22)

as $L\to\infty$ . This convergence is confirmed in figure 5 for $L=5\sigma$ for the lowest three modes of the structure factors. We can also observe from this plot that while the differences between simulation and the PY closure are smaller than in the quasi-confined fluids with higher degrees of confinement, they still grow with increasing packing fraction, consistent with bulk behavior Hansen and McDonald (2013).

Given the non-monotonic behaviour of the pair-distribution function and static structure factor, we would also expect to see oscillations in the pressure of the fluid as a function of confinement length. For a hard-sphere fluid without rotational symmetry, the pressure is measured by the diagonal components of the pressure tensor in the transverse and lateral directions De Miguel and Jackson (2006). Its form can be derived starting from thermodynamics, and was presented already by Irving and Kirkwood Irving and Kirkwood (1950). In terms of inter-particle forces, this expression is given by

\displaystyle P_{zz}

\displaystyle=nk_{B}T-\frac{n^{2}}{2}\int_{A}\int_{-L/2}^{L/2}\frac{z^{2}}{\sqrt{r^{2}+z^{2}}}u^{\prime}(\sqrt{r^{2}+z^{2}})g(r,z)\textrm{d}z\textrm{d}^{2}r

(23)

for the transverse component, in our quasi-confined geometry. Here, $u$ is the pair potential and $g(r,z)$ is the pair-distribution function, which depends only on the in-plane and transverse distances $r$ and $z$ . Following the approach outlined in Hansen and McDonald Hansen and McDonald (2013), this can be written in a form suitable for calculation in a hard-sphere fluid as

\displaystyle P_{zz}

\displaystyle=nk_{B}T+k_{B}T\pi n^{2}\int_{-\sigma}^{\sigma}g(\sqrt{\sigma^{2}-z^{2}},z)z^{2}\textrm{d}z.

(24)

Here the pair-distribution function has to evaluated at contact, while the $z$ -integral corresponds to averaging over the surface of the sphere at contact. Similarly, for the lateral pressure we find

\displaystyle P_{xx}

\displaystyle=nk_{B}T+k_{B}T\frac{\pi n^{2}}{2}\int_{-\sigma}^{\sigma}g(\sqrt{\sigma^{2}-z^{2}},z)(\sigma^{2}-z^{2})\textrm{d}z.

(25)

Both pressures can be handled numerically using Simpson’s rule with the pair-distribution function obtained within PY. The transverse component of the pressure is rather sensitive to the $z$ -discretization of the $g$ calculation, so we use a finer grid spacing of $0.001\sigma$ for $\varphi=0.183$ and $0.002\sigma$ for $\varphi=0.37$ .

In our simulations, the components of the pressure tensor are calculated from the collisions of the spheres Henderson and van Swol (1984) as

\displaystyle P_{zz}

\displaystyle=nk_{B}T-\frac{m}{AL}\left\langle\frac{1}{t_{m}}\sum_{c}\frac{z_{ij}^{2}}{\sigma^{2}}(\mathbf{r}_{ij}\cdot\mathbf{v}_{ij}+z_{ij}v_{zij})\right\rangle

(26)

for the transverse component. The sum runs over all collisions $c$ occurring in a time interval $t_{m}$ . The indices $i$ and $j$ refer to the particles involved in the collision, with the difference in their $x$ – $y$ position at the instant of collision denoted by $\mathbf{r}_{ij}$ , the difference in the $z$ -component of the positions by $z_{ij}$ , the difference in their $x$ – $y$ velocity by $\mathbf{v}_{ij}$ , and the difference in the $z$ component of their velocity by $v_{zij}$ . The lateral component $P_{xx}$ is calculated similarly. Both components oscillate as a function of $L$ , plotted in figure 6. The confinement has a larger effect on the transverse component of the pressure, evident by the larger oscillations, seen for both the simulation and PY result. At large values of $L$ both the lateral and transverse components approach the same constant value, consistent with the fluid approaching the bulk state. The simulation value is approximately $0.3\%$ above the one calculated with the Carnahan-Starling equation of state Carnahan and Starling (1969) for a bulk hard-sphere fluid, which is within its expected accuracy Hansen and McDonald (2013). The PY value, for both packing fractions considered, does not quantitatively match the simulation result. However, qualitatively the curves look similar; the oscillations in each component occur at the same $L$ values, for both packing fractions.

Notably, the oscillations in the two components of the pressure tensor are out of phase. As the lateral pressure spikes, it is partially compensated for by a dip in the transverse pressure. The location of the peaks is dependent on the packing fraction, but the period of the oscillations is in both cases equal to the particle diameter. This matches the oscillations seen in the static structure factor, and indicates a dependence on commensurate or incommensurate packing.

IV Conclusions

We have investigated the effect of confining a fluid by applying extremely small periodic boundary conditions in one direction. Through the use of event-driven simulations and integral equation theory within the Percus-Yevick closure we have characterized the structure of the fluid through the generalized static structure factor. We see that the correlations in the structure show a non-monotonic dependence on the confining length scale, in a very similar manner to a liquid confined between two plates, even though the density profile of our fluid is constant. This indicates that the correlations in confined liquids may not be intrinsically related to the oscillating density profile they exhibit.

Additionally, we have seen that the non-monotonic behavior of the fluid can be observed in the pair-distribution function, and the components of the pressure tensor. As the confinement length is increased, the bulk isotropic fluid is approached, indicated by both components of the pressure tensor approaching the same value, and all elements of the structure factor overlapping.

While the Percus-Yevick approximation gives qualitatively the same results as the simulations, there are significant discrepancies between the structure factors calculated. One may seek to improve the solution empirically with the simulation data, in a manner similar to the Verlet-Weis correction Verlet and Weis (1972) used in the bulk. However, this method relies on the Carnahan-Starling equation of state, for which we have no equivalent in the quasi-confined fluid. As such, deriving the viral expansion of the equation of state in this system could be of value. Based on our calculations of the pressure, it will clearly be a complicated function of confinement length.

The present findings motivate further investigations into quasi-confined fluids. In particular, the dynamics close to the glass transition are of interest, as we expect this system provides a means to investigate the interplay of relevant length scales and incommensurability effects on the glass transition. Such an investigation is in a similar vein to investigations into the glass transition that considered the effects of curvature on local particle packing, by confining liquids to the surface of a sphere Vest et al. (2015, 2014, 2018) or hypersphere Turci et al. (2017). There, non-monotonic behavior of the correlations in the structure was observed for particles confined to the surface of a hypersphere with a circumference between 10 and 20 particle diameter lengths. A further extension to our work could include confining fluids to hyperspheres of smaller size. This would have the advantage of even higher symmetry than the periodic quasi-confinement, and could reveal an interplay between the effects of quasi-confinement and curvature. In order to study fluids at a higher packing fraction, in the vicinity of the glass transition, it would be necessary to use polydisperse hard spheres in the simulations. We expect that mode-coupling theory dynamics of the monodisperse hard spheres could be compared to simulations of slightly polydisperse systems van Megen and Underwood (1994), as is done successfully in the bulk Weysser et al. (2010). The results presented here have demonstrated that the non-monotonic behavior of pair correlations in confined liquids is not simply a reflection of the layering. This extension would seek to disentangle the non-monotonic behavior of the glass transition from the inhomogeneous density profile of confined fluids.

Acknowledgements.

We thank Martin Oettel for introducing us to the numerics of integral equation theory. CFP is supported by the Austrian Science Fund (FWF): M 2471. TF and LS are supported by the FWF: I 2887. The computational results presented have been achieved in part using the HPC infrastructure LEO of the University of Innsbruck.

References

Beebe et al. (2002) D. J. Beebe, G. A. Mensing, and G. M. Walker, Annu. Rev. Biomed. Eng. 4, 261 (2002).
Urbakh et al. (2004) M. Urbakh, J. Klafter, D. Gourdon, and J. Israelachvili, Nature 430, 525 (2004).
Pinilla et al. (2005) C. Pinilla, M. G. Del Pópolo, R. M. Lynden-Bell, and J. Kohanoff, J. Phys. Chem. B 109, 17922 (2005).
Mittal et al. (2008) J. Mittal, T. M. Truskett, J. R. Errington, and G. Hummer, Phys. Rev. Lett. 100, 145901 (2008).
Varnik and Franosch (2016) F. Varnik and T. Franosch, J. Phys.: Condens. Matter 28, 133001 (2016).
Lang et al. (2010) S. Lang, V. Boţan, M. Oettel, D. Hajnal, T. Franosch, and R. Schilling, Phys. Rev. Lett. 105, 125701 (2010).
Krakoviack (2005a) V. Krakoviack, J. Phys.: Condens. Matter 17, S3565 (2005a).
Krakoviack (2005b) V. Krakoviack, Phys. Rev. Lett. 94, 065703 (2005b).
Schmidt and Löwen (1996) M. Schmidt and H. Löwen, Phys. Rev. Lett. 76, 4552 (1996).
Schmidt and Löwen (1997) M. Schmidt and H. Löwen, Phys. Rev. E 55, 7228 (1997).
Hansen and McDonald (2013) J.-P. Hansen and I. R. McDonald, Theory of simple liquids: with applications to soft matter (Academic Press, 2013).
Bunk et al. (2007) O. Bunk, A. Diaz, F. Pfeiffer, C. David, C. Padeste, H. Keymeulen, P. R. Willmott, B. D. Patterson, B. Schmitt, D. K. Satapathy, J. F. van der Veen, H. Guo, and G. H. Wegdam, Phys. Rev. E 75, 021501 (2007).
Goel et al. (2008) G. Goel, W. P. Krekelberg, J. R. Errington, and T. M. Truskett, Phys. Rev. Lett. 100, 106001 (2008).
Goel et al. (2009) G. Goel, W. P. Krekelberg, M. J. Pond, J. Mittal, V. K. Shen, J. R. Errington, and T. M. Truskett, J. Stat. Mech. 2009, P04006 (2009).
Bollinger et al. (2015) J. A. Bollinger, A. Jain, J. Carmer, and T. M. Truskett, J. Chem. Phys. 142, 161102 (2015).
Nugent et al. (2007) C. R. Nugent, K. V. Edmond, H. N. Patel, and E. R. Weeks, Phys. Rev. Lett. 99, 025702 (2007).
Nygård et al. (2012) K. Nygård, R. Kjellander, S. Sarman, S. Chodankar, E. Perret, J. F. Buitenhuis, and J. van der Veen, Phys. Rev. Lett. 108, 037802 (2012).
Nygård et al. (2013) K. Nygård, S. Sarman, and R. Kjellander, J. Chem. Phys. 139, 164701 (2013).
Nygård et al. (2016) K. Nygård, J. Buitenhuis, M. Kagias, K. Jefimovs, F. Zontone, and Y. Chushkin, Phys. Rev. Lett. 116, 167801 (2016).
Nygård et al. (2017) K. Nygård, J. Buitenhuis, M. Kagias, K. Jefimovs, F. Zontone, and Y. Chushkin, Phys. Rev. E 95, 062601 (2017).
Nygård (2017) K. Nygård, Phys. Chem. Chem. Phys. 19, 23632 (2017).
Deb et al. (2011) D. Deb, A. Winkler, M. H. Yamani, M. Oettel, P. Virnau, and K. Binder, J. Chem. Phys. 134, 214706 (2011).
Mandal et al. (2014) S. Mandal, S. Lang, M. Gross, M. Oettel, D. Raabe, T. Franosch, and F. Varnik, Nat. Commun. 5, 4435 (2014).
Mandal et al. (2017) S. Mandal, S. Lang, V. Boţan, and T. Franosch, Soft Matter 13, 6167 (2017).
Mandal and Franosch (2017) S. Mandal and T. Franosch, Phys. Rev. Lett. 118, 065901 (2017).
Geigenfeind et al. (2015) T. Geigenfeind, S. Rosenzweig, M. Schmidt, and D. de Las Heras, J. Chem. Phys 142, 174701 (2015).
Scheidler et al. (2004) P. Scheidler, W. Kob, and K. Binder, J. Phys. Chem. B 108, 6673 (2004).
Gokhale et al. (2016) S. Gokhale, R. Ganapathy, K. H. Nagamanasa, and A. K. Sood, Phys. Rev. Lett. 116, 068305 (2016).
Eral et al. (2009) H. B. Eral, D. van den Ende, F. Mugele, and M. H. G. Duits, Phys. Rev. E 80, 061403 (2009).
Lang et al. (2012) S. Lang, R. Schilling, V. Krakoviack, and T. Franosch, Phys. Rev. E 86, 021502 (2012).
Lang et al. (2013) S. Lang, R. Schilling, and T. Franosch, J. Stat. Mech. 2013, P12007 (2013).
Lang et al. (2014a) S. Lang, R. Schilling, and T. Franosch, Phys. Rev. E 90, 062126 (2014a).
Ingebrigtsen et al. (2013) T. S. Ingebrigtsen, J. R. Errington, T. M. Truskett, and J. C. Dyre, Phys. Rev. Lett. 111, 235901 (2013).
Götze (2009) W. Götze, Complex dynamics of glass-forming liquids: A mode-coupling theory, Vol. 143 (OUP Oxford, 2009).
Percus and Yevick (1958) J. K. Percus and G. J. Yevick, Phys. Rev. 110, 1 (1958).
Thiele (1963) E. Thiele, J. Chem. Phys. 39, 474 (1963).
Wertheim (1963) M. S. Wertheim, Phys. Rev. Lett. 10, 321 (1963).
Freasier and Isbister (1981) B. C. Freasier and D. J. Isbister, Mol. Phys. 42, 927 (1981).
Leutheusser (1984) E. Leutheusser, Phys. A 127, 667 (1984).
Wertheim (1964) M. S. Wertheim, J. Math. Phys. 5, 643 (1964).
González et al. (1990) D. J. González, L. E. González, and M. Silbert, Phys. Chem. Liq. 22, 95 (1990).
Robles et al. (2004) M. Robles, M. López de Haro, and A. Santos, J. Chem. Phys. 120, 9113 (2004).
Lado (1968) F. Lado, J. Chem. Phys. 49, 3092 (1968).
Leutheusser (1986) E. Leutheusser, J. Chem. Phys. 84, 1050 (1986).
Adda-Bedia et al. (2008) M. Adda-Bedia, E. Katzav, and D. Vella, J. Chem. Phys. 128, 184508 (2008).
Sokołowski (1980) S. Sokołowski, J. Chem. Phys. 73, 3507 (1980).
Kjellander and Sarman (1991) R. Kjellander and S. Sarman, J. Chem. Soc., Faraday Trans. 87, 1869 (1991).
Boţan et al. (2009) V. Boţan, F. Pesth, T. Schilling, and M. Oettel, Phys. Rev. E 79, 061402 (2009).
Bénichou et al. (2013) O. Bénichou, A. Bodrova, D. Chakraborty, P. Illien, A. Law, C. Mejía-Monasterio, G. Oshanin, and R. Voituriez, Phys. Rev. Lett. 111, 260601 (2013).
Bénichou et al. (2016) O. Bénichou, P. Illien, G. Oshanin, A. Sarracino, and R. Voituriez, Phys. Rev. E 93, 032128 (2016).
Grossmann et al. (1996) S. Grossmann, D. Lohse, and A. Reeh, Phys. Rev. Lett. 77, 5369 (1996).
Meneguzzi et al. (1996) M. Meneguzzi, H. Politano, A. Pouquet, and M. Zolver, J. Comput. Phys. 123, 32 (1996).
De Lillo and Eckhardt (2007) F. De Lillo and B. Eckhardt, Phys. Rev. E 76, 016301 (2007).
Laval et al. (2001) J. Laval, B. Dubrulle, and S. Nazarenko, Phys. Fluids 13, 1995 (2001).
Karmakar et al. (2009) S. Karmakar, C. Dasgupta, and S. Sastry, Proc. Natl. Acad. Sci. 106, 3675 (2009).
Berthier et al. (2012) L. Berthier, G. Biroli, D. Coslovich, W. Kob, and C. Toninelli, Phys. Rev. E 86, 031502 (2012).
Puosi et al. (2016) F. Puosi, D. L. Cardozo, S. Ciliberto, and P. C. Holdsworth, Phys. Rev. E 94, 040102 (2016).
Denton and Egelstaff (1997) A. Denton and P. Egelstaff, Z. Phys. B 103, 343 (1997).
Pratt and Haan (1981a) L. R. Pratt and S. W. Haan, J. Chem. Phys. 74, 1864 (1981a).
Pratt and Haan (1981b) L. R. Pratt and S. W. Haan, J. Chem. Phys. 74, 1873 (1981b).
González-Melchor et al. (2005) M. González-Melchor, P. Orea, J. López-Lemus, F. Bresme, and J. Alejandre, J. Chem. Phys. 122, 094503 (2005).
Bernardi et al. (2015) S. Bernardi, S. J. Brookes, and D. J. Searles, Chem. Eng. Sci. 121, 236 (2015).
Petravic and Evans (1998) J. Petravic and D. J. Evans, Int. J. Thermophys. 19, 1049 (1998).
Roth (2010) R. Roth, J. Phys.: Condens. Matter 22, 063102 (2010).
Bollinger et al. (2014) J. A. Bollinger, A. Jain, and T. M. Truskett, Langmuir 30, 8247 (2014).
Saw and Dasgupta (2016) S. Saw and C. Dasgupta, J. Chem. Phys. 145, 054707 (2016).
Herrera-Velarde and Castañeda-Priego (2007) S. Herrera-Velarde and R. Castañeda-Priego, J. Phys.: Condens. Matter 19, 226215 (2007).
Evers et al. (2013) F. Evers, R. D. L. Hanes, C. Zunke, R. F. Capellmann, J. Bewerunge, C. Dalle-Ferrier, M. C. Jenkins, I. Ladadwa, A. Heuer, R. Castañeda-Priego, and S. U. Egelhaaf, Eur. Phys. J. Spec. Top. 222, 2995 (2013).
Capellmann et al. (2018) R. F. Capellmann, A. Khisameeva, F. Platten, and S. U. Egelhaaf, J. Chem. Phys. 148, 114903 (2018).
Lang et al. (2014b) S. Lang, T. Franosch, and R. Schilling, J. Chem. Phys. 140, 104506 (2014b).
Henderson (1992) D. Henderson, Fundamentals of inhomogeneous fluids (Dekker, New York, 1992).
Talman (1978) J. D. Talman, J. Comput. Phys. 29, 35 (1978).
Hamilton (2000) A. J. S. Hamilton, Mon. Not. R. Astron. Soc. 312, 257 (2000).
Rapaport (2004) D. C. Rapaport, The art of molecular dynamics simulation (Cambridge university press, 2004).
De Miguel and Jackson (2006) E. De Miguel and G. Jackson, Mol. Phys. 104, 3717 (2006).
Irving and Kirkwood (1950) J. Irving and J. G. Kirkwood, J. Chem. Phys. 18, 817 (1950).
Henderson and van Swol (1984) J. R. Henderson and F. van Swol, Mol. Phys. 51, 991 (1984).
Carnahan and Starling (1969) N. F. Carnahan and K. E. Starling, J. Chem. Phys. 51, 635 (1969).
Verlet and Weis (1972) L. Verlet and J.-J. Weis, Phys. Rev. A 5, 939 (1972).
Vest et al. (2015) J.-P. Vest, G. Tarjus, and P. Viot, J. Chem. Phys. 143, 084505 (2015).
Vest et al. (2014) J.-P. Vest, G. Tarjus, and P. Viot, Mol. Phys. 112, 1330 (2014).
Vest et al. (2018) J.-P. Vest, G. Tarjus, and P. Viot, J. Chem. Phys. 148, 164501 (2018).
Turci et al. (2017) F. Turci, G. Tarjus, and C. P. Royall, Phys. Rev. Lett. 118, 215501 (2017).
van Megen and Underwood (1994) W. van Megen and S. M. Underwood, Phys. Rev. E 49, 4206 (1994).
Weysser et al. (2010) F. Weysser, A. M. Puertas, M. Fuchs, and Th. Voigtmann, Phys. Rev. E 82, 011504 (2010).