Particle dynamics and transport enhancement in a confined channel with position-dependent diffusivity

This work focuses on the dynamics of particles in a confined geometry with position-dependent diffusivity, where the confinement is modelled by a periodic channel consisting of unit cells connected by narrow passage ways. We consider three functional forms for the diffusivity, corresponding to the scenarios of a constant (D0), as well as a low (Dm) and a high (Dd) mobility diffusion in cell centre of the longitudinally symmetric cells. Due to the interaction among the diffusivity, channel shape and external force, the system exhibits complex and interesting phenomena. By calculating the probability density function, mean velocity and mean first exit time with the Itô calculus form, we find that in the absence of external forces the diffusivity Dd will redistribute particles near the channel wall, while the diffusivity Dm will trap them near the cell centre. The superposition of external forces will break their static distributions. Besides, our results demonstrate that for the diffusivity Dd, a high dependence on the x coordinate (parallel with the central channel line) will improve the mean velocity of the particles. In contrast, for the diffusivity Dm, a weak dependence on the x coordinate will dramatically accelerate the moving speed. In addition, it shows that a large external force can weaken the influences of different diffusivities; inversely, for a small external force, the types of diffusivity affect significantly the particle dynamics. In practice, one can apply these results to achieve a prominent enhancement of the particle transport in two- or three-dimensional channels by modulating the local tracer diffusivity via an engineered gel of varying porosity or by adding a cold tube to cool down the diffusivity along the central line, which may be a relevant effect in engineering applications. Effects of different stochastic calculi in the evaluation of the underlying multiplicative stochastic equation for different physical scenarios are discussed.


Introduction
The dynamics of many practical systems are restricted within the confines of structured and inhomogeneous environments. Examples include the motion of particles in specific Lorentz gases, e.g., a triangular lattice of immobile reflecting disks [1,2], protein folding [3][4][5] and polymer looping [6,7] in was investigated. Moreover, diffusion coefficients with a power-law form [44,67], or modelled by a Feller process [68] were studied. Cherstvy used a diffusion coefficient depending on the radial distance to a cell centre, corresponding to a radially dependent mobility [36] mimicking the local diffusivity variation in biological cells [11]. In fact, being opposite to a high-mobility centre, the diffusion with a low-mobility centre is also meaningful and applicable in practical systems. In addition, instead of a dependence on the purely geometric distance r = x 2 + y 2 , it is also reasonable to add a modulation parameter to adjust the dependence on the two coordinates, i.e. r = σx 2 + y 2 . In this paper, we consider two kinds of position-dependent diffusivities with distance modulation parameter, corresponding to high-mobility and low-mobility unit cell centre (or a heat source and a cold flow) respectively. The dynamics of particles in a confined channel under such diffusivities and the longitudinal transport behaviours are analysed in detail.
This paper consists of the following sections. In section 2, the model of channel diffusion with varying diffusivity is introduced. Two forms of position-dependent diffusion coefficients as well as their major properties are described. We also outline the numerical simulation method. In section 3, the probability density functions for different diffusivities are calculated. The influences of diffusivity, channel shape and external forces are presented and explained. In section 4, the transport properties of particles and influences of the position-dependence are studied by calculating the mean velocity. In section 5, the mean first exit time and first exit time distribution from a cell centre to either bottleneck with vicinal unit cells is determined to analyse the escaping behaviours. Finally, conclusions are drawn in section 6.

Model description
The model we study here consists of a confining channel in which the test particle experiences a position-dependent diffusion coefficient. The channel has a deformable structure, whose upper bound is defined as Correspondingly, the lower bound is −w (x). Here L is the period of w (x), L 1 = 1 + q L/2 and L 2 = 1 − q L/2. By adjusting the parameter q, one obtains different shapes of channels, such as left-skewed (q < 0), symmetric (q = 0) and right-skewed (q > 0) structures, as shown in figure 1. The maximum width of the channel is 2(b + 2h), and the width of the bottleneck to vicinal chambers is 2b. When q = 0, the widest position is located in the middle of the unit cell. This leads to a symmetric form of w (x). For any nonzero q, the widest position is always biased either to the left or right with respect to the cell mid-point. This leads to an asymmetric w (x). In the following, we denote the part labelled L 1 as the left-hand side (lhs) of a unit cell. Similarly, the part labelled L 2 is the right-hand side (rhs) of a unit cell, see figure 1(a). Generally, the movement of a particle confined by the channel ±w(x) can be modelled by the second-order Langevin equation where m is the mass of the test particle, r = (x, y) is the coordinate, γ is the friction coefficient, F (r) = [F, 0] T is the potential force, D(r) is the position-dependent diffusion coefficient, and η (t) = η x (t) , η y (t) T is a vector of independent Gaussian white noise components with zero mean and correlation function η i (t) η j (s) = δ ij δ (t − s) , i, j = x, y. In this work, we assume an environment with a constant friction coefficient γ determined by the ambient liquid. F 0 is a constant external force, under which particles within the channel will perform a net flow for nonzero external force, such that the mean particle velocity is non-zero. D(r) is a position-dependent diffusion coefficient. As mentioned above, we envisage this spatial modulation to stem from the local porosity of the medium determining the local mobility of the particle, or, alternatively, from a temperature field in the sense of weak gradients such that the thermophoretic effects can be neglected. When we consider an overdamped particles, i.e. γ r, the inertial term can be neglected, namely, the system (2) can be reduced to the overdamped forṁ x = F + 2D x, y η x (t) , y = 2D x, y η y (t) .
(3) Equation (3) can be also obtained by assuming m → 0, which is treated by the Itô interpretation [51,52] (see above). Thus, we evaluate equation (3) within the Itô integral in this work, unless noted otherwie. However, we also provide results with isothermal integral interpretation in the appendix for comparison. In [44], Cherstvy et al employed a power-law form D (x) = D 0 |x| a of the diffusivity in a one-dimensional system, where a can be either positive or negative. Whereafter, in the two-dimensional case, they extended it to D x, y = D 0 A/ A + x 2 + y 2 , where A is a small positive number to avoid 'death' or 'explosion' of D x, y at the origin [36]. For such a D x, y , the diffusivity depends symmetrically on both coordinates. However, in complex environments [11], the diffusivity may be coordinate-sensitive, namely, the gradient of D x, y is asymmetric. In this work, due to the intrinsic asymmetry of the channel geometry we take the asymmetry into account. Without loss of generality, we discuss the following two forms of D x, y where D 0 is a constant, and A is set A = 0.01 here. r σ x, y describes the distances between particles and the cell centre, which is defined as r σ x, y = σx 2 + y 2 . c = 1/ max r σ x, y is the inverse of the maximum of r σ x, y , which is superimposed to scale the distance weights.x is determined bȳ where the operator · means rounding up to an integer. Namely, |x| is the distance to the cell centre in the x axis. σ is to modulate the dependence of r σ x, y on the two coordinates. When σ < 1, r σ x, y depends more on the y direction than the x direction, while for σ > 1, r σ x, y depends more on the x direction than the y direction. Namely, a small change on x will lead to a large difference of r σ x, y . In this case, r σ x, y is x-sensitive. When σ = 1, r σ x, y is the classical distance function, which symmetrically depends on both coordinates [36].
Generally speaking, D d x, y is inversely proportional to the distance between the current position and the cell centre (kL, 0) (k is an integer). With the increase of the distance to the centre, the diffusivity increases, see figure 1(b). At the centre's position x, y = (kL, 0), the diffusion coefficient D d x, y = D 0 is the largest. Particles located here are the most mobile ones. For a given position x, y = (kL, 0), D d x, y < D 0 , i.e., the diffusion coefficient here is smaller compared with the cell centre. For the case D x, y = D m x, y , the situation is opposite. D m x, y is proportional to the distance between x, y and the cell centre (kL, 0). With the increase of the distance to the centre, the diffusivity increases, see figure 1(b). At the cell centre, D m x, y takes its minimum value AD 0 , and the particles here have the lowest activity. At other positions x, y = (kL, 0), larger r σ x, y implies larger D m x, y . Correspondingly, the higher the mobility of the particle, and the more likely it is for particles to move away from their original positions.
Some related studies have shown that when the diffusion coefficient is a constant D 0 , in the case of no bias, particles at each position are equally mobile, and homogeneously fill the entire space of the channel [2,19]. When we consider a non-constant diffusion coefficient, it is clear that the system may exhibit new interesting phenomena, such as inhomogeneous distribution of the particles. In this work, we intend to investigate the transport of particles driven by an external force in a position-dependent diffusive environment in the Itô sense of the stochastic equation. For ease of argument, we neglect wall-induced diffusion effects and hydrodynamics. Moreover, a low particle concentration and a large constant friction coefficient are assumed to make sure that thermophoretic effects are very weak when the diffusivity results from the temperature gradient. However, we also provide results with respect to the isothermal integral which is relevant to thermophoresis [55,56] in the appendix. To explore and characterize the effects of the position-dependent D x, y , some statistical indicators are calculated, including the steady-state probability density function (PDF), the mean velocity, and first exit time density (FETD) and mean first exit time (MFET). Their definitions are given below.
To estimate these quantities, we apply a numerical method, in which stepwise iterations are performed to reflect the influences of the position-dependent diffusivity. Considering the superposition of confined channels, we assume a reflective boundary condition here. In the calculation, we add an alternative noise-induced term corresponding to the Itô and isothermal integral form respectively. The detailed algorithm applied here reads [15,26]: Step 1. According to the discrete form of equation (3), iterate the position r i = (x i , y i ) at the ith step, as where Δt is the time step, and ζ x i and ζ y i are independent Gaussian random numbers obeying N(0, 1). α is a parameter taking the value 0 or 1. When α = 0, equation (3) corresponds to the Itô interpretation; α = 1 to the isothermal interpretation.
Step 2. If r i is inside w(x), update i to i + 1 and update D x i−1 , y i−1 to D x i , y i , and then return to step 1 to calculate the (i + 1)th step's position r i+1 = (x i+1 , y i+1 ); Step 3. Else if r i is outside w(x), the particle will collide with w(x). Then, the dichotomy and elastic reflection are used to deal with this case as follow: (a) Find the intersection (denote as point O(x p , y p )) between w(x) and the line from r i−1 = (x i−1 , y i−1 ) to r i by using the dichotomy. (b) Calculate the slope k of w(x) in the point O(x p , y p ). (c) Based on the elastic reflection condition, calculate the position r i = (x i , y i ) of the particle after the collision, as (d) Determine the relationship between r i and w(x). If r i is inside w(x), let r i = r i , return to step 2; otherwise, considering O and r i as points having the same meaning as r i and r i−1 , and repeat step 3.

Steady state probability density function (PDF)
The PDF is a useful quantity to evaluate how particles spread in the channel. In an open environment, the movement of particles will not be limited by boundaries. There, the probability evolution can be described by a Fokker-Planck equation (FPE) with natural boundary conditions P |x| → ∞, |y| → ∞, t = 0, in the sense of Itô integral ∂P x, y, t ∂t where Δ is the Laplace operator. As known, when F = 0 and D x, y = D 0 , the solution is P x, y, t = 1 4πD 0 t exp − x 2 +y 2 4D 0 t according to Einstein's results. With increasing t, for a fixed point, P x, y, t decreases quickly to zero, which indicates that particles will spread widely on the plane. In this way, the probability that particles stay within a certain region will decrease. However, when we add a closed reflective boundary in such a region, particles will move inside without flowing out, and the total probability will always be unity. Taking the periodic boundary w(x) into account, one is mainly concerned with the dynamics in the x direction. In the y direction, due to the confined channel, it can be regarded to be homogenous, which means in the steady state ∂ 2 ∂y 2 P x, y = 0. Based on this assumption, the FJ equation is developed, starting from the FPE, to describe the PDF in confined channels, which under the condition of D x, y = D 0 takes the form [19,20,[23][24][25][26] where is the effective potential combining the influences of F and w(x). The steady-state PDF can be solved with a periodic boundary condition, i.e. P (X) = P (X + L). However, when D x, y = D 0 , the corresponding FJ equation has not been derived yet. Thus, in the following, we evaluate the PDFs with the numerical technique given above. The PDFs are calculated by in which I {·} is an indicator function and N is a large positive integer to ensure enough sample data. The overline has the same meaning as previously. Δx and Δy are the statistical intervals in the two directions, where x i and y j are the middle values of the ith and jth statistical interval, respectively. Next, we show our results and discussions with respect to the parameters q, F and σ. Figure 2 is the result for a constant diffusion coefficient, i.e. D x, y = D 0 . Under this condition, the dynamics of the particles induced by D 0 is essentially the same at any position. When F = 0, regardless of q, D 0 is the only factor affecting the movement of the particles, and causes a homogeneous distribution in the channel, as shown in figure 2(a). This phenomenon can be observed through the uniform colour in figure 2(a). To be precise, we here calculate the density of the marginal PDF, which is defined as the ratio between P(x) and the corresponding channel width. The horizontal dashed lines in figure 2(a) demonstrate the uniform distribution of particles in the y direction. This result is different from systems in open environments, but coincides with the assumptions and results of the FJ equation. When F = 0, both F and D 0 affect the transport of the particles, and lead to a completely different PDF as presented in figure 2(b). It shows that an F > 0 will destroy the original homogeneous distribution and pushes particles to the right. However, because of the bottleneck effect, particles are hampered to passage from one cell to another, and accumulate ahead of the right-hand side bottleneck. This increases the PDF of the right-hand side cell larger than that of the left-hand side part. We note that in practical simulations, the FJ equation is also valid for small F. This can be verified by the changes of colours in the y direction and the red dash lines. Although the red lines are not horizontal and straight, the difference is not excessively large. In this case, the FJ equation is approximately applicable.
Another interesting phenomenon is that the PDF is also related to the shape of the channel. When q < 0, the lhs of the cell is wider and the rhs is relatively narrower. When moving to the right, particles will likely collide with the walls closed to narrow passage ways. Thus, particles crowd ahead of the long narrow part, which results in the large red colour region in figure 2(b1). When q = 0, the rhs of the cell has a larger space than that of q < 0. The number of collisions is less than for q < 0. Hence, the red colour region becomes more concentrated and smaller. Conversely, for q > 0, the cell has an even larger space on the right side. This leads to a more concentrated PDF on the right.
The impact of D d (x, y) on the PDF is detailed in figure 3. The maximum of D d x, y locates at the cell centre (KL, 0) and gradually decreases with increasing r σ x, y . It is clear that σ affects the dependence and sensitivity of D d (x, y) along both axes: when σ > 1, D d (x, y) decreases faster in the x direction; when σ = 1, D d (x, y) decreases equally in both directions; when σ < 1, D d (x, y) decreases faster in the y direction. Thus, with increasing σ, the shape of D d (x, y) in the x direction becomes progressively narrower, while that in the y direction becomes wider and wider. As we discussed previously, particles tend to leave the region with large diffusivity and accumulate in the region with small diffusivity. Correspondingly, for F = 0 in figure 3(a), particles prefer to remain in the upper and lower corners for σ < 1. When we increase σ, their distributions gradually change from the vertical corners towards the bottlenecks. Moreover, with growing σ, the particles become more focussed at the bottlenecks. Similarly, for decreasing σ, they are more focussed at the vertical corners. These variations can be seen clearly by the marginal PDFs in figure 3(a4). With the increase of σ, the marginal PDFs change from unimodal to bimodal shapes. It is clear that when particles are focussed at the vertical corners, their passage from one cell to another will be impeded. However, even though particles gather at the bottleneck, it does not always accelerate the transport. This phenomenon will be demonstrated in detail later. When we calculate equation (3) with isothermal interpretation (α = 1), the results are significantly different. As shown in figure A1(a), although the PDFs still have the conditioned position preference, they are almost uniformly distributed with slight variations. This phenomenon is consistent with previous work at thermodynamic equilibrium. Related work involved fraction Brownian motion with the reflecting boundary condition can be found in [69,70]. However, it should be noted that the observed diffusion gradients come from particle-wall effects or the preciseboundary structure, which has a different regime with the temperature gradient. This is in contrast to our model, we assume a constant friction coefficient, where fluctuation-dissipation theorem does not holds. For non-equilibrium systems such as living biological cells only experiments can determine the correct calculus. In the limits considered here the Itô calculus appears appropriate.
When we add an external force, say F = 5 in our dimensionless units, the particles will be forced to move within the channel. Due to the narrow space at the bottlenecks, collisions with the walls will be extensive, and particles will be hindered to get through the bottleneck. Therefore, many particles will likely gather around there. Moreover, D d (x, y) has large values near the cell centre. This increases their tendency to move near the channel walls, while avoids crossing the cell centre. As shown in figure 3(a), when σ < 1, the particles gather at the vertical corners, and close to the channels. For F > 0, the interaction of the external forces and the gradient of diffusivity tends to push particles to collide with the walls, which results in large amounts of particles near the right-hand side bottleneck. This can be seen from the shift of the peak position of the marginal PDF, and the tendency will be enhanced by the sticky wall effects [71]. When σ = 1, the forced particles will have a higher probability to move away from the channels than for σ < 1. Therefore, it is easier for particles to pass through the bottleneck, and results in higher particles concentrations on the lhs of the cell. Correspondingly, the marginal PDF changes from a symmetric to a right-skewed shape. For even larger σ > 1 and F = 0, particles gather around the bottleneck, being separated by a high diffusivity region. Under the driving of F, one might expect that σ > 1 is a better condition than σ = 1 to promote the transport. In fact, when particles move from the left to the right, their trajectories will steer clear of the centre region (see the baby blue and white areas in figure 3). For σ < 1, the trajectories from the left to the right distribute around the cell centre with two clearly separating branches. For σ = 1, the trajectories concentrate near the cell centre with two very close branches, which is indeed more beneficial for transport than σ < 1. For σ > 1, although the branches merge into one, the trajectories spread widely, which is not certainly better than σ = 1. In next section, we show that the mean velocity for σ = 1 is indeed larger.
In figure 2(b), the rectification effect of the channel shape is shown. Here, we take q = 0.4 as an example to illustrate the PDF for D d (x, y). In this case, when the particles move to the right, the collisions with the channel walls may reverse the direction of mean motion, while for a flat wall, the collisions may only induce a difference in the moving direction without changing the sign of the velocity. Thus, in figure 3(c), we find that more particles gather along the rhs of the channels as compared to that of q = 0. Similar to figure 3(b), with increasing σ, more particles appear near the right bottleneck, and the PDF gets more condensed.
Comparing with the isothermal case, particles will not steer clear of the central area for F > 0. This is caused by the interaction between the noise-induced drift and the diffusion term, which will counteract the effect of D d (x, y), and leads to a relatively homogenous distribution for F = 0. Therefore, under the driving of an external force, particles are more concentrated in figures A1(b) and (c) than in 3(b) and (c).
In figure 4, we consider D m (x, y), which has a shape opposite to that of D d (x, y). According to the definition, the minimum of D m x, y locates at the cell centre, and increases with r σ x, y . For σ = 1, D m x, y has equal sensitivity to both coordinates. For σ < 1, D m x, y shows less sensitivity to the x coordinate, while for σ > 1, D m x, y depends more on the x coordinate. As the particles are more easily trapped at positions with low diffusivity, the PDFs have a peak at the centre with different shapes ( figure 4(a)). With increasing σ, the top views of the PDFs change from a horizontal ellipse-like shape (corresponding to σ < 1) to a vertical ellipse-like shape (corresponding to σ > 1). According to the marginal PDF in figure 4(a4), particles are more focussed around x = 0 for σ > 1. Decreasing σ will induce broader and broader tails of the marginal PDF. However, for the isothermal interpretation, although particles still prefer to stay in the central area, the PDFs vary slightly. They are close to an uniform distribution, due to the interaction between the noise-induced drift and the diffusion term, see figure A2(a).
In figure 4(b), we analyse the influences of F. For σ = 0.1 in figure 4(b1), F completely redistributes the particles and drives them to move from near the centre to the bottleneck. This causes the PDF to appear like a 'stick', spreading along the whole cell length L and getting restricted within a narrow band in the y direction. In this case, the particles mainly stay away from the walls, which will reduce the collision probability and benefit particles passing through the bottleneck. For σ = 1, F pushes the particles to the right, and shift the peak position of the PDF to near the right-hand side bottleneck. For σ = 10, we get a similar phenomenon as for σ = 1, but the PDF is less concentrated. From figure 4(b4), we find that a larger σ leads to a thinner left tail, whose peak position is also further away from the right-hand side bottleneck. This indicates that with increasing σ, particles gather near the right and wait to crowd out the bottleneck. Comparing these three scenarios, it shows that σ = 0.1 responds strongly to F, and such a parameter is convenient for particle transport. The PDFs for the isothermal case (see figure A2(b)) are generally similar to figure 4(b), because a large F = 5 overwhelms the influences of noise-induced drift. For a very small F, however, the PDF will be significantly different (not shown).
Basing on the influences of F, q will further reshape the PDF. When σ = 0.1, a broad right-hand side cell space and a horizontal stick-like distribution of particles make more particles move towards the right-hand side bottleneck. However, the steep right-hand side walls will prevent particles from passing through the bottleneck fluently, which induces particles to stay close to the walls. For large σ, the distributed area of the particles almost does not change, but the colour in this area becomes lighter, which is induced by a wider distribution in the y direction. For a comparison with the isothermal case, see figure A2(c).
In this section, the influences of position-independent D(x, y) on the PDF are studied with different parameters, including F, q and σ. In summary, all these parameters can modulate the shape of PDF. In the absence of F, particles tend to stay along the walls for D d (x, y), while near the cell centre for D m (x, y). A positive external force will drive the particles to the right direction further compared to the force-free case. Taking q into account, it will significantly reshape the PDF.

Transport properties
The steady-state PDF describes the static distribution of particles, and thus provides only very limited information about the dynamical evolution of the particles. To explore the transport behaviours of the particles in the presence of an external driving force F, we calculate the mean velocity v in this section. This is an important and commonly used indicator to measure the mean velocity of particles, and it is Here, x n (t) is the position of the nth realization at time t. N is the total number of trajectory. v is calculated here for the three kinds of D(x, y) with respect to the parameters F, σ and q. First, we consider the simplest case D x, y = D 0 . By solving the given FJ equation (7), one can derive the steady-state probability current J (x, t). As known, in a ratchet potential, the mean velocity can be obtained by J = v L. Thus, we get an explicit expression for D x, y = D 0 , which has proved to be accurate for small external forces in [23] However, to be consistent with the other two cases, here we use numerical simulations to obtain v to avoid inaccuracies, especially for large external driving forces. As shown in figure 5, for a fixed q, v increases monotonously with F as expected. This is intuitive. In addition, the effects of q on v are presented. In the absence of F, v is practically zero as it should. In fact, when q = 0, for F = 0 the particles move to either direction with equal probability, which induces a zero net current. When q = 0, however v is nonzero due to the ratchet effect, but with a very small magnitude not shown in detail in figure 5. Ratchet effects are always used to rectify the system current by setting an asymmetric periodic potential. When F > 0, particles are forced to move to the right. Besides, with increasing q, the channel changes from left-skewed to right-skewed shapes. Namely, the right-hand side space in a cell relatively increases. Thus, particles are increasingly directed to the right. However, a steep right-hand side channel boundary will increase the collisions between particles and walls, and make it difficult to locate the bottleneck. Consequently, a decreasing trend of v is generated for large q.
Secondly, v for D d x, y are studied in figure 6. It is obvious in figures 6(a) and (b) that the influence of F has the same trend as that of D x, y = D 0 . Namely, with increasing F, v increases monotonously, and for large F, it is almost linearly increasing. Moreover, for a fixed F, with increasing q, v decreases. In figure 6(b), the black line corresponds to the v for the case D 0 which is smaller than that of D d x, y . Generally, in the absence of external forces, a small diffusivity will induce a long time for the particles to move significantly away, and is thus detrimental for particle transport. However, under external forces, the particles are forced to move preferentially in one direction. A large diffusivity will result in extensive fluctuations and disturb the movement. Thus, a small diffusivity is preferred in this situation to get a larger v . Because D d x, y D 0 , it is no surprize that D d x, y induces a larger v . Considering σ, we find that v firstly increases for 0 σ 1,and then decreases very slowly with σ > 1. These slight changes are quite different from the significant differences of the PDFs in figure 3. Moreover, when σ = 0, v shows a sharp decline compared to σ > 1. For σ = 0, D d x, y only depends on the y variable, which leads to a ridge-like shape of the diffusivity along the x direction and induces a horizontal groove-like gap separating particles to the vertical corners (refer to figure 3). This will greatly increase the collision probability with the cell walls and induce a small v . Recalling figure 3(b), the trajectories from the left to the right have two separating branches for small σ; become two close-by paths for σ = 1; and spread widely for a large σ. Therefore, v reduces sharply for σ = 0, and decreases slowly for σ > 1.
In addition, we calculate v in the q − σ plane for F = 1 and F = 10, respectively. In figure 6(c), v generally increases with decreasing q. It is small for small σ and large q, while large for small q and σ ≈ 1. In fact, when q = 1, each cell has an upright wall at the right side. Such a channel shape will reflect the particles back repeatedly, and is detrimental for particles to locate the bottleneck. The black line corresponds to v of D 0 , which is smaller than v of D d x, y for small q, but larger for large q. Thus, even if D d x, y D 0 , for large q, its v can be smaller. In the inset, we calculate the ratio v σ=1 /v D 0 between v σ=1 and v of D 0 to examine the enhancement effect of D d x, y on v . It shows that for F = 1, the ratio increases slightly first and then decrease fast with q, which indicates that D d x, y enhances the transport for a small q and suppresses it for a large q. For F = 10 in figure 6(d), the variation of the amplitude of v induced by σ becomes relatively small. Furthermore, as shown in the inset, the ratio first increases and then decreases. This shows that for large q, the influences of σ on v change from suppression to enhancement with F = 10. Thus, a larger F can weaken the enhancement effect of D d x, y for small q, but increases it for extremely large q.
For the isothermal case, due to the interplay between F and the noise-induced drift, the result is relatively similar with the Itô case for a large F. However, for a small F, the variation of v becomes quite complex. As shown in figure A3(b), for F = 0.4, with increasing σ, v goes through three periods, i.e. decreasing, increasing and decreasing. However, σ = 1 is still an optimal point to achieve the largest speed for q = 0 and fixed F. Additionally, for a small F, v varies totally differently for q near −1 and 1, however, for a large F, v will always increase for different q. A detailed discussion is provided in the caption of figure A3.
Next, we consider the diffusion coefficient D m x, y in figure 7. From figure 7(a), we see a similar shape of v as in figures 5 and 6(a). Note that v for large F and large q is larger than those for D 0 and D d x, y . However, v for very small q does not change significantly. This indicates that channels with a steep left wall are not sensitive to the types of diffusivity. In figure 7(b), we find that v decreases monotonously with σ for a fixed F. This phenomenon is opposite to that for D d x, y which shows an increasing trend. As discussed before, D m x, y has its minimal value at the cell centre. This facilitates the trapping of particles and keeps them near the cell centre. Especially with increasing σ, the PDF of the force-free case changes from a transverse distribution to a vertical one. Thus, for small σ value, it is profitable for particles to pass through the bottleneck. In contrast, a large σ will impede the process. Consequently v decreases with increasing σ. It should be pointed out that when σ = 0, D m x, y only depends on the y variable, and has a groove-like shape along the x direction. Thus, particles can easily gather along the x axis and effectively avoid colliding with channel walls. This will significantly accelerate the directed movement of particles, which is an ideal condition for transport.
In figures 7(c) and (d), the influence of the parameters σ and q on v is analysed for F = 1 and F = 10, respectively. As presented in both figures, for a fixed q, v decreases monotonously with increasing σ. A small increase around σ ≈ 0 will induce an abrupt reduction of v , and the loss amplitude increases with q. In figure 7(c), for small σ, v keeps increasing with decreasing q; for large σ, v first decreases slightly and then increases. Additionally, from the black line we find that v values for large σ and large q are smaller than those for D x, y = D 0 . However, for small σ and small q, v is larger. Therefore, even if D m x, y D 0 , the channel shape can significantly rectify the current to reduce v significantly for large σ. In the inset, all ratios are larger than 2.5, which means D m x, y accelerates the particle transport more than D d x, y . Moreover, a moderately large q will cause a higher enhancement. We note that the case σ = 0 is quite special, as it does not vary perceivably with changing q. This is also true for F = 10 in figure 7(d). But there are also some differences between small and large F. In figure 7(d), v increases with decreasing q, and it is always larger than for D x, y = D 0 . The inset shows that the ratio increases monotonously. This illustrates that for small q a large F suppresses the enhancement effect of D m x, y , while for a large q it amplifies the enhancement effect. Thus, with such a diffusivity, one can achieve fast particle transport by changing the channel shape for given external force.
Results for the isothermal case are provided in figure A4 for comparison. A major difference is that with increasing σ, the amplitude variation of v is relatively small compared to the Itô case. For small F, with different q, v first decreases and then increases, while in figure 7(c) it will decrease. A detailed analysis is provided in the supplementary figures in appendix A.
In summary, the influence of the parameters F, q and σ is studied by calculating the mean velocity v for different D x, y . Contrasting to D 0 , v of both D d x, y and D m x, y can become larger. For the case D d x, y , a small σ is detrimental for particles to pass through the bottleneck, while for D m x, y , a small σ promotes the process.

Mean first exit time (MFET) and first exit time density (FETD)
In section 4, we explored the influence of D x, y on the mean velocity. Because particles may move to the left or to the right, v measures the overall speed of the particles. It is hard to deduce from v how precisely particles exit a cell. Thus, in this section, we study the time that particles first pass through a bottleneck from either side. One of the most powerful measurements for such times is the FETD, which contains rich information, such as the MFET, the variance and the most probable exit time. Numerical and analytical techniques have been proposed to study the FETD in an open environment [72][73][74][75]. Generally, the FETD p (τ ) can be derived by counting the number of exit times that fall in a series of intervals. To do this, we obtain a large sample set of the first exit time as where the set T contains the times that particles first exit the cell with initial position at (0, 0). Then, counting the number of times locating in different intervals ( τ k − δt/2, τ k + δt/2 ], and dividing the total number and length of the interval, we obtain the discrete FETD in which T j is the jth element in set T. δt is the statistical interval. N is a sufficiently large integer to ensure sufficient statistics. I {·} is the indicator function. Based on the FETD, we also obtain the MFET by the expectation formula A given particle may exit from either side of a cell. Therefore, we also take the splitting probability (SP) into account to analyse the exiting process. For simplicity, we calculate the SP to the right (the SP to the left equals to one minus the SP to the right), which can be obtained by All symbols have the same meaning with equation (11). Firstly, the MFET for D x, y = D 0 is presented. We uncover that the MFET decreases as F increases. Namely, a larger F leads to shorter times particles take to reach the bottleneck. Additionally, it shows that for a very small F the MFET is almost symmetric around q = 0; for large F, the MFET increases with q. When F = 0, particles have equal probability to move to right and left for q = 0 (see the inset in figure 8(c)). For q > 0, the cell is right-skewed and has a wide space at the right side. Thus, such a channel shape is convenient for particles to diffuse to the right and then exit. This can be verified by the large SP for q > 0 shown in figure 8(c). Similarly, for q < 0, the cell has a large left-hand side space, which is easy for particles to pass through the left-hand side bottleneck, and induces a small SP to the right. Thus, the MFET is exactly symmetric in dependence on q when F = 0. From figure 8(c), we can clearly see the FETDs for different q. For large |q|, the FETD has a higher peak, which indicates a thinner tail and a smaller MFET. This coincides with figure 8(a).
When F is relatively large, particles are forced to move to the right. In this case particles will likely exit from the right-hand side bottleneck, rather than the left one as in the case of F = 0. The large SP in the inset in figure 8(b) verifies this intuition. Moreover, with increasing F, the SP increases and reaches its maximum at SP = 1, which implies that all particles consistently exit from the right-hand side bottleneck. Figure 8(b) shows that with increasing F the FETD is concentrated increasingly at small t, namely, the peak grows progressively higher. This corresponds to the small MFET in figure 8(a). In addition, for a large F, with increasing q the MFET increases. As stated above, a large q will induce a wide space at the rhs of the cell. Although this makes it convenient to diffuse to the right, a steep right-hand side wall will make it harder than a flat wall for particles to locate the bottleneck and exit. Therefore, the MFET increases monotonously with q for large F.
When D x, y = D d x, y , the variations of the MFET with the system parameters F, q and σ are shown in figure 9. In figure 9(a), although the increase of F still results in a reduction of the MFET, the magnitude of the MFET changes dramatically. On the one hand, when F is small, the MFET of D 0 is around 3. However, for D d x, y , the MFET increases to 60. On the other hand, with increasing F, the MFET reduces sharply rather than the slowly decreasing behaviour in figure 8. The influence of F on the FETD and SP are similar to those of D 0 , so the figure is not presented here.
In addition, there are some other interesting phenomena to report. In figure 9(b), when F = 0, the MFET first increases for 0 σ 1 but decreases for σ > 1. The MFET attains its maximum at σ = 1. To get a clear insight into the distribution of first exit times, we present the FETD with respect to σ for F = 0 in figure 9(f). When σ = 0 the FETD has a maximum very close to t = 0 and a relatively thin tail. With increasing σ the maximum becomes smaller and its location moves away from t = 0. In addition, the tail of the FETD is very broad. Namely, there exist extremely long first exit times, which induces a large MFET. This phenomenon corresponds to region the 0 σ 1 in figure 9(b). With further increase of σ the maximum of the FETD increases, which indicates progressively shorter first exit times. Thus, the MFET decreases, corresponding to the region σ > 1 in figure 9(b). Moreover, when F > 0 in figure 9(b) the MFET decreases monotonously. This implies that σ shows different degrees of impacts for F = 0 and F = 0. Considering the isothermal case, when F = 0 in figure A5(b) the MFET first increases fast and then remains almost unchanged with increasing σ, which is very different from the Itô case. Specially, for σ = 0, the MFET is much smaller than the Itô case, however, for large σ, the MFET is much larger. Additionally, when F is small the MFET does not always decrease with σ. Conversely, the MFET will first increase, then decrease and then increase again, resulting in a minimal MFET at σ = 1. This is consistent with figure A3.
To further explore the interplay of q and σ, the MFETs for F = 0 and F = 0.1 are displayed in figures 9(c) and (d). When F = 0, for a small |q| the MFET first increases and then decreases with σ, leading to a maximum for each q (for q = 0 figure 9(f) shows the corresponding FETD). However, for extremely large |q| the MFET increases almost monotonically. Moreover, the σ of the location of the maximum increases with |q|. It is worth noting that for q = 0 the location is at σ = 1, which is the smallest value for different q. In fact, for q = 0, the diffusion coefficient becomes x, y is the smallest for all σ. For |q| = 0, the σ of the maximal D σ d x, y changes with the channel shapes, i.e., it is increasing with |q|. All these factors result in the quite complex graph in figure 9(c). We present the typical FETD for σ = 1 in figure 9(e), in which the maximum of the FETD increases with |q|. Namely, the FETD concentrates more at small t for large |q| than for small |q|. Conversely, for σ > 1, the maximum of the FETD will decrease with |q|, corresponding to figure 9(c) (not shown here). In the isothermal case, except for the extremely short MFET for small σ, we find that a smaller |q| will trap particles longer, and this is opposite to the Itô case, see figure A5(c). Moreover, the MFET generally increases with σ for large |q|.
When F = 0.1 in figure 9(d), the situation becomes totally different compared to figure 9(c). It shows that for σ = 0 the MFET decreases with reducing q; for large σ, in contrast the MFET first decreases and then increases. Namely, the superposition of F eliminates the maximums of MFET for F = 0, which indicates that F is able to weaken the influence of σ. This phenomenon is verified in figure 9(b), in which for a large F the MFET is almost unchanged for different σ. In other words, the influence of σ is very sensitive to F, such that already small F values can effect a dramatic change. In the isothermal case, we provide figures A5(d) and (e) for F = 0.1 and F = 1 to show the variation of MFET with respect to different force strengths.
Next, we consider D x, y = D m x, y . Figure 10(a) shows that the overall trend is similar to figure 9(a), but the MFET is much larger than for D 0 , and smaller than for D d x, y . Moreover, the MFET varies slightly with respect to q for F > 0.
In figure 10(b), for F = 0, the MFET first decreases and then increases with σ. This can be verified from the FETD in figure 10(f). For σ = 0, the tail of the FETD is broad, and its maximum is relatively small, which induces a large MFET. With increasing σ the maximum of the FETD first increases and then decreases (see the bold red and yellow lines). Correspondingly, the tail first gets thinner and then becomes broader. Thus, we obtain the shape of the MFET marked by the red line in figure 10(b). In addition, for F > 0 in figure 10(b) the MFET increases monotonously. Namely, a small change of F will lead to a significant difference of the MFET. We note that when σ = 0 the MFET shows a quite dramatic descent with only a tiny increase of F = 0 (see the red and green dots, the inset shows the large differences of the FETDs for F = 0 and F = 0.5). This indicates that the MFET is very sensitive to F at σ = 0. An abrupt decrease of the MFET is also observed for F = 0 in the isothermal case-but with further increase of σ the MFET almost stay unchanged in figure A6(b) rather than increasing, as in figure 10(b). Finally, for relatively small F a large σ value will induce a small MFET, but for large F and σ = 0 it is more likely to effect the smallest MFET.
In figure 10(c), the behaviours of the MFET is almost opposite to figure 9(c). The MFET first decreases and then increases, and it decreases faster for small |q|. For moderate |q|, there is a minimal MFET, located at σ = 1 for |q| = 0. The location increases with |q| (for extremely large |q|, the minimum is located at small σ). In fact, for q = 0, the diffusion coefficient becomes D σ=0 x, y . Thus, we get the minimal MFET at σ = 1. For |q| = 0, the corresponding σ increases, and so does the position of the minimum. The corresponding FETD with respect to q is displayed in figure 10(e), which shows a similar property with figure 9(e), but the maxima have a small amplitude variation. In the isothermal case, figure A6(c) shows that the MFET for σ = 0 is much larger than that in figure 10(c). This means the noise-induced drift will trap particles at the central area and play a negative role to the diffusion of particles. Without external forces, a large σ is preferred to achieve fast exiting.
When F = 0.1 in figure 10(d), except for σ = 0 the variation of the MFET coincides with figure 10(c). For σ = 0 a very small F will induce a dramatic reduction of the MFET. In fact, when σ = 0 the shape of D m x, y appears like a horizonal groove along the x axis. Such a structure is greatly beneficial for particles transport under external forces. The reason is that particles are trapped near the x axis, which greatly avoids collisions with channel walls. On driving by an external force, the number of collisions will be further reduced, and particles will be rectified fast and significantly accelerated. In this sense, D σ=0 m x, y is the best condition for particle transport. In the isothermal case, figure A6(b) shows the different influence of small and large F on the MFET. Here we provide two figures illustrating the differences in the q-σ plane with respect to F = 0.1 and F = 1 in figures A6(d) and (e).
In this section, the MFET of particles is investigated. The MFETs for different cases indicate that the diffusivity of the environment has a great influence on the effective transport of particles. Generally, compared to the constant diffusivity D 0 both D d x, y and D m x, y increase the difficulty for particles to exit. However, the effects are very sensitive to external forces, with which the MFET changes dramatically, and features non-trivial maxima and minima for different parameters.

Conclusions
The transport of particles in a confined channel with three scenarios for position-dependent diffusion coefficients is studied in this work. The PDF, the mean velocity, and the MFET are calculated to illustrate the influences of these D x, y . We find that, in the absence of external forces, for a constant diffusion coefficient D 0 , we obtain the expected homogeneous PDF of particles. For the diffusivity D d x, y corresponding to a heat flow or mobility gradient, the PDF mainly focuses near the channel walls. By adjusting its dependence on the x coordinate, the preferred gathering position of particles changes from vertical corners to the bottleneck. Conversely, for the diffusivity D m x, y corresponding to a cold flow, the PDF has a peak at the cell centre. By modulating its dependence on the x coordinate, the contour plot of the PDF changes from a horizonal ellipse-like shape to a vertical ellipse-like shape. In the presence of external forces, the PDFs change significantly.
For the mean particle velocity, both D m x, y and D d x, y can lead to a larger velocity than for a constant D 0 . For the case D d x, y , a high dependence on the x coordinate will accelerate the transport speed. In contrast, for D m x, y , a low dependence on the x coordinate will remarkably enhance the transport. Moreover, a large external force can weaken the influence of the position-dependent diffusivity. In addition, we find that although D m x, y and D d x, y increase the mean velocity compared to D 0 , particles take longer to exit a cell from the centre. Generally, D d x, y can trap particles within a cell longer than for the other two cases under the same external forces and channel shapes.
From our results, we conclude that in the presence of external forces, the diffusivity of D m x, y is more beneficial for particles to move with a high velocity in periodic channels. In particular, a zero dependence on the x coordinate will extremely promote the process. The results can be applied to accelerate particle transport by adding a cold fine tube along the central line of two-and three-dimensional channels in real systems, or by introducing a gel with modulated porosity.
However, we also point out that the results in the isothermal interpretation show considerable differences when the external force is zero or small, while for large external forces the dynamics are almost the same. It still needs continued effort from both experiments and theories to confirm the applications of Itô and isothermal integral for given systems. Generally, the particles still tend to stay at low diffusivity areas, for example, the upper and lower corner in (a1), all around in (a2) and left and right sides in (a3). However, the PDFs are almost uniform such that the differences in the colourmaps are very small. Moreover, we see that the marginal PDFs vary slightly. (b) F = 5.0, q = 0. As particles can distribute almost uniformly in the channel, they do not steer clear of the central area as seen in figure 3(b). Further, more particles gather at the right bottleneck for σ = 0.1 than in the other cases. (c) F = 5.0, q = 0.4. We see that more particles crowd at the right bottleneck than for q = 0, which means that the steep wall on the right is not beneficial for particles to cross the bottleneck.

Appendix A
Here, we provide some additional figures discussed in the main text. Generally, the PDFs are almost uniform distributions such that the differences in the colourmaps are very small. But we also find that the central area has a higher probability. Additionally, we see that the marginal PDFs are almost the same. (b) F = 5.0, q = 0. For σ = 0.1, more particles distribute along the central line than in the other cases. For larger σ the probability is higher at the right side. (c) F = 5.0, q = 0.4. For a right-skewed shape more particles crowd at the right bottleneck than for q = 0. This means that a steep right side wall makes it difficult for particles to locate the bottleneck and prevents them from crossing to the next cell. Figure A3. v for D d x, y in the isothermal interpretation. (a) σ = 0.1. With increasing F, v almost increases linearly for different q. But v increases faster for a small q than for large q. (b) q = 0. For large F, with increasing σ, v first increases and then hardly varies. But for a small F, v first decreases, then increases and then decreases. (c) F = 1.0. We see that v first decreases and then increases with σ for very small q. When q is around zero, v first increases and then decreases. However, for large q, v always decreases. (d) F = 10.0. For large F, the effects of the noise-induced drift are very limited. v will increase with σ when it is small. Generally speaking, σ ≈ 1 is an optimal choice to achieve a large v . Figure A4. v for D m x, y in the isothermal interpretation. (a) σ = 0.1. For a fixed q, v increases monotonously with F, and it increases faster for a smaller q. (b) q = 0. v will first decrease, then increase with σ, and remain almost constant for large σ. (c) F = 1.0. For fixed σ, with increasing q, v will decrease monotonously. For fixed q, v will first decrease and then increase with σ. (d) F = 10.0. When |q| is large, v will decrease monotonically. But when |q| is small, v will first decrease and then increase. Finally, for a fixed σ, v will increase with decreasing q. Figure A5. MFETs for D d x, y in the isothermal interpretation. (a) σ = 0.1 when F = 0, the MFET decreases with |q|. For large F, the MFET will increase with q. (b) q = 0. For F = 0, the MFET will first increase fast with σ, and then it will remain almost constant. For small F = 0, the MFET will first increase, then decrease and then increase, resulting in minimal MFET at σ = 1. This is consistent with figure A3. For large F, the MFET will first decrease and then increase. (c) F = 0. The MFET is much smaller for a small σ than for a large σ. For fixed q, with increasing σ the MFET first increases fast and then levels off. (d) F = 0.1. Due to the interactions with the noise-induced drift, the variation of MFET with respect to q and σ is complex. It will further change with increasing F. However, we see that a small change for F will cause a large reduction for the MFET. (e) F = 1. With a larger F than that of (d), the MFET becomes totally different. For large q, the MFET increases almost monotonically with σ, however, for a q > 0.9, the MFET first increases, then decreases and increases. Finally, for a fixed σ, the MFET will decrease with decreasing q. Figure A6. MFETs for D m x, y in the isothermal interpretation. (a) σ = 0.1 when F = 0, the MFET decreases with |q|. For large F, the MFET will increase with q. (b) q = 0. For F = 0, the MFET has an extremely large value at σ = 0. With increasing σ, the MFET will decrease fast and then remain almost constant. For a relatively small F, the MFET will first increase and then decrease, but for a large F, the MFET will increase and then level off, where σ = 0 induces the smallest MFET. (c) F = 0. When |q| is large, the MFET decreases monotonously, but for small |q|, the MFET first decreases and then increases. (d) F = 0.1. A small external force will lead to an abrupt reduction of MFET for small σ, while it decreases less for a large σ. (e) F = 1. With a larger F than (d), the MFET becomes totally different. For large q, the MFET firstly increases abruptly and then remains almost constant with σ. For small q, the variation of σ impacts slightly on the MFET. Generally, for a fixed σ, with increasing q, the MFET increases monotonously.