Interacting particles crossing barriers: spatial statistics, choice of route taken and role of interactions

Introduction

Escaping from a potential well and crossing a barrier is a remarkably important problem in physics and chemistry. It underpins the dynamics of atoms and molecules and processes ranging from Brownian particles in a potential well to enzyme catalysis to Josephson junctions and more [1–3].

While most studies focus on the dynamics of a single particle trapped in a potential well, it is also interesting to study the effect of interactions when multiple atoms or molecules are present and simultaneously attempting to escape the well, driven by thermal fluctuations. The dynamics of permanently bound particles (dimers) in a deterministic, anharmonic potential were investigated [4], and in the presence of noise [5], where a non-monotonic dependence of the escape rate on the binding constant was observed. Brownian motors (ratchets) have also been studied in the context of interparticle interactions, which, depending on the coupling, can even lead to a reversal of probability currents [6, 7]. The problem of long polymers confined in the metastable side of a double-well potential has also been addressed [8]. More recently, Singh studied coupled particles in one-dimensional potentials [9, 10].

In this work, we are interested in slightly more realistic potential energy landscapes, that can often be complex and rugged, featuring numerous saddle points, barriers of different heights, and multiple paths connecting stable states [11, 12], which arise in condensed matter systems, such as glassy systems, colloidal fluids, or gels. At the most basic level, a complex potential energy landscape exhibits competing pathways and metastable minima.

Specifically, we investigate how a pair of particles explore and cross the barriers in the energy landscape previously studied in Refs. [12, 13], where the dynamics of a single particle was investigated. We give all the necessary details of the potential below, but the key features of the potential are that it has a choice of two routes for the particles to reach the global minimum in the potential. One route is longer, but goes over a lower barrier (i.e. lower saddle point), while the other route is shorter, but is over a higher barrier. One of the main findings of Refs. [12, 13] is that if the time available for barrier crossing is short, then if barrier crossing is to occur, the most likely route from one minimum to the other is the path over the higher barrier. In contrast, if the time available is longer, then the most likely route is over the lower barrier. Here, we consider a pair of particles confined within this potential and how the particle interactions influence which particle arrives in the destination minimum first and the spatial statistics of the second particle, at the moment when the first arrives. We find that the nature of the particle interactions greatly influences these quantities, with order of magnitude changes in arrival times possible, depending on the nature of the interactions.

We investigate here the effect of an interaction potential $V(r)$ between two Brownian particles separated by a distance $r$ and navigating a background potential energy landscape $\phi(\mathbf{x})$. We compare two very different interparticle potentials $V(r)$. First, we consider the soft Gaussian core model $V_{\mathrm{G}}(r) \propto\exp[-r^2/\lambda_{\mathrm{G}}^2]$, with range parameter $\lambda_{\mathrm{G}}$, originally introduced by Stillinger in the context of crystallization of charged colloidal particles or polymer solutions [14], and is now widely used for complex fluids [15, 16]. Second, we consider a Lennard-Jones pair potential $V_{\mathrm{LJ}}(r)$.

In recent years, transition path times have garnered considerable attention [17–32] due to their ability to offer valuable information on structural transformations, such as those occurring during protein folding and phase transitions [33]. However, the majority of studies have predominantly focused on one-dimensional energy landscapes or two-dimensional systems involving a single particle, leaving the exploration of more complex, multiparticle systems largely under-explored.

The rest of the paper is organized as follows. Section II describes the two simple interaction potentials used in this work, and the equations governing the dynamics of our system. In Sec. III we present results for spatial distributions about particle positions, for the two particles, under varying settings, and discussion on time taken for the particle to reach the destination. Finally, Sec. IV, concludes the paper with some final remarks and proposed further work.

Model

Potential energy landscape

We study a system of two particles moving in two dimensions with coordinates $\boldsymbol{x}=(x,y)$, that experience a complex potential energy landscape [13], displayed in Figure 1. The dimensionless potential is defined as

where $\beta_{\mathrm{ref}}\equiv(k_{\mathrm{B}}T_{\mathrm{ref}})^{-1}$, $T_{\mathrm{ref}}$ is a reference temperature, $k_{\mathrm{B}}$ Boltzmann’s constant, and $\ell$ is a reference length-scale in our model. In the following, we use a dimensionless temperature $T^{*}=T/T_{\mathrm{ref}}$ and set the unit of length $\ell=1$. The first term in Eq. (1) is a stretched Mexican-hat type potential

The second term generates a small constant force parallel to the $x$-axis, and the third term

consists of four Gaussian contributions to generate local maxima and minima in the potential at specific points. The final term in $\phi_{\mathrm{G}}$ is added to make the local maximum near the origin even higher and thus suppress the probability of the system passing near this point $C$, and instead more often taking either of the two paths near the saddle points $B$ and $D$; see Figure 1.

The potential energy $\phi(x,y)$ exhibits five stationary points with coordinates

The point $A$ is a local minimum, while $E$ is the global minimum. The point $C$ is a local maximum, and the points $B$ and $D$ are saddle points.

Given these features of the potential in Eq. (1), a particle moving under the influence of $\phi(x,y)$ and starting in the potential well around point $A$ will typically follow one of two pathways from $A$ to $E$. The first pathway involves passing in the vicinity of the saddle point $D$, whereas the second pathway goes through the vicinity of point $B$, the other saddle point. Passing near $B$ requires the system to travel a longer distance, resulting in a longer pathway time.

A particle starting from point $A$ and going through the saddle at point $D$ has to overcome a barrier height $\beta_{\mathrm{ref}}\Delta\phi_{AD}\equiv\beta_{\mathrm{ref}}(\phi(x_D,y_D)-\phi(x_A,y_A))=2.7$. In contrast, the barrier height for going up to the saddle at point $B$ from initial point $A$ is $\beta_{\mathrm{ref}}\Delta\phi_{AB}\equiv\beta_{\mathrm{ref}}(\phi(x_B,y_B)-\phi(x_A,y_A))=1.3$. Thus, the system has two choices: either (i) go over the higher saddle at point $D$, with a relatively short distance to travel in state space, or (ii) go over the lower saddle point $B$, which however involves a longer journey.

Interparticle potential

To investigate the role of interactions between pairs of Brownian particles in the barrier crossing dynamics, we consider two paradigmatic interparticle potentials. The first is the Gaussian potential

where the parameter $\lambda_{\mathrm{G}}$ determines the range, $\varepsilon_{\mathrm{G}}$ determines the interaction strength, and $r \equiv|\mathbf{x}_{1}-\mathbf{x}_{2}|$, where $x_{1}$ and $x_{2}$ are the positions of the two Brownian particles. When $\varepsilon_{\mathrm{G}}>0$, this potential is often referred to as the Gaussian core model and it provides a good approximation for the effective interaction between the centres of mass of polymers in solution, under good solvent conditions [16]. In this case, $\lambda_{\mathrm{G}}$ corresponds to the polymer’s radius of gyration, $\varepsilon_{\mathrm{G}}$ is an energy scale on the order of $k_{\mathrm{B}}T$ [16] representing the strength of the interparticle interaction. Here we consider both the case $\varepsilon_{\mathrm{G}}>0$, corresponding to repulsive interactions, and also the case $\varepsilon_{\mathrm{G}}<0$, corresponding to an attractive interaction.

We also investigate the Lennard-Jones (LJ) potential

where $\varepsilon_{\mathrm{LJ}}$ is the Lennard-Jones energy scale (depth of the potential well), and $\sigma_{\mathrm{LJ}}$ the particle diameter. The potential minimum is at $r_{\mathrm{eq}}=2^{1/6}\sigma_{\mathrm{LJ}}$.

Stochastic dynamics

We study a system of two Brownian particles governed by the following overdamped Langevin equations

where $k=\{1,2\}$, $\gamma$ is the friction coefficient and the random force $\boldsymbol{\xi}_{k}$ represent independent Gaussian white noise. Henceforth, $\gamma$ is set equal to one, i.e. is absorbed into the timescale. Thus, the natural Brownian timescale for the system $\tau_{\mathrm{B}}=\beta_{\mathrm{ref}}\gamma\ell^{2}$, is also set equal to one. The random components of $\boldsymbol{\xi}(t)$ (we dropped the subscript $k$ here deliberately) obey

where $\delta(t-t')$ is a Dirac delta distribution, $\delta_{ij}$ is the Kronecker delta, with $i,j=\{x,y\}$, and $\langle\cdot\rangle$ denotes the ensemble average. The second term on the right-hand side of Eq. (7) is the force stemming from the interparticle potential. In the following, this term will be set equal to either $V_{\mathrm{G}}(r)$ or $V_{\mathrm{LJ}}(r)$ – see below Secs. III A and III B, respectively.

Equations (7) and (8) are based on the assumption of a Markovian process, which means that the system has no memory. However, this is not always a trivial matter. The interested reader can consult the Markovianity test along a reaction coordinate proposed in [34].

Simulation initialization and methods

We initialize the system with both particles at or in proximity of point $A$ (a local minimum of $\phi(x,y)$) and then we numerically integrate Eqs. (7) and (8) using the standard Euler–Maruyama finite difference method [35], with a time step $\Delta t = 10^{-3}$, for a total of $2 \times10^{7}$ time steps. to increase statistical accuracy we performed $5 \times10^{5}$ independent runs. Spatial derivatives are computed on a square grid, with spacing $\Delta x = \Delta y = 0.01$. In order to compute the first passage time to reach the global minimum at $E$, we define ‘reaching $E$’ as the first time when a particle $k$ satisfies the condition $\lvert x_k - x_E\rvert< 0.1$.

For the case when the particles interact via the Gaussian potential $V_G(r)$, we set the initial position of both particles equal to point $A \equiv(x_A,y_A)$. The Gaussian potential in Eq. (5) does not diverge as the particle separation $r \to0$, unlike typical interatomic or colloidal potentials, allowing configurations in which the centers of the two particles overlap. When this potential is used to model polymers in solution, this represents the centers of mass of two polymer chains overlapping, reflecting the non-compact nature of polymers. Also, even with overlapping centers of mass, chains can fluctuate without intersecting due to their soft, interpenetrable structure [16]. The sets of parameters for the Gaussian potential (5) explored in this work are $\epsilon_G = \{-10,-5,0,5\}$ and $\lambda_G = \{0.2,0.5,0.8\}$. We thus explore both attractive and repulsive interparticle potentials.

The LJ potential diverges $V_{LJ}(r) \to\infty$ as $r \to0$. Therefore, to avoid the singularity that arises from the overlapping of particles, we initialize the particles with a separation of $r_{eq}$ (potential minimum distance) around point $(x_A,y_A)$. That is, we assign the initial positions as $(x_A + 0.5r_{eq},y_A)$ and $(x_A - 0.5r_{eq},y_A)$ and label these two points as $P_1$ and $P_2$, respectively. Likewise, in this LJ case we refer to the particle starting at $P_1$ as ‘particle 1’ and the other as ‘particle 2’. During the simulation, to avoid excessively large forces that arise at small interaction distances $r$, causing numerical instability, we replace the interaction force for $r < r_c$, where $r_c = 0.9\sigma_{LJ}$, with the value of the repulsive force at $r = r_c$. The value of $r_c$ is chosen to strike the right balance between avoiding singularity and preserving the dynamics of the pair potential and is consistent with best practices in molecular dynamics simulations, where similar short-range cut-offs are used to ensure numerical stability [36]. The sets of parameters for the LJ potential (6) that we consider are: $\epsilon_{LJ} = \{0,3.33,6.67,10\}$ and $\sigma_{LJ} = \{0.33,0.67,1.0\}$.

To study the features of spatial and temporal probability distributions of the two interacting particles, we considered both the Jensen–Shannon divergence (JSD) [37] and order statistics for reaching the minimum at $E$. The JSD $J(p_1(x)\Vert p_2(x))$ is a symmetrized and smoothed version of the Kullback–Leibler divergence $D(p_1(x)\Vert p_2(x))$ [38] and is defined as

where $M(x) = \frac{1}{2}(p_1(x) + p_2(x))$ is a mixture distribution of $p_1(x)$ and $p_2(x)$ and

A value of JSD close to zero indicates a high similarity in the two distributions, while the larger the value, the more different the distributions.

We investigate the distribution of the arrival time $t$ for either of the two particles to be the first to reach the destination. This problem is closely related to the first-passage time problem [39]. The minimum of individual arrival times is an ordered statistics [40]

If the arrival times of each of the particles are independent events (which is the case when the interaction parameters $\epsilon_G = 0$ or $\epsilon_{LJ} = 0$, for Gaussian or LJ interactions, respectively), then the probability density function (PDF) of $t_{(1)}$ reads

where $p_1(t)$ is the PDF of the arrival time for a single particle, and $G_1(t)$ is the cumulative distribution function (CDF) of $p_1(t)$. As a result, the arrival time distribution for pairs of noninteracting particles is shifted toward earlier times, with a higher concentration of probability mass in the lower $t$ range. In other words, $p_{t_{(1)}}(t)$, which is henceforth denoted as $p(t)$, is expected to be narrower than the arrival time distribution for a single particle, assuming the particles are noninteracting.

For example, if the arrival time for a single particle starting at point $A$ follows an exponential distribution, with decay parameter $\tau$, defined by

for all $t>0$, $\tau>0$, this then gives $G_1(t)=1-e^{-\frac{t}{\tau}}$, and so Eq. (12) gives

Note that in Ref. [13] it was demonstrated that, except at very early times, the expression in Eq. (13) provides a good approximation for the single-particle arrival time distribution. Equation (14) therefore implies that the probability density for the arrival time of the first of a pair of particles to arrive at the destination $E$ should exhibit a faster decay over time than the single-particle case. As we show below, for interacting particles, this is not always the case.

Results and discussion

This section presents our results for the dynamics defined in the previous section. We discuss first the Gaussian pair interaction case, followed by the LJ potential case.

Gaussian pair potential

Location

Figure 2 shows the two-dimensional distribution of the positions of one particle when the other has reached the destination. Operationally, reaching the destination point $E$ is defined as the time when the particle $k$’s distance from $E$ satisfies $\lvert x_k-x_E\rvert<0.1$. We observe that the position distribution of the particle not at point $E$ is highly sensitive to changes in the values of the pair interaction parameters, namely $\varepsilon_G$, the interaction strength, and $\lambda_G$, the interaction range, i.e. the ‘size’ of the particles. The values of these define the dynamics of particles interacting through the Gaussian potential (5).

Vector (PDF)

In the strongly attractive regime, e.g. the $\varepsilon_G=-10$ cases in Figure 2, we see that the second particle is predominantly localized in the region close to the destination point $E$, at the time when the first particle reaches the destination. This is because the strong attractive interaction has effectively ‘glued’ the particles together, so that they predominantly move as a bound pair. The density distribution of the second particle exhibits two prominent local maxima, one either side of point $E$, and a much smaller local maximum at point $A$, the start point, corresponding to the rare cases when one of the particles gets left behind at the start point by the other particle.

For the smallest value of the interaction range that we explored $\lambda_G=0.2$, the second particle’s distribution is very tightly localized around the destination point $E$ when $\varepsilon_G<0$, sticking close to the first particle due to the strong attraction. Upon increasing $\lambda_G$, we observe that the distribution becomes broader around the destination point $E$. This can be attributed to the fact that although the pair of particles are tightly bound to one another, a larger interaction range (i.e. particle size) means that the distance between the pair of particles can be much greater.

In moderately attractive cases, such as the $\varepsilon_G=-5$ case in Figure 2, the second particle can detach from the first particle and remain more often in the vicinity of the start position (the local minimum $A$), while the first particle crosses over one of the barriers to point $E$. As a result, the distributions in Figure 2 exhibits several prominent local maxima: one close to the start point $A$ and then a pair either side of the destination point $E$. We can infer from this that the weaker attractive force allows the second particle to maintain greater independence, and its motion is less impacted by the dynamics of first particle.

In the absence of interaction, $\varepsilon_G=0$, Figure 2 shows that the second particle behaves completely independently of the first and exhibits a very broad distribution, but mainly remains confined to the vicinity of its initial position $A$, regardless of the fact that the first particle has moved to the destination at the point $E$. This behavior is expected because the lack of interaction. Note that further details of the time evolution in the single particle case can be found in Refs. [12, 13].

When the Gaussian interaction between the particles is repulsive, such as the $\varepsilon_G = 5$ case, displayed in bottom row of Figure 2, we see that the distribution of the second particle is largely independent of the position of the first, with a broad distribution similar to the non-interacting case, except for in the case of large particle sizes $\lambda_G$. In this repulsive and large interaction range scenario, we see an increased density around the saddle point $B$, defining the barrier with a lower height at a longer distance. This suggests that the preferred path of the second particle to reach the destination is through the saddle point $B$. We note that the time to reach the destination for large particles is very small [see Figure 3(a)]. This finding suggests that repulsive forces push the particles apart such that the first particle to reach the destination at point $E$ effectively prevents the second particle from approaching the destination, particularly when the particle size is large.

These observations illustrate that when one particle reaches the destination first, the spatial distribution of the second particle changes in a systematic way from its initial position $A$ to the destination point $E$. As a result, strong repulsive forces drive separation, while attractive forces encourage particle cohesion, with the effect becoming more pronounced with an increase in particle size.

In order to check the stability of the numerical method used for the time integration of Eq. (7), we also produced some results with the smaller time step of $\Delta t = 10^{-4}$, which is an order of magnitude smaller than the value used to produce Figure 2. This smaller value of $\Delta t$ results in spatial distributions that are visually identical to those displayed in Figure 2.

Arrival time

We now discuss the arrival time distributions, i.e. the time taken by the first particle to reach the destination, for the case where the particles interact via a Gaussian pair potential. Figure 3 shows the probability density of arrival times $p(t)$ for different values of $\varepsilon_G$ and $\lambda_G$. Note that we do not display results for the noninteracting case where $\varepsilon_G = 0$, because these (effectively single-particle) results are already presented in Ref. [13]. It is immediately apparent from comparing the ranges of the three different $t$-axes that the typical time for the first particle to travel from $A$ to $E$ increases as we move from repulsive $\varepsilon_G > 0$ to attractive $\varepsilon_G < 0$. This is due to the particles sticking together when the attraction is strong enough. The resulting pair take a longer time to reach the destination than a single particle on its own.

The insets of Figure 3 show the temporal dependence of $\ln(p(t))$. Except at very short times, the distributions $p(t)$ have the following exponential form

where the constant $\chi$ (obtained by fitting to the data) is close, but not exactly equal, in value to $1/\tau$. This is because of the deviations at short times from the exponential distribution. The exponential form in Eq. (15) would be exact if the rates $k_{AB}$ and $k_{AD}$ for crossing each of the barriers at $B$ and $D$, respectively, were constant in time. However, as discussed in Refs. [12, 13], at early times these take different values from their long-time limiting values, so there are deviations from the exponential form in Eq. (15) at early times $t$. In Figure 3 we also display exponential fits to the tails of the data. We observe that for each $\varepsilon_G$ value, upon increasing $\lambda_G$, the reciprocal of the characteristic time $\tau^{-1}$ increases (see inset legends).

Lennard-Jones potential

We now discuss results for pairs of Brownian particles interacting via the LJ pair potential in Eq. (6). Unlike for the Gaussian potential, we cannot initialize both particles fully overlapping at the same start position $A$, due to the core repulsion of the LJ potential. As discussed in Sec. II D, we initialize particle 1 at point $P_1$ and particle 2 at point $P_2$, with coordinates $(x_A+0.5r_{eq},y_A)$ and $(x_A-0.5r_{eq},y_A)$, respectively. Since the LJ equilibrium distance $r_{eq}=2^{1/6}\sigma_{LJ}$ depends on the interaction potential parameter $\sigma_{LJ}$, this raises the natural question of how the dynamics of the system is affected by the choice of initial conditions. In the following, we vary the value of $\sigma_{LJ}$, with the largest value considered being $\sigma_{LJ}=1$. This sets the maximum distance between the particle positions either side of the initial point $A$ to be $r_{eq}=2^{1/6}\sigma_{LJ}\approx1.12$. Thus, we ensure that the size of the particles is not large enough to cover the entire attractive basin around point $A$ in the external potential.

To investigate the impact of the initial particle separation $r_{eq}$ on the probability of particle 1 (starting at $P_1$) to reach $E$ first, we calculated this probability for the range of values $\sigma_{LJ}=\{0,0.1,0.2\cdots,1\}$ and for $\varepsilon_{LJ}=\{0,0.5,2,5,10\}$. Our results are displayed in Figure 4. When both $\sigma_{LJ}=0$ and $\varepsilon_{LJ}=0$, each particle has the same probability (50%) of reaching the destination; this case corresponds to two free point-like particles. Upon increasing $\sigma_{LJ}$, the probability of particle 1 (starting at $P_1$, located to the right of the point $A$) to reach the destination increases. This is because an increase in the value of $\sigma_{LJ}$ makes the start point $P_1$ slightly closer to the saddle point $D$ that defines the higher energy barrier, and much farther away from the saddle point $B$ that defines the lower energy barrier. This makes it easier to climb the high energy barrier to reach the destination $E$. Strikingly, the probability increases further with higher values of $\varepsilon_{LJ}$, but becoming non-monotonic as a function of $\sigma_{LJ}$. We believe this is due to the fact that for higher values of $\varepsilon_{LJ}$, the pair strongly bind, forming a dimer of length $\sigma_{LJ}$, that is relatively inflexible. For some dimer lengths, the resulting ‘rod’ is commensurate with the features (minima and valleys) of the external potential, while for other lengths it does not fit well, strongly affecting which way the dimer travels through the potential energy landscape away from the start point.

Vector (PDF)

Because the initial start points $P_1$ and $P_2$ have an impact on the probability of the respective particles making their way to point $E$, we must discuss the location and arrival time statistics for each particle separately. Thus, in the next section, we provide a discussion of the location statistics of particles 1 and 2 at the moment if/when the other particle has reached point $E$ first, followed by a discussion of their respective arrival times, for various values of $\varepsilon_{LJ}$ and $\sigma_{LJ}$.

Location

Figure 5 shows the probability distributions for the position of one particle, when the other has reached the destination point $E$. Columns 2, 4, and 6 of Figure 5 refer to the position of particle 1 (starting at $P_1$), when particle 2 has reached $E$ for four values of $\varepsilon_{LJ}$ (rows) and three values of $\sigma_{LJ}$ (columns). Columns 1, 3, and 5 correspond to the position distributions of particle 2, when particle 1 has reached $E$ first, for the corresponding values of $\varepsilon_{LJ}$ and $\sigma_{LJ}$.

Vector (PDF)

In the $\varepsilon_{LJ}=0$ (noninteracting particles) case, the distribution of particle 1, when particle 2 has arrived $E$, and the distribution of particle 2, when particle 1 has reached the point $E$, are very similar for the smaller values of $\sigma_{LJ}$, only becoming noticeably different for the largest value ($\sigma_{LJ}=1$). This is due to the fact that the start points of the two particles are different, so the distributions are different, even though they are noninteracting. In the absence of interaction, the distributions of each particle spreads throughout the landscape of $\phi(x)$ starting from the local minima at point $A$, making their way to the saddle point $B$; see Refs. [12, 13] for further discussion of this effectively single independent dynamics case. When the LJ diameter is small ($\sigma_{LJ}=0.33$), the distributions for the three cases $\varepsilon_{LJ}=\{3.33,6.67,10\}$ do not drastically differ from the non-interacting case $\varepsilon_{LJ}=0$. This is on account of the particle’s diameter being rather small so that the interaction quickly decays to zero. Furthermore, the concentration of the particle that did not reach the destination is mostly concentrated around the local minimum point $A$ and above the saddle point $D$.

For $\sigma_{LJ}=0.67$, we observe that the presence of a non-zero interaction causes a drastic drop in particle concentration around the saddle point $D$ and the appearance of local maximum in probability density in the potential valley between saddle point $B$ and the global minimum point $E$. An increase in the value of $\varepsilon_{LJ}$ makes this structural change in the probability density more dramatic.

For $\sigma_{LJ}=1$ (and nonzero $\varepsilon_{LJ}$), we observe that the probability density concentrates mostly around the local minimum $A$ in a crescent shape skewed towards the saddle point $B$.

Figure 6 shows the JSD metric (9) for the difference between the probability distributions shown in Figure 5 of the positions of the particles. It is clear that for $\sigma_{LJ}=0.33$ the probability distributions of the two particles are rather similar, reflected by the fact that the value of JSD is relatively small. In contrast, for a moderate value of $\sigma_{LJ}=0.67$, we observe that upon increasing $\varepsilon_{LJ}$ the JSD grows, reflecting the fact that the probability distributions grow increasingly different. However, for the largest value of $\sigma_{LJ}$ explored here, we observe a maximum in JSD at $\varepsilon_{LJ}=1.67$ followed by a slowly decreasing trend, suggesting a maximum dissimilarity in the probability distribution of two particles as $\varepsilon_{LJ}$ and $\sigma_{LJ}$ are varied. These observations made through the lens of the JSD in Figure 6, supplement our finding from Figure 5. A particularly striking aspect of Figure 6 is that the JSD is a non-monotonic function of $\sigma_{LJ}$ for $\varepsilon_{LJ}>7$, highlighting again the non-trivial dynamics of dimers in the potential $\phi(x)$, as a function of dimer length $r_{\mathrm{eq}}$.

Vector (PDF)

Arrival time

Next, we examine the temporal evolution of our two-particle system, illustrating how LJ interactions induce changes in the statistics of the arrival times. Figure 7 shows the arrival time distribution $p_1(t)$ of particle 1, when particle 2 has reached $E$ first, whereas Figure 8 shows the arrival time distribution $p_2(t)$ of particle 2, when particle 1 has reached $E$ first.

From Figure 7, we observe that the arrival time distribution $p_1(t)$ is strongly peaked around small times for all $\varepsilon_{LJ}$. However, this peak increases in magnitude upon increasing $\sigma_{LJ}$. The strong peak at relatively small times is to be expected, since particle 1 starts slightly closer to the destination $E$, allowing it to reach it earlier.

The insets of Figure 7 show semilogarithmic plots of $p_1(t)$ for the same data in the primary plots. Except at very short times, the distributions are well fitted by straight lines, suggesting that the tails of $p_1(t)$ can be modeled by exponential distributions, similar to Eq. (15)

where the constant $\chi_1$ is not exactly equal to $1/\tau_1$ but is $\approx1/\tau_1$, due to deviations at short times. These exponential fits are shown as straight lines in the Figure 7 insets.

Given the behavior of the spatial distributions in Figure 5, one might expect $p_1(t)$ to exhibit a strong dependence on $\sigma_{LJ}$, but infact this is not the case. An increase in $\sigma_{LJ}$ for any given $\varepsilon_{LJ}$ does give an increase in $\tau_1^{-1}$, as shown in the insets, but the variations between the various panels (corresponding to changing the value of $\varepsilon_{LJ}$) leads only to small variations.

Figure 8 shows that the arrival time distribution $p_2(t)$ is broader than the distribution of $p_1(t)$ for all $\varepsilon_{LJ}$ explored here. Furthermore, also in this case there is a dependence of the $p_2(t)$ distribution on $\sigma_{LJ}$. In the insets of Figure 8, we show the semilogarithmic plot of $p_2(t)$, for the same data set utilized in the primary plots. These plots again appear to be straight lines, except at very short times, of similar nature of exponential distribution defined in Eq. (15) and hence the updated distribution $p_2(t)$ for particle 2 reads

where the constant $\chi_2$, $\tau_2$ are defined in the same way.

Comparing the arrival time distributions $p_1(t)$ and $p_2(t)$, for particles 1 and 2 respectively, displayed in Figure 7 and Figure 8, indicates that particle 1 reaches the destination more quickly than particle 2, which is to be expected, given the starting position $P_1$ is closer to the destination point $E$.

Figure 9 shows the dependence of the time relaxation coefficients $\tau_1$ and $\tau_2$ defined in Eqs. (16) and (17) on the interaction parameter $\varepsilon_{LJ}$ of the LJ-pair potential. When $\sigma_{LJ}$ is small compared to the size of the wells in the external potential $\phi(x)$ and the distance between the wells, we see from Figure 9 that particle interactions increase the overall time for particles to reach the destination $E$. However, there is a nonmonotonic effect, and they are the slowest when

Vector (PDF)

$\varepsilon_{LJ} \approx6$.

When $\sigma_{LJ}=1$, the values of both $\tau$ remain largely unaffected by variations in $\varepsilon_{LJ}$, that is, as $\sigma_{LJ} \to\lvert x_C-x_B\rvert$, the values of $\tau_k$ are relatively insensitive to changes in $\varepsilon_{LJ}$, potentially implying underlying stability or saturation in the system response. In fact, when $\sigma_{LJ}$ is large, particle interactions speed up the particles to reach the destination.

Concluding remarks

We have studied the (noise) activated barriers-crossing process, and specifically the influence of particle interactions, of a pair of particles crossing barriers in systems with competing pathways. The external potential $\phi(x)$ considered is of interest because standard reaction rate theory (RRT) predicts that the most probable route is to cross the lowest-energy barrier [12, 13]. Thus, the expectation from RRT is that this longer route should dominate the dynamics of the system. However, even when there is a single particle, the standard RRT framework can fail to accurately capture the dynamics, becoming critically dependent on the timescale over which the system is observed or sampled [12, 13]. Such barrier-crossing problems involving competing pathways have broad implications.

Here, we have considered the case when there are a pair of particles available to cross the barriers and focussed on the influence of particle interactions on the (rare event) barrier-crossing processes. We have explored how interactions affect the choice of route taken and the time taken. Two well-known interaction potentials, the Gaussian pair potential and the LJ potential, were considered. In the LJ case, the interacting particles cannot overlap, so cannot start at exactly the same point, as was considered in the Gaussian pair potential case. Taken all together, our results show that particle interactions can play a crucial role in determining which particle crosses one of the barriers first and how long it takes. When there are strong attractive interactions, the pair of particles bind together and can sometimes be treated as a single larger particle. However, the presence of a repulsive core to the interaction, which introduces an additional length-scale, can have a major additional influence. When such particles bind together, the length of the resulting dimer and how that fits (or not) within the potential wells of the confining potential has a strong impact on the resulting first-passage time statistics. When interactions are purely repulsive, we generally find barrier-crossing occurs faster, since one particle helps to ‘push’ the other over the barrier.

The main findings of our results on barrier-crossing of pairs of particles are that: (i) the initial conditions of the particles in the external potential matter. (ii) Particle interactions introduce additional physics, which can greatly change the barrier crossing rates and choice of routes. (iii) Particle interactions not only introduce additional (binding) energy considerations, they also introduce additional length-scales, resulting in a complex interplay with the geometry of the confining potential. We expect these observations to have important implications for our understanding of chemical and biological reactions, for example stochastic target-search and confined-transport processes in which several particles simultaneously explore a structured energy landscape [41, 42]. Competing pathways and transition times are important in molecular recognition and drug–target binding, where molecular crowding can modify recognition kinetics [43].

Acknowledgments

GA gratefully acknowledges Loughborough University Doctoral College UK and HEC Pakistan for funding support under “Overseas Scholarships for PhD in Selected Fields Phase-III Batch 3”.

References