Diffusion in Quasi-One Dimensional Channels: A Small System $n,p,T,$ Transition State Theory for Hopping Times

TL;DR

This paper develops a transition state theory to calculate hopping times for particles in quasi-one-dimensional channels, bridging single file and Fickian diffusion regimes, validated by simulations but limited by interaction neglect in wider channels.

Contribution

It introduces a transition state theory approach using the small system isobaric-isothermal ensemble to accurately compute hopping times in confined particle systems.

Findings

Theory accurately predicts hopping times for ideal gases and hard discs in narrow channels.

Simulation scheme effectively calculates free energy barriers for particle hopping.

Method underestimates barriers in wider channels due to interaction neglect.

Abstract

Particles confined to a single file, in a narrow quasi-one dimensional channel, exhibit a dynamic crossover from single file diffusion to Fickian diffusion as the channel radius increases and the particles can begin to pass each other. The long time diffusion coefficient for a system in the crossover regime can be described in terms of a hopping time, which measures the time it takes for a particle to escape the cage formed by its neighbours. In this paper, we develop a transition state theory approach to the calculation of the hopping time, using the small system isobaric--isothermal ensemble to rigorously account for the volume fluctuations associated with the size of the cage. We also describe a Monte Carlo simulation scheme that can be used to calculate the free energy barrier for particle hopping. The theory and simulation method correctly predict the hopping times for a two dimensional confined ideal gas system and a system of confined hard discs over a range of channel radii, but the method breaks down for wide channels in the hard discs case, underestimating the height of the hopping barrier due to the neglect of interactions between the small system and its surroundings.

Tables

Table 1. Table 1: Selected channel radii and system densities for the idea gas and hard sphere systems at β ​ P = 1.0 𝛽 𝑃 1.0 \beta P=1.0 .

Systems	$R_p$	$\rho$
2D-ideal gas	1.1-1.2	1.00
2D-Hard discs	1.01-1.2	0.38-0.41

Equations

$$⟨(x_t-x_0{)}^2⟩=⟨\Delta x_t^2⟩=2D_{x}t,$$

$$⟨\Delta x_t^2⟩=2F_x{t}^{1/2},$$

$$D_x\sim 1/{\tau }_{hop}^{\gamma }.$$

$$\Delta (N,p,T)=(V_o{)}^{-1}{\int }_{V}Q(N,V,T)e^{-pV/kT}\text{mbox},$$

$$Q_{n,v}^{\ast }dv=\frac{dv}{(n-1)!{\Lambda }^{Dn}}{\int }_{v+dv}{e}^{-\beta (U_n+U_{\sigma })}(dr{)}^{n-1}\text{mbox},$$

$$Q(N,V,T)={\int }_0^{V}Q(N-n,V-v,T){⟨Q_{n,v}^{\ast }⟩}_{0}dv\text{mbox},$$

$${⟨Q_{n,v}^{\ast }⟩}_0=\frac{{\int }_{V-v}{Q}_{n,v}^{\ast }{e}^{-\beta U_{N-n}}(dr{)}^{N-n}}{{\int }_{V-v}{e}^{-\beta U_{N-n}}(dr{)}^{N-n}}\text{mbox},$$

$$\begin{array}{rl}{P}_n(v)dv& =\frac{Q(N-n,V-v,T){⟨Q_{n,v}^{\ast }⟩}_{0}dv}{Q(N,V,T)}\\ & =e^{\beta n\mu }{e}^{-\beta W(v)}{⟨Q_{n,v}^{\ast }⟩}_{0}dv\text{mbox},\end{array}$$

$$\frac{Q(N-n,V-v,T)}{Q(N,V,T)}=e^{\beta n\mu }{e}^{-\beta W(v)}\text{mbox}.$$

$$\Delta =\int {⟨Q_{n,v}^{\ast }⟩}_0{e}^{-\beta W(v)}dv\text{mbox},$$

$$\Delta =\int Q_{n,v}^{\ast }{e}^{-\beta pv}dv\text{mbox}.$$

$$G=⟨U^{\ast }⟩+p⟨v⟩-TS\text{mbox},$$

$$G=-kT\ln \Delta \text{mbox}.$$

$${\tau }_{hop}=\frac{\kappa }{P^{\ast }}=\kappa e^{\beta \Delta G^{\ast }}\text{mbox},$$

$$dv=2R_{p}dL\text{mbox}.$$

$$\Delta =\int Q_{n,v}^{\ast }{e}^{-\beta p_{l}2R_{p}L}2R_{p}dL\text{mbox},$$

$${\Delta }_{x_c^{\prime }}d{x}_c={\int }_L{Q}_{n,v,T}^{\ast }{e}^{-\beta p_{l}2R_{p}L}\delta (x_c^{\prime }-x_c)2R_{p}dLd{x}_c\text{mbox},$$

$$\Delta ={\int }_{x_c}{\Delta }_{x_c}d{x}_c\text{mbox}.$$

$$P(x_c^{\prime })d{x}_c=\frac{{\Delta }_{x_c^{\prime }}d{x}_c}{\Delta }\text{mbox},$$

$${\int }_0^{\infty }P(x_c)d{x}_c=1\text{mbox}.$$

$$\beta \Delta G^{\ast }=-\ln P^{\ast }$$

$$P^{\ast }={\int }_0^{x_c^{\ast }}P(x_c)d{x}_c\text{mbox},$$

$$\begin{array}{rl}\Delta =& \frac{1}{(N-1)!{\Lambda }^{2n}}\times \\ & \int \int \int \int e^{-\beta U_n}{e}^{-\beta p_{l}2R_{p}L}dLd{y}_{1}d{x}^{n-1}d{y}^{n-1}\text{mbox},\end{array}$$

$$\begin{array}{rl}\Delta =& \frac{1}{(N-1)!{\Lambda }^{2n}}\times \\ & \int \int \int \int L^{n-1}{e}^{-\beta U_n}{e}^{-\beta p_{l}2R_{p}L}dLd{y}_{1}d{\bar{x}}^{n-1}d{y}^{n-1}\text{mbox},\end{array}$$

$$P({\bar{x}}^n,y^n;L)dL\propto L^{n-1}{e}^{-\beta U_n}{e}^{-\beta p_{l}2R_{p}L}dL\text{mbox},$$

$$\begin{array}{rl}\text{mbox}acc& \text{mbox}(old\to \text{mbox}new)=\text{mbox}min(1,\exp \{-\beta [U_n^{\text{mbox}new}-U_n^{\text{mbox}old}]\\ & -\beta p_{l}2R_p[L^{\text{mbox}new}-L^{\text{mbox}old}]+(n-1)\ln [L^{\text{mbox}new}/L^{\text{mbox}old}]\})\text{mbox}.\end{array}$$

$${\Delta }_{x_c}d{x}_c=\frac{d{x}_c}{{\Lambda }^4}{\int }_{-R_p}^{+R_p}{\int }_{-R_p}^{+R_p}{\int }_{x_c}^{\infty }{e}^{-\beta p_{l}2R_{p}L}dLd{y}_{1}d{y}_2$$

$${\Delta }_{x_c}d{x}_c=\frac{2R_p}{{\Lambda }^4}\frac{e^{-\beta p_{l}2R_p{x}_c}}{\beta p_l}d{x}_c.$$

$$\begin{array}{rl}\Delta & =\frac{2R}{\beta p_l{\Lambda }^4}{\int }_0^{\infty }{e}^{-\beta p_{l}2R_p{x}_c}d{x}_c\\ & =\frac{1}{{\beta }^2{p}_l^2{\Lambda }^4}\text{mbox}.\end{array}$$

$$V=2R_p⟨L⟩=-\frac{1}{\beta }{\left(\frac{\partial \ln \Delta }{\partial p_l}\right)}_{n,T}=\frac{2}{\beta p_l}\text{mbox},$$

Full text

I Introduction

Particle motion in bulk systems is described by Fickian (normal) diffusion, where the mean square displacement (MSD) of a labeled particle increases linearly with time in the long time limit, and is given by the Einstein relation,

[TABLE]

where $D_{x}$ is the diffusion coefficient along an arbitrary direction $x$, $x_{t}$ is the position of the particle at time, $t$, and $x_{0}$ is the initial particle position. However, the fundamental nature of diffusion can change dramatically if the movement of the particles is restricted to a quasi-one dimensional narrow channel in which the particles cannot pass each other Gelb et al. (1999); Gubbins et al. (2011). The geometric restriction, when combined with the presence of independent stochastic forces Levitt (1973); Mon and Percus (2003) or a Brownian background Percus (1974) acting on the particles, results in a form of anomalous diffusion known as single file diffusion (SFD), because the relative order of the particles is preserved and a particle must wait for its neighbors to move before being able to diffuse. This effect reduces the translational mobility of the particles in a way that causes the MSD to increase proportionally to the square root of time in the long time limit and is described by an Einstein–like relation,

[TABLE]

where the diffusion coefficient has been replaced by a mobility factor, $F_{x}$ Hahn and Kärger (1998); Lin et al. (2005).

Single file diffusion occurs in a variety of natural and engineered systems, including ion and water transport through biological membranes Hodgkins and Keynes (1955); Finkelstein (1987); Jensen et al. (2002); Rasaiah, Garde, and Hummer (2008), diffusion in molecular sieves such as zeolites Kärger et al. (1992); Kärger (2015); Kumar (2014), carbon nanotubes Mukherjee et al. (2010), or metal organic frameworks Salles et al. (2011); Jobic (2016), and charge-carrier diffusion in one dimensional channels Kharkyanen, Yesylevskyy, and Berezetskaya (2010). At the molecular level, pulsed field gradient nuclear magnetic resonance (NMR) techniques have been used to study SFD in adsorbate molecules confined in zeolites Gupta et al. (1995); Kukla et al. (1996); Hahn, Kärger, and Kukla (1996) and gas mixtures confined to polycrystalline dipeptide channels Dutta et al. (2016). However, the observation of SFD by NMR is subject to some uncertainty in systems where alternative diffusion mechanisms could be responsible Das et al. (2010); Valiullin and Kärger (2010). Single file diffusion has also been observed in colloidal systems Wei, Bechinger, and Leiderer (2000); Lutz, Kollmann, and Bechinger (2004), where it is being used in nano- and micro-fluidic devices Siems et al. (2012); Locatelli et al. (2016), as well as to control flow rates in the delivery of drugs Yang et al. (2010).

Fluids in quasi-one dimensional channels also exhibit an interesting dynamical crossover regime from SFD to normal diffusion as the channel radius increases beyond a passing threshold that allows the particles to hop past each other Mon and Percus (2002); Chen et al. (2010); Sané, Padding, and Louis (2009); Lucena et al. (2012); Herrera-Velarde, Pérez-Angel, and Castañeda-Priego (2016). In this crossover regime, the MSD in the long time limit is described by Eq. 1 because the particles can pass each other. However, the particles are trapped between their neighbors for long periods of time, during which they perform SFD, because passing is difficult. The resulting MSD at intermediate times is then described by Eq. 2. In the case of mixtures formed from particles of different sizes, it is also possible to select a channel radius that induces dual–mode diffusion, where the large particle component performs SFD while the smaller components exhibit normal diffusion Ball et al. (2009); Liu et al. (2010). The difference in particle mobilities can potentially be used for separation Adhangale and Keffer (2003); Wanasundara, Spiteri, and Bowles (2012).

Mon and Percus Mon and Percus (2002) used a simple phenomenological model to show that the long time diffusion coefficient of a system in the crossover regime can be written,

[TABLE]

where $\tau_{hop}$ is the average time for a particle to escape the cage formed by its two nearest neighbors and the exponent $\gamma$ depends on the nature of the hopping event. If the system is near the passing threshold for the particles, so that hopping events are rare and the particles are trapped long enough to perform SFD, then $\gamma$ is predicted to be $1/2$, and this has been confirmed by simulation for quasi-one dimensional hard discs Mon and Percus (2006, 2007), hard spheres Mon and Percus (2002) and soft-potential systems such as Lennard–Jones Wanasundara, Spiteri, and Bowles (2014), as well as for mixtures Ball et al. (2009). The hopping time approach is useful because the important details that influence the long time dynamics of the particles, such as the particle–particle and particle–wall interactions, as well as the density or pressure, are captured in a single, short time, phenomenological parameter that is also accessible by simulation and theory.

The process by which particles pass each other in a narrow quasi-one dimensional system can be treated in terms of a free energy barrier crossing process because the particle–particle and particle–wall interactions combine to increase the energy of interaction, and restrict the available configuration space for two particles, as they attempt to pass. This suggests that transition state theory (TST) could provide a useful description of the hopping time. In the case of hard sphere systems, $\tau_{hop}$ diverges as the radius of the channel decreases towards the passing threshold from above because the infinite nature of the interaction potential eventually prevents the particles from ever passing and TST correctly predicts the exponent for the power law behavior Bowles, Mon, and Percus (2004). Transition state theory also provides a qualitative description of $\tau_{hop}$ for soft potential systems, where the passing threshold is not as clearly defined Wanasundara, Spiteri, and Bowles (2014). However, these approaches are not quantitative. It is also worth noting that TST has been used to study diffusion in zeolite systems where the diffusing particles move between large structural cages formed in the channel by hoping through narrow windows Chmelik and Kärger (2016). In these systems, the cage volumes remained fixed and traditional TST is effective. However, in the case of structureless channels, such as the ones considered here, the particles are trapped by their neighbours so that the cages are dynamic and have fluctuating volumes. The goal of this paper is to develop a rigorous approach to the calculation of the hopping time, within the framework of transition state theory, that can provide the basis for the quantitative prediction of $\tau_{hop}$. To achieve this, we develop the theory and associated simulation method using the small system isobaric–isothermal ensemble to account for the fluctuating nature of the nearest neighbour cage that traps a particle in a single file fluid.

The remainder of the paper is organized as follows: Section II provides a brief summary of the background and general derivation of the small system $n,p,T$ ensemble in order to highlight the key elements of the method. This is followed by the description of how the small system $n,p,T$ ensemble can be used to develop a transition state theory for hopping times in confined single file fluids. Section III outlines the Monte Carlo (MC) simulation methods used to calculate the transition state free energy barriers associated with particle hopping and the direct measurement of the hopping times in a large system. Section IV describes the application of the theory and simulations for two simple confined systems: a two dimensional system of ideal gas particles confined between hard walls and a two dimensional system of hard discs confined between hard walls. We have chosen to focus on these two dimensional systems because they are analytically tractable, which allows us to check our simulation method, and our results can be compared to earlier hopping time studies performed in two dimensions. However, the two dimensional geometry is also relevant to the dynamics of colloidal particles confined to lithographically formed channels Wei, Bechinger, and Leiderer (2000), where the particles are effectively restricted to two dimensions by the gravitational potential well within the channel. Section V contains our discussion and our conclusions are summarized in Section VI.

II Background and Theory

II.1 The Small System $n,p,T$ Ensemble

In the thermodynamic limit, fluctuations away from the most probable microstates are small and the thermodynamics of the system can be described by the properties of the states associated with the maximum term in the probability distribution. The maximum term for one ensemble can be shown to be degenerate with the partition function of another so that all the ensembles are equivalent and selecting a particular ensemble for the calculation of a system property becomes a matter of mathematical convenience. However, as the system size decreases, fluctuations become more important and in the nano-scale regime it is essential to consider the types of fluctuations available to the system in selecting an appropriate ensemble. In the case of the isobaric-isothermal, $(N,P,T)$ ensemble, where the volume fluctuates, it also becomes necessary to consider the details of how the ensemble itself is formulated Brown (1958); Münster (1959); Sack (1959); Attard (1995); Koper and Reiss (1996); Corti and Soto-Campos (1998); Corti (2001, 2002). The original formulation of the ensemble by Guggenhiem Guggenheim (1939) can be expressed in integral form Hill (1987) as

[TABLE]

where $Q(N,V,T)$ is the canonical partition function for $N$ particles contained in a volume $V$ at temperature $T$, $k$ is Boltzmann’s constant and $p$ is the external pressure being applied to the system. The prefactor contains a volume scale, $V_{o}$, to ensure the partition function remains dimensionless. In the thermodynamic limit the choice of $V_{o}$ is arbitrary and it is usually taken to be $kT/p$, but for small systems, the value of the volume scale depends on the properties of the system such as the number of particles and the details of the boundary defining the volume.

To address these concerns Corti Corti (2001) developed a rigorous approach to the isobaric-isothermal ensemble for small systems that avoids the need to specify a volume scale. Here, we will provide a brief summary of Corti’s derivation for a small system, $n,p,T$ ensemble, highlighting the key elements of the method and how they pertain to our particular application. The analysis begins by considering the division of the canonical partition function for $N$ particles, in a volume $V$ and at a temperature $T$, into a small subsystem consisting of $n$ particles in a sub–volume $v$ that is located at a point $r_{0}$ and surrounded by a system of $N-n$ particles in the remaining volume, $V-v$ (see Fig. 1(a)). Separating the subsystem from its surroundings is a mathematical boundary that has no mass or momentum and is assumed to be spherical. To avoid an over counting of configurations caused by the fluctuation of the volume that does not involve a change in particle coordinates, it is necessary to tie the location of the boundary to a shell molecule. This involves requiring that at least one of the $n$ particles of the subsystem is located in a shell, $dv$, at the boundary. Now when the volume fluctuates, the shell molecule moves with the boundary to ensure the system samples a new configuration.

A canonical partition function for the $n$ particles in the volume $v+dv$ can be written as,

[TABLE]

where $\Lambda$ is the de Broglie wavelength, $D$ is the dimensionality, $U_{n}$ is the potential energy of the $n$ particles in the small volume interacting amongst themselves, $U_{\sigma}$ is the interaction between the particles in the subsystem and the $N-n$ particles in the surroundings and $(dr)^{n-1}=dr_{1,2}\cdots dr_{1,n}$ are the particle volume elements with positions relative to the location of the shell particle. Furthermore, the integration of the shell particle, over the volume in the shell, $dv$, and the other $n-1$ particles over $v+dv$, is performed with the positions of the surrounding $N-n$ particles held fixed so the canonical partition function for the complete $N,V,T$ system can be expressed,

[TABLE]

where,

[TABLE]

is an average over the configurations of the surrounding $N-n$ particles interacting amongst themselves with a potential energy $U_{N-n}$.

Now, the probability of finding the system with $n$ particles in a volume $v+dv$ and the remaining particles outside can be written,

[TABLE]

where $\mu$ is the chemical potential of the surroundings, $W(v)$ is the work required to form an empty cavity in the system and the following relation has been used:

[TABLE]

The small system $n,p,T$ partition function can then be written,

[TABLE]

by using the normalization condition, $\int_{v}P_{n}(v)dv=1$ and noting $\Delta=\exp(-\beta n\mu)$. Depending on the size and shape of the cavity, $W(v)$ may contain work in addition to the usual volume work, such as surface work, so that the pressure inside the system, $\hat{p}$, does not necessarily equal the imposed pressure. However, if the interface is planar, or the small subsystem becomes macroscopic, then $\hat{p}=p$ and $W(v)=pv$. If the interactions between the subsystem and its surroundings are also negligible, Eq. 10 becomes,

[TABLE]

Most notably, Eq. 11 contains no volume scale because the uncertainty associated with identifying new configurations in phase space as the volume of the container fluctuates is avoided by the use of the shell molecule. However, considerable care needs to be taken when implementing the new small system $N,p,T$. For example, the integral over the shell particle position performed to obtain Eq. 11 relies on the underlying spherical symmetry of the system volume. The same approach can be taken for systems with alternative shapes such as cubes etc, but this is not always possible when considering non-isotropic systems, as we will do in this work, and alternative approaches to evaluating degrees of freedom associated with the shell particle may be needed.

Corti Corti (2002) also showed that the small system isobaric-isothermal partition function is related to the thermodynamic potential,

[TABLE]

where $\left<\cdots\right>$ denotes the ensemble average and $U^{*}$ is the internal energy in the presence of a shell particle, through the expression

[TABLE]

It is important to note that $G$ does not necessarily equal the Gibbs free energy of the system, except in the thermodynamic limit, because it contains information regarding the surroundings, as is evident from Eq. 8 where the chemical potential of the surrounding is introduced, and because the small system samples a different set of volume fluctuations compared to a thermodynamically sized system. This is the case even when $U_{\sigma}=0$, and there is no interaction between the system and the surroundings. Consequently, a change in $G$, defined by Eq. 13, represents the reversible work of moving the system between two states, including contributions from the surroundings. A more complete discussion of these subtleties can be found in references Corti and Soto-Campos (1998); Corti (2001, 2002).

II.2 A Transition State Theory for Hopping Times

Transition state theory is an equilibrium based approach to the calculation of transition rates in activated processes, originally developed by Eyring Eyring (1935) to describe rates in chemical reactions. While the approach neglects the effect of barrier recrossing, it provides an upper bound to the rate and more sophisticated methods, such as those proposed by Bennet Bennett (1975) and Chandler Chandler (1978), reduce to TST in the regime where the barrier is high relative to the thermal energy in the system and recrossing is rare.

According to TST, the time it takes for a barrier crossing event like hopping to occur is inversely proportional to the rate and is given by,

[TABLE]

where $\kappa$ is a kinetic prefactor, $P^{*}$ is the probability of finding the system in the transition state and $\Delta G^{*}$ is the height of the free energy barrier for the activated process.

Figure 1(b) shows the construction of the small $n,p,T$ system used here to study the TST free energy barrier associated with a particle escaping its cage in a quasi-one dimensional channel by hopping past one of its neighboring particles. We consider an $n=2$ particle system confined to a long, uniform, two dimensional (2D) channel that extends longitudinally along the $x$-axis and has a radius $R_{p}$ in the $y$-axis. The cell center is located at a point $r_{0}$ and the shell particle is located at $+L/2$ so that the total length of the cell is $L$ and its two dimensional volume is $v=2R_{p}L$. The channel radius remains fixed so that fluctuations in the volume are given by,

[TABLE]

The small system isobaric–isothermal partition function given by Eq. 11 then becomes,

[TABLE]

where $p_{l}$ is the longitudinal pressure acting on the end of the channel Kofke and Post (1993). Particle 1 is the shell particle placed on the boundary at one end and the $2R_{p}$ factor in $dv$ arises from the integration of the shell particle over its $y$-axis degree of freedom.

The second particle in our system represents the caged particle. If we define a reaction coordinate for the hopping process as the axial separation of the caged particle from the shell particle, $x_{c}=x_{1}-x_{2}$, then the transition state occurs at $x_{c}=0$ as shown in Fig. 1(c). The configuration space associated with specific points along the reaction coordinate can be measured by applying a delta function, $\delta(x_{c}^{\prime}-x_{c})dx_{c}$, that is one when $x_{c}^{\prime}-x_{c}=0$ and zero otherwise, to Eq. 16, to obtain a reaction coordinate partition function,

[TABLE]

which satisfies,

[TABLE]

The probability of having the caged particle at $x_{c}^{\prime}$ is then,

[TABLE]

where $P(x_{c}^{\prime})$ is a probability density and

[TABLE]

Finally, the Gibbs free energy barrier for particle hopping, $\Delta G^{*}$, is directly related to the probability of finding the caged particle in the transition state by,

[TABLE]

where

[TABLE]

and the range of the reaction coordinate, $x_{c}=0$ to $x_{c}^{*}$, is the transition state region. Defined in this way, $\Delta G^{*}$ represents the work required to bring the caged particle into the transition state from anywhere in the cage and, as such, it will always be positive because the configuration space of the transition state is a restricted subset of all the possible configurations available to the system. Identifying the entire configuration space of the “reactive region” as the appropriate reference state leads to a significant improvement in the predicted rates using TST in activated process as shown in the case of heterogeneous nucleation Scheifele et al. (2013).

The boundary separating two neighbouring cages along the reaction coordinate occurs at $x_{c}=0$. However, it is necessary to define a transition region with $x_{c}^{*}$ to ensure the probability $P^{*}$ in Eq. 21 is dimensionless, rather than a probability density. We do not have a general method for selecting $x_{c}^{*}$, but we would expect the transition region to be small relative to size of the fundamental particle motion responsible for taking the particle from one cage to the other. If $x_{c}^{*}$ is chosen to be too large, then the transition state will contain contributions from configurations that have little or no chance of crossing the barrier in the small time limit, which would lead to poor estimates of the hopping time.

It is important to note that we have chosen to focus on a system of just two particles so that we have a single transition state to consider when a particle is attempting to hop. It could be argued that the cage of the central particle is actually formed by two particles, one on either side, giving rise to two transition states. However, each transition state involves two particles, and both of these particles escape their cages for each hopping event, which suggests we only need to consider one transition per caged particle. Explicitly including the third particle would also complicate the analysis as it would require the use of two shell molecules to define the volume. The reaction coordinate, which measures the distance of the cage particle from the transition state, would also take on two values for each configuration due to the symmetry of the system. Furthermore, the configurations of the particles outside the cell are formally included in the analysis when the full canonical partition function is divided into the small system to formal the small system $n,p,T$ ensemble. This include configurations where a particle from the surroundings sits on the boundary of the cell at $-L/2$, to form the other end of the cage.

III Simulation Methods

III.1 Free Energy Barrier Calculations

Corti Corti (2002) showed that the probability for accepting a trial MC move in the small system $n,p,T$ ensemble depends on the shape of the small system volume, which also constitutes the simulation cell. If we consider the case of $n$ particles confined to the two dimensional channel described in Fig. 1(b), and ignore interactions between the particles in the cell and those in the surroundings, then the partition function for the system can be written,

[TABLE]

where the integrals over the $dxdy$ coordinates of the $n-1$ particles within the cell represent $Q^{*}_{n,v,T}$, the integral over $dy_{1}$ represents the degrees of freedom of the shell particle over the diameter of the channel and the integral over $dL$ is the degrees of freedom of the shell particle associated with the volume of the cell. If a scale factor $x_{i}=L\bar{x}_{i}$ is introduced, then Eq 23 can be expressed as,

[TABLE]

where the $U_{n}$ is now a function of the scaled coordinate system. Note there is no need to rescale the $y$-coordinates because fluctuations in the volume only involve changes in the longitudinal dimension of the cell. The probability of finding the system in a configuration specified by $\bar{x}^{n},y^{n}$, in a volume with the cell length in the range $L$ to $L+dL$, is then,

[TABLE]

and the MC acceptance probability for a trial move from an old to a new configuration becomes,

[TABLE]

While our analysis is specific to the two dimensional channels studied here, extending Eq. 26 to the study of three dimensional, quasi-one dimensional channels, where fluctuations in the volume only involve changes in $L$, simply requires the appropriate geometric factor for the channel cross section to be included in the pressure-volume term.

The probability, $P(x_{c})dx$, of finding the caged particle in a volume $dx$ at a distance $x_{c}$ from the transition state can now be calculated directly using MC simulation. The small system $n,p,T$ simulations are carried out with $n=2$ particles confined to a channel with radius, $R_{p}$, and length, $L$. In the case of the hard disc model, the unit of length is taken relative to the diameter of a particle $\sigma=1$. Particle one is the shell particle that defines the sub-volume and is located at $L/2$ and the center of the sub-volume is located at the origin. Particle two is placed within the cell between $-L/2$ and $L/2$. Since the two models being studied are the two dimensional ideal gas and two dimensional hard discs, the MC moves proceed as follows: A particle is selected at random and moved randomly by a step $\delta x$ and $\delta y$, up to a maximum displacement of $|\Delta x|=|\Delta y|=0.05$. The move is immediately rejected if the trial displacement causes the particle to leave the channel or overlap with with another particle. If particle one was selected, the trial move also corresponds to a change in $L$ of $2\delta x$ and the position of particle two is rescaled to ensure it remains within the simulation cell. Eq. 26 is then used to determine the probability of accepting the move. It should also be noted that the position of the shell particle must remain positive so that $V>0$. For each state point studied, $2\times 10^{7}$ MC moves are used to obtain equilibrium and data is collected over the next the next $10^{8}-10^{9}$ moves, depending on the height of the free energy barrier. The probability is calculated by building a histogram of configurations along the reaction coordinate using bin sizes of $\Delta x_{c}=0.005$ and $0.03$ for the ideal gas and hard discs respectively, and sampling the system every 1000 MC moves to ensures configurations were not correlated. Simulation runs were sufficiently long to ensure the transition state region, $x_{c}<\Delta x_{c}$, was sampled at least 100 times and the probability density was obtained by dividing the probability by the bin size.

III.2 Hopping Time Calculations

The hopping time, which is the average time it takes for a particle to pass one of its caging nearest neighbors in the single file, is calculated using MC simulation in the canonical ensemble. The system consists of $N=100$ particles confined to a two dimensional channel of radius $R_{p}$ in the $y$–direction and length $L$ along the $x$–axis, where periodic boundaries are also employed. Particles are moved using the standard Metropolis MC scheme Frenkel and Smit (2002) as a simple approximation to Brownian motion Patti and Cuetos (2012). A MC trial move involves the random selection of a particle, followed by the random selection of a trial step in both the $x$ and $y$–directions up to a maximum step size of $|\Delta x|=|\Delta y|=0.05$. The trial move is rejected if it causes either a particle or wall overlap and the particle is returned to its original position, otherwise it is accepted. A unit of time is defined as a single MC cycle and consists of $N$ MC trial moves. The particles are initially spaced evenly along the tube and placed randomly across the channel such that they do not overlap then the system is equilibrated over $2\times 10^{7}$ MC cycles. At the start of data collection, the cages for the particles are identified as their immediate right and left neighbors, and their initial hopping time is set to zero. A particle hopping time is the number of MC cycles it takes for the particle to escape its cage. After each hopping event the hopping time for the particle is recorded, then reset to zero, and the new caging neighbors are identified. Data is collected over $5\times 10^{7}$ MC cycles and the average hopping time, $\tau_{hop}$, is calculated over all hopping times recorded for all particles.

The hopping times for systems with different radii must be compared at the same applied external pressure. However, the calculations of $\tau_{hop}$ are carried out in the $N,V,T$ ensemble because volume fluctuations in the $N,p,T$ ensemble cause the particle positions to be rescaled, leading to possible errors. To ensure $\tau_{hop}$ is obtained at the correct state point, the density for each system with a different radius is calculated using the $N,p,T$ ensemble where $\beta P=1.0$, $N=1000$, and $10^{8}$ MC cycles are used after the system reaches equilibrium. A MC trial move involves the random selection of a particle, followed by the random selection of a trial step in both the $x$ and $y$–directions. The maximum step size of $|\Delta x|=|\Delta y|$ is varied during the simulation in the range of $0.05-0.1\sigma$ to ensure $80\%$ of the trial moves are accepted. The shell particle approach is used. Data are sampled at every $10,000$ time steps to obtain the average volume, $\left<V\right>$, which gives the density $\rho=N/\left<V\right>$. Table 1 shows the range densities, $\rho=N/V$, obtained for the different radii studied in both the ideal gas and hard disc systems. The density range for hard disc two particle system is $\rho=0.54-0.65$, which is significantly higher than that obtained for the large system because the lack of interactions with the surrounding allows the system to small volume configurations that would be inaccessible if more particles were present. When a third particle is included in the small system $N,p,T$, the density range reduces to $\rho=0.47-0.53$. Nevertheless, as noted earlier, it is necessary to calculated $\tau_{hop}$ using the two particle system but this result highlights the need to ultimately account for the interactions between the cell and the surrounding, i.e. the inclusion of $U_{\sigma}$.

IV Applications

IV.1 Two Dimensional Ideal Gas

To illustrate our approach, we first calculate the hopping barrier in a two dimensional ideal gas system, where there are no interactions between particles, so that $U_{\sigma}=U_{n}=0$, and the hard interaction of the particles with walls keeps them confined between $-R_{p}$ and $R_{p}$. The reaction coordinate partition function for this system can be written,

[TABLE]

where the length of the cell is bounded in its lower limit by the position of the caged particle, and the integrals $y_{1}$ and $y_{2}$ can be performed independently to give,

[TABLE]

The partition function for the system can then be obtained by integrating over the remaining $x_{c}$ coordinate to yield,

[TABLE]

To demonstrate that we have the correct partition function, we calculate the system volume,

[TABLE]

which is the expected equation of state for an ideal gas Corti (2001).

Using Eqs. 28 and 29 in Eq. 19 yields the probability of finding the caged particle at a point $x_{c}$ along the reaction coordinate as,

[TABLE]

The probability of finding the caged particle within a region $dx_{c}$ of the transition state, where $x_{c}=0$, is then,

[TABLE]

Figure 2 shows that the free energy density, $-\ln P(x_{c})$, increases linearly as a function of $x_{c}$, with a slope of $\beta p_{l}2R_{p}$, which reflects the fact that increasing $x_{c}$ restricts the total number of configurations available to the system by increasing the lower bound on the accessible volume fluctuations. However, it also seems to imply that the hopping process for an ideal gas particle is “barrierless”, but this is not the case. The free energy barrier $\Delta G^{*}$, given by Eq. 21, is necessarily positive because it involves an integral of the normalized probability over a small region of the reaction coordinate in the transition region. This indicates that it takes work to restrict the two particles to the transition state region relative to having the caged particle located anywhere within the small system volume. Figure 2 also shows that the probability density obtained from our simulations matches the results obtained from our analytical analysis. While this is not surprising for such a simple problem, it provides a proof of principle and suggests that our simulation approach would be applicable to systems involving more complicated particle–particle and particle–wall interactions.

Figure 3 shows a plot of $\ln\tau_{hop}$, obtained from our hopping time simulations versus the free energy barrier to hopping obtained from the analytical results over a range of channel radii, where we have used Eqs. 31 and 21, with $x_{c}^{*}=0.005$. The slope is one, which suggests that our TST approach, using the small system isobaric–isothermal ensemble, provides a potentially useful method for the calculation of hopping times in these quasi-one dimensional systems. The intercept yields the kinetic prefactor $\kappa=0.90$ and essentially measures the time it takes for the particles to move out of the transition state region. Equation 32, with Eq. 14, suggests that $\tau_{hop}$ should vary linearly as $R_{p}^{-1}$, with a slope $\sim\kappa/2\beta p_{l}x_{c}^{*}$ (see Fig. 4). This dependence on channel radius arises because our analysis is performed in the small system $n,p,T$ ensemble, where it is necessary to integrate over all possible volume fluctuations (see Eq.28). It differs from the scaling law observed for hard disks in the narrow channel regime and provides a limiting case for the hopping time in wide channels where the particles can easily pass.

IV.2 Two Dimensional Hard Discs

We now calculate the hopping barrier for a system of two dimensional hard discs confined to a channel by hard walls, where the particle-particle exclusion interaction makes it more difficult for particles to escape their cage as the channel radius decreases. For the hard disc system the particle–particle interaction is defined as,

[TABLE]

where $\sigma$ is the particle diameter and $r_{ij}=\lvert\textup{r}_{i}-\textup{r}_{j}\rvert$ is the the distance between two particles. The particle–wall interaction is given by,

[TABLE]

where $y_{i}$ is the $y$–coordinate of the particle. To make the solution to the small system $n,p,T$ partition function tractable, we consider the case where the pressure is low. The interactions between the particles inside the sub-volume and those in the surroundings are then negligible and we can assume $U_{\sigma}=0$. The reaction coordinate partition function for this system is then given by,

[TABLE]

where,

[TABLE]

and the factor of 2 appears because we need to include the configurations in which the order of particles 1 and 2 in the $y$-coordinate are reversed. Taking into account the piecewise nature of $\sigma_{2}$ and integrating over $x_{c}$ yields,

[TABLE]

where $I_{n}(z)$ is the modified Bessel function of the first kind, and $L_{n}(z)$ is the modified Stuve function, both with $n=1$ and $z=\beta p_{l}\sigma(2R_{p}-\sigma)$ Abramowitz and Stegun (1965).

Figure 5 shows the free energy density, $-\ln P(x_{c})$, as a function of $x_{c}$ for channel radii in the range $R_{p}/\sigma=1.01-3$. For $x_{c}>\sigma$, the two discs do not exclude volume with respect to each other ($\sigma_{2}=0$) and they effectively behave as ideal gas particles in a channel with a reduced radius of $R_{p}-\sigma/2$, which gives rise to the linear increase in $-\ln P(x_{c})$. When $x_{c}<\sigma$, the two discs begin to exclude volume, leading to a restriction in configuration space. However, when the channel diameter is wide, the volume exclusion effect is small and the free energy density remains close to being linear. As the channel becomes narrower, we see the appearance of a maximum that begins at $x_{c}=1$ and moves towards $x_{c}=0$ as $R_{p}\rightarrow\sigma$. This unusual feature results from the convolution of two competing effects. The excluded volume interaction between the two discs reduces the number of accessible configurations available to the system as the transition state is approached. Conversely, the pressure–volume term increases the number of accessible configurations because the range of allowed volume fluctuations increases as $x_{c}\rightarrow 0$. The free energy density profile obtained here, in the small system isobaric-isothermal ensemble, differs significantly from that obtained in the canonical ensemble where the free energy maximum is always located at $x_{c}=0$ Wanasundara, Spiteri, and Bowles (2014).

In the context of particle motion, the transition state for hopping is still located at $x_{c}=0$, despite the presence of a maximum in the free energy density at values of $x_{c}>0$, because it marks the point where the particles exchange position, signifying the crossover from single–file to normal diffusion. The probability of finding the caged particle within a region $dx_{c}$ of the transition state can be obtained using Eqs. 35 and 37 in Eq. 19, and noting $\sigma_{2}=\sigma$ (Eq. 36) when $x_{c}=0$, to give,

[TABLE]

The hopping time in the hard disc model is expected to diverge as a power law, $\sim(R_{p}-\sigma)^{\alpha}$, as the channel radius approaches the diameter of a disc because the excluded volume interaction prevents the two particles from reaching the transition state. Once the channel radius becomes smaller than the passing threshold, $R_{p}=\sigma$, the discs become permanently caged between their neighbours. A Ficks–Jacobs analysis Zwanzig (1992) that projects the diffusion of the two discs onto a one-dimensional reaction coordinate predicts Kalinay and Percus (2005); Kalinay (2007) $\alpha=-3/2$, while TST predicts Bowles, Mon, and Percus (2004) $\alpha=-2$. Simulation studies of the hopping time agree with the TST result, but the scaling only becomes apparent for very narrow channels Mon (2008); Mon and Percus (2007). Transition state theory suggests $\tau_{hop}\sim 1/P(0)$ and Eq 38 shows that the hopping time should diverge with the exponent $\alpha=-2$, as expected. Figure 6 also shows how $1/P(0)$ crosses over from the passing threshold limit for narrow channels to the wide channel scaling behaviour $\sim(2R_{p}-\sigma)^{-1}$, which is consistent with our ideal gas results. Finally, Fig. 7 shows that Eq. 14 works well for the hard disc system when the channel radius is narrow, giving rise to high barriers. This is also the regime where the hopping time is directly connected to the diffusion coefficient through Eq. 3 because the particle must be trapped long enough to perform SFD between hops. However, our TST approach begins to break down for wider channels where we appear to underestimate the height of the barrier needed to predict the hopping time.

V Discussion

The phenomenological hopping time approach to understanding diffusion in single file fluids confined within quasi–one dimensional channels is attractive because all the effects that influence particle motion in the long time limit are contained within a single parameter, $\tau_{hop}$, that describes the short time event of a particle hopping past one of its neighbours. The activated nature of the hopping process also makes the calculation of $\tau_{hop}$ amenable to the application of transition state theory, but the challenge is to provide a rigorous approach for the calculation of the probability of being in the transition state, i.e. the height of the hopping free energy barrier, for a nanoscale system that exhibits volume fluctuations. While the isobaric–isothermal ensemble is the appropriate choice of ensemble, a direct application of the bulk thermodynamic results to a system in the nanoscale is problematic due to the system size dependence of the necessary volume scale.

We have shown that this problem can be overcome by using the small system isobaric–isothermal ensemble proposed by Corti, where the use of a shell particle to define the volume prevents the over counting of configurations when considering volume fluctuations and eliminates the need for an explicit volume scale. Figures 3 and 7 highlight the key results of our work and show that $\ln\tau_{hop}$ varies linearly with $\Delta G^{*}$ with the slope of unity for both systems studied here, over a range of different radii, as predicted by TST. One of the main advantages of our approach is that it only requires the simulation of two fluid particles and the wall interactions, while the direct measurement of the diffusion coefficient and $\tau_{hop}$ involves the simulation of hundreds of particles and requires a long simulation time because $\tau_{hop}$ is slow to converge. The small number of particles also means that the probability density along the reaction coordinate can be sampled by “brute force” because simulations with a large number of MC cycles can be performed quickly, but the method can also be adapted to include more advanced techniques such as biased umbrella sampling Frenkel and Smit (2002) in cases where the barrier is large.

Choosing the entire partition function of the caged particle as the appropriate reference state to define $\Delta G^{*}$ is critical to the implementation of our approach. As noted earlier this means that $\Delta G^{*}$ is necessarily positive because the configuration space of the transition state is smaller than that of the entire system, even for the case of the ideal gas where there is no volume exclusion by the particles as they approach each other. This, combined with the use of the small system $n,p,T$ ensemble, leads to the unexpected prediction that $\tau\sim 1/R_{p}$ which was confirmed through our simulations.

It is interesting to note that our TST approach breaks down for wider channels in the hard disc model, but works well over all channel radii in the ideal gas model even when the barrier is only a few $kT$. The main difference between the two models is the exclusion interaction of the hard discs and the effect this has on our decision to neglect the interactions between the small system and its surroundings. Setting $U_{\sigma}=0$ is exact for the ideal gas system. In the hard discs model, this allows the two particles in small system to sample configurations in the transition state that would otherwise be ruled out by the exclusion interaction with particles in the surroundings sitting near the boundary. As a result, our method over estimates the partition function of the transition state which leads to a lower predicted barrier. Neglecting these interactions becomes less important as the channel becomes narrower. This assumption also restricts the application of our method to low pressures, where again, interactions between the small system and the surroundings are minimal.

The hopping time is only connected to the diffusion coefficient through Eq. 3, with $\gamma=1/2$, when the particles are caged long enough between hops that they perform single file diffusion during this time, which should occur when the barrier is high. This suggests that our method should be effective in predicting the hopping time and dynamic properties of single file fluids in the crossover regimes and there are a wide range of important practical challenges in nano-fluidic systems, including the separation of mixtures, that can be addressed in the regime where our method is valid and simple enough to provide useful insight. However, Mon and Percus Mon and Percus (2002) also showed that Eq. 3 may be applicable in the wide channel regime, but with a different hopping time exponent. When the channel is wide, the particles can pass each other easily because the barrier is low. The hopping time is then short so that the distance traveled between hops becomes independent of the hopping time, and proportional to the average distance between nearest neighbors in the channel, which leads to an exponent of $\gamma=1$. This has not been tested, but suggests the hopping time approach may still be applicable even when the channel is wide. To make our TST approach work in the wide channel regime, it would be essential to account for the interactions between the small system and the surroundings in our simulation method.

VI Conclusions

In this work, we have shown that the small system isobaric–isothermal ensemble provides a rigorous framework for the development of a transition state theory for the calculation of the hopping time and we have applied it to two simple systems. While the simulation method is currently restricted to the narrow channel, low pressure regime, where interactions between the small system and its surroundings can be neglected, it is simple to implement and is potentially applicable to the study of a wide range of systems.

Acknowledgements.

We would like to thank Natural Sciences and Engineering Research Council of Canada (NSERC) and POS Bio-Sciences for financial support. Computing resources and support were provided by Compute Canada and WestGrid.