{'fraction_non_alphanumeric': 0.04255700284054559, 'fraction_numerical': 0.013477663757879742, 'mean_word_length': 4.127361791462561, 'pattern_counts': {'":': 1, '<': 31, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Monte Carlo Approaches for Simulating a Particle at a Diffusivity Interface and the"Ito--Stratonovich Dilemma"

The possibility of different interpretations of the stochastic term (or calculi) in the overdamped Langevin equation for the motion of a particle in an inhomogeneous medium is often referred to as the"Ito--Stratonovich dilemma,"although there is, in fact, a continuum of choices. We introduce two Monte Carlo (MC) simulation approaches for studying such systems, with both approaches giving the choice between different calculi (in particular, Ito, Stratonovich, and"isothermal"). To demonstrate these approaches, we study the diffusion on a 1D interval of a particle released at an interface (in the middle of the system) between two media where this particle has different diffusivities (e.g., two fluids with different viscosities). For reflecting boundary conditions at the ends of the interval, a discontinuity at the interface in the stationary-state particle distribution is found, except for the isothermal case, as expected. We also study the first-passage problem using absorbing boundary conditions. Good agreement is found when comparing the MC approaches to theory as well as Brownian and Langevin dynamics simulations. Also, the results themselves turn out to be interesting. For instance: 1) for some calculi, there are more particles on the low-viscosity side (LVS) at earlier times and then fewer particles at later times; 2) there is no preference to end up on a particular wall for the Ito variant, but a bias towards the wall on the LVS in all other cases; 3) the mean first-passage time to the wall on the LVS grows as the viscosity on the other side is increased, except for the isothermal case where it approaches a constant; 4) when the viscosity ratio is high, the first-passage-time distribution for the wall on the LVS is much broader than for the other wall, with a power-law dependence of the former in a certain time interval whose exponent depends on the calculus.


I. INTRODUCTION
While deceptively simple to consider conceptually, the diffusion of a particle in an inhomogeneous medium is a subtle problem [1][2][3][4]. For a system consisting of a high viscosity and a low viscosity region, will a tracer particle be found more often on the "thicker" side, the "thinner" side, or with an equal probability at all locations? While an assumption of a Boltzmann distribution would suggest an equal probability at all locations, a simple physical argument suggests that the particle may be more likely to be found on the high viscosity side. If the only inhomogeneity in the system is the viscosity, then statistically, the dynamics are the same at all points with the exception that time is rescaled on the high viscosity side where the particle moves slower. Hence the particle will spend a greater amount of time on that side and should be more likely to be found there.
While this discussion can go back and forth, in truth each answer corresponds to a different formulation of the problem -both with their own mathematical and physical assumptions. In essence, the question boils down to considering a particle at the viscosity interface: will there be an additional "push" at the interface to the low viscosity side that will compensate the "trapping" in the high viscosity region due to the slower dynamics there? Let us assume that the fluids can be considered implicitly, so their effect on the diffusing particle is reduced to a random force and a friction force that the particle feels. We further assume that the interface is sharp and that there are no interfacial effects such as surface tension that would make the particle prefer or avoid the interface. Finally, we assume the particle motion is overdamped such that it can be described by the overdamped Langevin equation [5]. The situation where the assumption of equivalent dynamics (apart from the time scale) holds and particles accumulate on the high viscosity side corresponds to interpreting the stochastic term in the overdamped Langevin equation according to the rules of the Ito calculus [6]. In other interpretations, of which the most well-known corresponds to Stratonovich calculus [7], but perhaps of most interest to physicists is the so-called isothermal formulation [8][9][10][11], the particle "feels" the interface and has a tendency to move to the low viscosity side. In particular, in the isothermal case (also referred to in the literature as Hänggi-Klimontovich or kinetic [12]) this bias at the interface cancels the trapping effect and, in the stationary state, the particle has an equal probability of being found at all locations.
There have been several experimental explorations of systems with nonuniform viscosity or, more generally, particle diffusivity [1,2]. Most relevantly, van Milligen et al. [1] performed experiments in which the relative concentration of dye was monitored in a system comprised of two viscosity regions created by adding gelatine to one half of it. In this setup, the dye was found to accumulate on the high viscosity side. As demonstrated in their paper, the concentration distributions they observed can be approximated by solutions of the Fokker-Planck equation, which corresponds to the Ito formulation (see Sec. II). On the other hand, in experiments by Lançon et al. [2], particles diffused between two nearly, but not perfectly, parallel walls, which, as the authors showed, is equivalent to diffusion with a variable diffusion coefficient. In contrast to Ref. [1], while individual particles were observed to drift to regions of higher diffusivity, there was no overall flux when starting with a uniform concentration of particles and hence the uniform concentration was preserved, which is consistent with the isothermal case. On the theoretical side, the general conclusion has been that different formulations can be appropriate depending on the details of the system. For instance, in a simple model of diffusion in a slowly modulated periodic potential by Sokolov [12], Ito, Stratonovich, and isothermal rules correspond to different variants of the model.
Given the experimental and theoretical relevance of diffusion in inhomogeneous media, it is desirable to be able to model it computationally. One popular computational approach to modeling the dynamics of particles in a fluid medium is molecular dynamics (MD) [13]. Two of the most common MD approaches assuming an implicit fluid are Brownian dynamics (BD) and Langevin dynamics (LD). The BD approach consists in solving numerically the overdamped Langevin equation in the Ito formulation. LD drops the overdamped assumption and adds an inertial term to the equation of motion that is then solved numerically. Interestingly, even in the limit of an infinitely large damping (but decreasing the time step so it is always much smaller than the corresponding relaxation time), LD is not reduced to BD and instead corresponds to the isothermal formulation of the overdamped Langevin equation. Thus the choice between the two methods is not arbitrary when an inhomogeneous system is considered.
Monte Carlo (MC) methods are another frequently used group of computational approaches [14]. Unlike MD, MC methods generally sacrifice at least some aspects of the dynamics of the system to achieve a computational speedup compared to MD. How much is sacrificed can vary. The popular Metropolis algorithm [15] is designed solely to reproduce the equilibrium state and is generally unsuitable for studying the dynamics, even qualitatively. For instance, in the problem of inhomogeneous diffusion considered here, a particle would diffuse equally rapidly regardless of the viscosity, unless care is taken to adjust the step length and/or the time increment at every MC step (the time step) appropriately, which is not a part of the basic algorithm. However, there are also methods (dynamical MC ) which do reproduce the principal aspects of the dynamics and even the correct time scale of the process. This includes generally faster, but sometimes less accurate and/or reliable methods where the particles are restricted to pre-defined points in space forming a lattice (lattice MC, or LMC), as well as off-lattice methods. Dynamical LMC methods include kinetic Monte Carlo [16][17][18] and methods designed to reproduce the correct dynamics in an arbitrarily strong uniform field [19][20][21].
In this article, we present two dynamical LMC algorithms which can be used to generate results consistent with different formulations of the inhomogeneous diffusion problem. We consider a simple system consisting of a single random walker diffusing in a box which contains two separated fluids in equal proportions: a fluid of viscosity η L on the left side and viscosity η R on the right side, similar to the experimental setup in Ref. [1]. Hence, there is a viscosity interface at the center of the box. Projecting the motion of the random walker onto the direction perpendicular to the interface reduces the problem to a 1D one. The particular choice of the problem with a sharp interface has been made for two (in a way, opposite) reasons. On the one hand, a priori it is the simplest setup in which to study diffusion in an inhomogeneous medium; it lets us define simple quantities characterizing the process, such as the fraction of the time spent on a particular side of the interface. On the other hand, when it comes to actually solving the problem, it turns out there are difficulties specifically in the case of a sharp interface, as discussed in Sec. II, so in that sense this problem is, in fact, more complicated than that of a gradual change in viscosity; our algorithms should generalize straightforwardly to the latter case. In addition, in future work we intend to apply the algorithms developed here to studying translocation of a polymer through a nanopore in a thin membrane, in which case having different viscosities on the two sides of the membrane with a sharp boundary between them corresponding to the pore is natural. Overall, the system of a single random walker in a fluid with a single, sharp viscosity interface is chosen as a simple "toy problem" such that we can clearly demonstrate the validity of our LMC algorithms. However, as will be shown in the results section, many interesting features arise for the dynamics of a particle in such a system -particularly when contrasting results for different calculi. This manuscript thus not only introduces LMC approaches for simulating diffusion in a system containing regions of differing diffusivity, but it demonstrates the complexity which arises for even the simplest case and highlights the importance of choosing the calculus that is appropriate for the problem under study.
In the next section, we briefly review the relevant parts of the theory of the overdamped Langevin equation for the motion of a Brownian particle and obtain the corresponding equation for the particle concentration. We then introduce Ito, Stratonovich, and isothermal calculi giving some examples of their physical relevance. We also point out some problems with using the overdamped Langevin equation in the case of a sharp interface. In the following two sections, we consider different simulation techniques that can be used for the problem of particle diffusion. We first note that the well-known BD and LD methods correspond to Ito and isothermal calculi, respectively, and then introduce two MC algorithms. The details of the simulation protocol are outlined in Sec. V, after which in Sec. VI we present comparisons between different methods (BD, LD, MC) and results of the simulations. We end the paper with a discussion of possible uses of the methods developed here and their future generalizations.

A. Diffusion in a liquid
Diffusion of a pointlike particle along the x axis under the influence of a potential U (x) can be described by the Langevin equation [22], where m is the particle mass, k is the Boltzmann constant, T is the temperature, and R(t) is a random function that satisfies R(t) = 0 and R(t)R(t ′ ) = δ(t − t ′ ), where . . . denotes the ensemble average. The coefficient of the last term follows from the fluctuation-dissipation theorem. It is assumed that the drag force F f (the first term on the right-hand side) is proportional to the particle velocity v, where ζ is the friction coefficient. For a particle moving in a fluid, ζ is proportional to the fluid viscosity η: where κ is constant for a given particle. For instance, for a spherical particle of radius R, according to the Stokes' law (which assumes no-slip boundary conditions for the fluid), for other particle shapes and boundary conditions the expression for κ will be different, but Eq. (3) is still valid. In many cases of interest, the left-hand side (or the inertial term) in Eq. (1) is small. If this term is neglected (the overdamped limit ), we obtain the overdamped Langevin equation It is this limit that we will consider in what follows. However, neglecting the inertial term can bring about subtle problems, as discussed below. Taking the case of free diffusion [U (x) = 0] we can rewrite Eq. (5) as The particle diffusion coefficient (or diffusivity) D is defined by the expression where ∆x is the particle displacement over time ∆t. Based on Eq. (6), it can be shown that We then obtain B. Inhomogeneous media The consideration below and in Section II C is similar to Ref. [2], but somewhat more detailed.
If we consider diffusion in a stationary inhomogeneous medium, D will be dependent on x such that To proceed, we must determine how the multiplicative noise term in Eq. (10) is to be evaluated. For a deterministic differential equation where f (x, t) is a nonrandom and well-behaved function, we can write down the derivative on the left-hand side as the limit of a difference and then In the last term, instead of the value of the function f at time t, we could have taken its value at time t+δt, or any linear combination of the two -this does not matter in the limit δt → 0, i.e., in this limit the equation can be used with any α between zero and one. Varying α changes the right-hand side of Eq. (13) by a negligible amount O(δt 2 ). However, for stochastic differential equations, like Eqs. (9) and (10), there are complications.
First, R(t) is not a well-behaved function. Instead of its instantaneous value, which is not well-defined, one needs to take its average over the interval δt, which depends on δt and is r/ √ δt, where r is a random number with r = 0, r 2 = 1 and no correlation between different time intervals. Then, for Eq. (9), Second, when D is variable [as in Eq. (10)], there is an additional uncertainty, since the result depends on where D is evaluated, that is, on the value of α in (15) and this dependence now remains even in the limit δt → 0. Indeed, let us assume for now that D(x) is a smooth function. From now until the end of this section, we will use x(t + δt) to denote the value calculated using Eq. (15), rather than the actual value of the coordinate at time t + δt, thus, the "≈" sign in Eq. (15) is replaced by the "=" sign. Expanding and using this in Eq. (15), we get In Eqs. (16), (17) and below, we have not specified explicitly the point at which the derivative dD/dx is taken. It is easy to check that replacing the value of dD/dx at the point x(t) with that at the point x(t + δt) or vice versa introduces an extra term of order δt 3/2 in the last line of Eq. (17) that can be neglected, so the exact point at which the derivative is taken does not matter. On the other hand, both the first and the second term in the last line of Eq. (17) are important in the limit δt → 0, even though they are of different orders in δt, because r = 0, but r 2 = 1 = 0. Indeed, over a fixed time interval ∆t such that ∆t = N δt, the contribution of the first term as a function of N and δt is ∝ √ N δt = √ ∆t, while that of the second term is ∝ N δt = ∆t, and both terms survive when δt → 0. In the second term, we can replace r 2 with its average value 1, since the contribution of the fluctuating part over the interval ∆t is ∝ √ N δt = √ ∆tδt, which is negligible when δt → 0. This gives x(t + δt) ≈ x(t) + r 2D(x(t))δt + α(dD/dx)δt. (18) It is easy to show that the same result is obtained, if in Eq. (15) the value of D is taken at a single point between the end points of the jump, as in Ref. [2], that is, Eq. (15) is replaced with The last term in Eq. (18) is deterministic and produces drift of the particles in the direction of the diffusivity gradient known as "spurious flow" [23], "spurious drift", or "noise-induced drift" [5]. The drift velocity is it vanishes when α = 0. Note that the same drift can be produced by a suitably chosen external deterministic force, so it is always possible to change the value of α and still get the same dynamics by introducing an additional fictitious force [12].

C. Particle fluxes and continuum diffusion equations
Although the drift velocity of particles is zero for α = 0 [Eq. (20)], the net flux through a fixed point x 0 does not vanish in this case, even without a particle concentration (or probability density) gradient. Indeed, assuming for simplicity that r = ±1 (with no influence on the final result), in a single jump the point with x = x 0 is crossed from the left by one half of all particles (namely, those with r = +1) with x between x 1 and x 0 , where x 1 obeys the equation or with the solution The number of such particles is where ρ is the (uniform) particle concentration. Likewise, the same point x 0 is crossed from the right by one half of all particles (those with r = −1) with x between x 0 and x 2 , with and the number of such particles is The net flux is thus When α = 0, the flux due to the drift, ρV , with V given by Eq. (20), has to be added to Eq. (27), which gives the total flux This flux is nonzero when the drift velocity is zero at α = 0, but it vanishes when α = 1. When a concentration gradient is present [ρ = ρ(x, t)], this adds the usual flux given by Fick's first law, −D ∂ρ ∂x : In two special cases this expression is simplified: for α = 0, while for α = 1, which conforms to normal Fickian diffusion rules. The rate of the concentration change is In particular, for α = 0, which is the standard form of the Fokker-Planck equation without the drift term [23]; for α = 1, which is Fick's second law. The stationary state can be obtained by setting ∂ρ/∂t = 0 in Eq. (32). This condition for no net change in the particle concentration is then If there is no net flux, as, for example, when the system is confined between two reflecting walls, then the constant in Eq. (35) is zero and which gives For α = 1, the particles can be found with equal probability anywhere (ρ = const), but for α < 1 they are more likely to be found in regions of lower diffusivity (or higher viscosity).

D. Ito, Stratonovich, and isothermal calculi
Even though the actual solution of Eq. (10) depends on the value of α in its discretized variant [Eq. (15)], its formal solution can always be written as Different α thus correspond to different interpretations of the formal integral in Eq. (38), or different calculi. Three cases have received special attention. The case α = 0 corresponds to Ito calculus [6]. In this case, the stationary distribution is given by The primary attractive quality of this formulation is that its implementation is "non-predictive": According to Eq. (15), the value of the diffusivity D is taken at the initial point of the jump only, so, as opposed to all other cases, one does not need to look ahead in time to evaluate the spatially dependent term. This property makes the Ito formulation a frequent choice for mathematical problems -including financial modeling. But it is also physically relevant: For instance, the overdamped limit of the Langevin equation (1) corresponds to the Ito calculus in the idealized situation where the temperature varies in space, but the friction coefficient is temperatureindependent and thus the same everywhere [24]. This is an inherently non-equilibrium situation from the thermodynamic point of view and so the fact that the stationary distribution does not coincide with the equilibrium Boltzmann distribution is not surprising. Another example would be a particle in a modulated quasiperiodic potential where the tops of the barriers are all at the same level, but the well depths are position-dependent [12]. In this case, the change in diffusivity is accompanied by the change in the free energy and the Boltzmann distribution is obeyed. The case α = 1/2 corresponds to the Stratonovich calculus [7]. In this case, the value of D at the midpoint of the jump (or the average of the values at the endpoints) is taken. Stratonovich calculus conforms to the "normal" rules of calculus. For example, if x obeys Eq. (10), then one would expect where the right-hand side is interpreted according to the rules of the given calculus. This is indeed true for Stratonovich, but not for other calculi. The Stratonovich calculus arises when considering the Langevin equation (1) with the noise R(t) having a finite and space-and time-independent correlation time, in the limit when both this correlation time and the relaxation time m/ζ go to zero, but the former is kept much larger than the latter [25]. In fact, even for the overdamped Langevin equation (10) (corresponding to the relaxation time equal to zero) with a finite noise correlation time the issue of interpretation of the noise term does not arise and the limit when the correlation time goes to zero corresponds unambiguously to the Stratonovich calculus [5,23]. This situation violates the fluctuation-dissipation theorem which requires that a finite noise correlation time be accompanied by a memory in the friction force, thus giving rise to the generalized Langevin equation [22] instead of the ordinary one. So the fact that the density in the stationary state is not uniform, but rather obeys is natural. Sokolov [12] has also proposed a variant of his modulated quasiperiodic potential model that corresponds to the Stratonovich calculus; again, the free energy in that model is position-dependent. Since Ito and Stratonovich are the two calculi most often discussed in the literature, the choice between different calculi has often been referred to as the "Ito-Stratonovich dilemma" [23]. But besides these two cases, the case α = 1 is also of special interest, because in that case Eq. (37) becomes which corresponds to the equilibrium Boltzmann distribution in the "purest" case when it is only the diffusivity that varies in space, but the particle free energy and the temperature are the same everywhere. There is no universally accepted name for the α = 1 calculus, with names such as Hänggi-Klimontovich [12], kinetic [11,12] and isothermal [2] appearing in the literature. We will use the latter in what follows, since it emphasizes the fact that this calculus corresponds to the case when the system is kept at a constant temperature, as opposed to the Ito case which can be achieved by varying the temperature across the system (still assuming the situation where the free energy would be constant for T = const).

E. Discontinuous D(x)
So far, we have assumed that D(x) varies smoothly in space. First of all, this allowed us to do the expansion in Eq. (16). Even more importantly, we note that if there is a jump in the diffusivity, then even the solution of Eq. (10) is not well-defined for α = 0. Consider α = 1 and let us assume that with D L > D R . Let the particle coordinate at time t be x 0 < 0 and suppose that the time step δt and the random number r chosen at the following step are such that x 0 + √ 2D L rδt > 0, but x 0 + √ 2D R rδt < 0. Note that the jump size depends on the final position at time t + δt. If we assume that the particle ends up to the right of the interface, then the jump size should be √ 2D R rδt, but then the particle coordinate after the jump is x 0 + √ 2D R rδt < 0, so we come to a contradiction. Conversely, if we assume that the particle ends up on the left side, then the jump size should be √ 2D L rδt, but x 0 + √ 2D L rδt > 0, so there is again a contradiction. Thus at this point the particle behavior is undefined. On the other hand, for x 0 > 0 it is possible that both assumptions lead to valid results, so there is ambiguity. Given that in principle |x 0 | can be arbitrarily close to zero, this will always be a potential problem, no matter how small δt is. In fact, this problem exists whenever the jump size depends not just on the initial position, that is, for any α = 0. Thus, in order to be able to use Eq. (10) [discretizing it as in Eq. (15)] for calculi other than Ito (and even for the very concept of different calculi to be meaningful), we need to assume that the function D(x) is smooth everywhere. In that case, what we mean by a "sharp" interface and a "jump" in diffusivity is that actually the diffusivity is smooth and thus the interface between two liquids has a finite width, but the interface is narrow compared to all other length scales (in our simple problem, there is only one such length scale, the size of the system). Physically, the smoothness of D means that it changes little on the length scale of the typical interatomic distance. This is very likely to be true in reality, and even if the liquids are completely immiscible and the interface is atomically sharp, but the diffusing particle is much larger than the interatomic distance, the properties of the system are effectively averaged out over the particle size and the condition is satisfied. In this case, Eq. (10) is still meaningful for all calculi and in principle can be solved numerically by choosing a very small time step when iterating Eq. (15), so that typical particle jumps in the simulation are smaller than the width of the interface. However, this would be a waste of computational resources; but choosing a larger step would lead to the indeterminacy problem described above. This is where LMC methods are helpful, as we will see below.
In the remainder of this paper, we will concentrate on this situation of a sharp, but not infinitely sharp interface. For simplicity, we will still use notation that may seem to imply that the interface is infinitely sharp, but will keep in mind that in reality it is not. So, for instance, ρ(+0, t) denotes the particle concentration immediately to the right of the interface, but still far enough from it that the diffusivity has already reached its asymptotic value D R ; likewise, D(+0) ≡ D R . Bearing this in mind, Eq. (37) for the stationary distribution and its particular cases, Eqs. (39), (41), and (42), should remain valid in our situation. Moreover, we can show that these relations are valid at all times, even in the transient state, in the vicinity of the interface, that is, at all times Indeed, with the exception of the case α = 1, the first term on the right-hand side of Eq. (29) for the particle flux is large at the interface and in fact infinite in the limit of an infinitely sharp interface. However, the flux itself, even though it can be nonzero when the stationary state is not reached, is expected to be finite. Therefore, the second term on the right-hand side of Eq. (29) should be large (infinite in the limit) as well, and the flux J can be neglected compared to either of these terms. Dropping J in Eq. (29) produces the stationary-state equation (36), from which Eq. (44) follows immediately. In the special case α = 1, the first term on the right-hand side of Eq. (29) is zero and therefore the second term is finite, thus there is no jump in ρ, which is likewise consistent with Eq. (44). Condition (44) describes the concentration jump across the interface that forms immediately and always exists, if α = 1. The full solution of the diffusion problem for the system consisting of two uniform regions separated by a "sharp" interface (where "sharp" is understood as above) can be obtained by solving two uniform diffusion problems in the two regions, with two additional matching conditions at the interface: Eq. (44) and the flux continuity condition, In addition, of course, there have to be boundary conditions at the walls, for example, the no-flux condition for a reflecting boundary or for an absorbing boundary, where x b is the coordinate of the boundary.

III. MOLECULAR DYNAMICS SIMULATION APPROACHES
In this section and the next one, the various simulation approaches used to explore diffusion in an inhomogeneous medium are presented. We start by describing two well-known molecular dynamics (MD) approaches, namely, Langevin dynamics (LD) and Brownian dynamics (BD).
MD is a simulation approach based on solving equations of motions for particles comprising the system. In the most straightforward, but very computationally intensive approach, both the probe particle(s) whose diffusion is studied and the fluid particles are included explicitly and Newton's Second Law equations are solved numerically, where m i and x i are the mass and the position of the ith particle and F i is the sum of all forces acting on it. The computational expense can be reduced by considering the fluid particles implicitly, as is done in the LD and BD approaches.

A. Langevin Dynamics
In the LD approach ( [26], Sec. 14.4), only the equation of motion of the diffusing particle is considered explicitly and the action of the fluid particles on it is replaced by the friction force and a random force, and the resulting equation of motion is given in 1D by Eq. (1) with U = 0: Dividing by ζ and assuming that ζ is a function of x, we obtain or Unlike in the case of the first-order Eq. (10), there is no ambiguity due to different interpretations of the random term even when ζ (or D) is space-dependent, with all interpretations and all "correct" numerical methods giving the same solution in the limit when the time step δt → 0. This is because, unlike in the case of Eq. (10), the particle velocity does not diverge, thus when the equation is discretized, the displacement during a single step is much smaller [O(δt) instead of O(δt 1/2 )], so in the limit δt → 0 it no longer matters whether the position at the beginning or at the end of the step is used when calculating D [24]. In practice, some variation of the Verlet method ( [26], Sec. 13.4) is usually used to solve Eq. (52) numerically. By construction, Eq. (52) obeys the fluctuationdissipation theorem, thus in the case T = const the stationary distribution should coincide with the Boltzmann distribution ρ = const. Since this fact does not depend on the value of the parameters entering the equation, it should remain valid in the limit m → 0 when the inertial term vanishes, and thus this overdamped limit should correspond to the isothermal interpretation of Eq. (10). In fact, this was proved in Ref. [27]. However, in order to reproduce this result numerically, it is important to retain the inertial term decreasing its coefficient, rather than simply dropping it, and to make sure that the time step is always much smaller than the relaxation time, m/ζ. Remarkably, unlike Eq. (10) with α = 0, Eq. (52) can still be solved straightforwardly in the case of a sharp interface, as we will see in Sec. VI.

B. Brownian Dynamics
The motion of microscopic particles in a liquid is nearly always overdamped, with the time of relaxation to the terminal velocity, m/ζ, usually much smaller than any other relevant time scale. As mentioned above, the time step of the LD algorithm should be much smaller than even this extremely small time. Even if the particle mass m is increased artificially to increase the relaxation time and thus the maximum allowed time step, this time step may still end up being too small and as a result, the simulation will be inefficient. A seemingly straightforward approach is eliminating the left-hand side of Eq. (52) (the inertial term) altogether and solving the resulting firstorder equation. This equation coincides with Eq. (10), so the discussion of the ambiguity of the stochastic term in Sec. II applies. Normally in MD simulations, to solve Eq. (10) numerically the simplest Euler algorithm (often referred to as Euler-Maruyama in the stochastic case [28]) is used, in which, in order to obtain the particle position at time t + δt, the right-hand side of the equation is evaluated at time t. This approach, called Brownian dynamics, obviously corresponds to the Ito interpretation of the stochastic differential equation (10). Thus, surprisingly, LD (even in the overdamped limit) and BD will give different results when the friction coefficient ζ varies in space, as they correspond to different calculi. Just as LD, BD can be used in the case of a sharp interface without any complications, since these complications only arise for α = 0 when solving Eq. (10). Some difficulties arise when the random number r is chosen to be discrete (e.g., ±1), as described in Sec. IV C, so using r with a continuous distribution is preferable.

IV. MONTE CARLO APPROACHES
The main goal of this paper is to develop and test LMC algorithms for simulating a random walker at a viscosity interface corresponding to different calculi [α in Eq. (15)], including the popular Ito, Stratonovich, and isothermal calculi. There are several reasons to develop such algorithms. First, LMC algorithms are popular for studying diffusion problems [29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45] and in some of these problems studying the case when the effective viscosity is inhomogeneous in space (e.g., different degrees of crowding in different parts of the system) can be of interest (see, e.g., Ref. [30]). Second, LMC algorithms can offer advantages compared to MD (or, in general, solving numerically the Langevin equation) specifically in the case of a sharp interface. As we have seen, sharp interfaces are algorithmically problematic when solving the overdamped Langevin equation, except in the case α = 0 corresponding to BD. While the α = 1 (isothermal) case can be treated by LD, this may require a very small time step for the solution to be reliable. All other cases are in principle equivalent to Ito plus an additional drift force, but that drift force would be infinite at the interface and zero elsewhere, which would cause obvious numerical problems. LMC algorithms, on the other hand, are straightforward in all of these cases, as we will demonstrate.
Two different approaches will be taken. In the first method, changes in viscosity values are introduced by altering the probability p s that a particle does not jump during a MC jump attempt: the particle has a higher chance of "staying put" in the high viscosity region. In this approach, the jump lengths on the low and high viscosity sides are the same. Conversely, in the second method, the change in viscosity across the interface is included by altering the jump length (i.e., the mesh step of the lattice): The particle jumps shorter distances when the viscosity is higher. Here, the probability of not jumping is constant (except perhaps for sites adjacent to the interface). In both methods the time step is the same in both regions and stays constant throughout the simulation. This is essential for so-called numerically exact variants of LMC methods [21,36,40,42,44,46,47] and is also very helpful when several particles are simulated at once.
Both LMC methods are derived based on three requirements: (1) the diffusion should be unbiased away from the interface, (2) the diffusion rates should be correct in both regions, and (3) the probabilities of the moves at the interface should ensure that the stationary distribution given by Eq. (37) is reproduced. Only the last of these requirements depends on the calculus (the value of α), and so this value will only affect the rules of the algorithms at the interface, that is, in the only place where the viscosity gradient is present, as expected.
When deriving both algorithms, the following consideration is used. Suppose the particle moves on a lattice with a constant mesh step a and the time step is τ . At every successful step (i.e., ignoring the time steps where the particle stays put), the particle moves left or right with an equal probability. The mean-square displacement after N successful steps is Since in time ∆t there are ∆t/τ total steps and thus on average ∆t(1 − p s )/τ successful steps, and comparing this to Eq. (7), A. MC I: Variable ps, fixed jump length a

The basic algorithm
In the first Monte Carlo approach, all "successful" jumps have the same length and the dependence on viscosity is introduced by altering the probability of remaining still during a time step (Fig. 1).
FIG. 1: A schematic of the MC I algorithm. Different shades of gray depict regions of different viscosity, with the darker region being more viscous. Vertical lines correspond to the locations of 1D lattice sites. The distance between sites is the same in both regions. Arrows show transitions between sites: the thicker the arrow, the higher the transition probability. The notation for the probabilities corresponds to that used in the text. The top part illustrates the behavior of a particle away from the interface, which is independent of the calculus, while the middle and bottom parts are for a particle at the interface site in the two particular cases of the isothermal and Ito calculi, respectively.
Given the diffusion coefficient D and the time step τ , from Eq. (55) the probabilities to move right and left should be and the probability to stay put, Obviously, τ and a need to be chosen so that all of these probabilities are between zero and one for the diffusivities both to the left (D = D L ) and to the right (D = D R ) of the interface. Since D is inversely proportional to the fluid viscosity η [see Eq. (8)], it follows from Eq. (56) that p + and p − should be inversely proportional to η, or, introducing the viscosity η I 0 at which these probabilities are equal to 1/2 and correspondingly p s = 0, where sub-and superscripts L and R denote the left and right side of the viscosity interface, respectively, and "I" shows that these parameters correspond to the MC I algorithm. The same probabilities were used previously in Refs. [48,49]. The parameter η I 0 should not be larger than the smaller of the two viscosities η L and η R to ensure that both p L Is and p R Is are non-negative. In fact, even the situation where η I 0 is equal to or only slightly smaller than either η L or η R should be avoided, since some simulation artifacts, such as spurious oscillations both in time and in space, are possible when the waiting time is zero [50], or, more generally, the variance of the time between successful steps is small. On the other hand, from Eqs. (8) and (56), so this limitation on the allowed values of η I 0 restricts possible choices of τ I and a I . The fact that the probabilities of the moves are different in the two parts of the system is depicted schematically in the top part of Fig. 1.
Equations (58) and (59) should be used for lattice sites away from the interface. They ensure that the first two of the three requirements mentioned above (no bias and the correct diffusivities) are satisfied. The third condition, the correct ratio of the particle concentrations on the two sides of the interface, should be achieved by choosing a special set of probabilities for sites near the interface. It is convenient to choose the lattice so that the interface is on a lattice site. Let us denote this site 0, and the sites immediately to the left and to the right of it, −1 and +1, respectively. Since we have a single condition to satisfy for several probabilities, we can make an arbitrary choice that the probabilities of moves out of sites −1 and +1 are still given by Eqs. (58) and (59) and that only the two probabilities for the moves out of site 0 (plus the rest probability at site 0) are special. Based on Eq. (37), the condition that needs to be satisfied is that where n i is the average number of particles at site i in the stationary state. Since from detailed balance, and where p 0 I+ and p 0 I− are the probabilities of the moves from site 0 to the right and to the left, respectively, we get or This means that for η R = η L and α = 0 there is a bias at the interface (p 0 I− = p 0 I+ ). Equation (65) defines the ratio p 0 I− /p 0 I+ , but not each of these probabilities separately. In the isothermal case (α = 1), ρ = const in the stationary state and thus n 1 = n −1 ; it is reasonable to require that n 0 is the same as n −1 and n 1 , and this gives unambiguously The probability of staying put is then where the reduced viscosity is the harmonic mean of the two viscosities. Since it makes sense for the rate of moving away from the interface to always be related to some sort of average between the two viscosities, we can require, somewhat arbitrarily, that the probability of staying put be given by Eq. (68) for all values of α and use this as an additional condition together with Eq. (65) to obtain In fact, we have checked for the Ito calculus that the dynamics are insensitive to the choice of p 0 Is provided that p 0 I+ = p 0 I− , which in this case follows directly from Eq. (65). Results essentially indistinguishable from those generated by the above formulation were obtained in simulations with p 0 Is = 0 and p 0 I± = 0.5. Equations (70) and (71) have simple interpretations in two cases. In the isothermal case (α = 1), p 0 I− = p L I± and p 0 I+ = p R I± . In other words, the probability of the move to the left from the interface site (the move that is within the left region) is the same as the probability of any other move (from a non-interface site) within the left region, and the same is true for the right region. This property is illustrated in the middle part of Fig. 1. It ensures explicitly the absense of a net flux through any surface when the particle concentration is the same everywhere, the key property of the α = 1 case mentioned in Sec. II, since for every surface drawn between two sites there will be as many particles moving from left to right as from right to left. In the Ito case (α = 0), the probabilities of moves to the left and to the right from the interface site are equal (the bottom part of Fig. 1). This is consistent with the description of the Ito formulation in Secs. I and II as the case where any bias at the interface is absent.

Variants with all move attempts successful
When the ratio of the two viscosities, η L and η R , is very large or very small, then, since η I 0 should be chosen at least as small as the smaller of the two, the probability of moving in the more viscous region will be very low, according to Eq. (58). This makes the algorithm very inefficient in this case, since many random numbers need to be generated per single successful move. However, one may note that the number of steps until the next successful move is exponentially distributed, with the probability that this number is n being where p Is is either p L Is , p R Is , or p 0 Is , depending on the location of the particle. Instead of waiting for many steps until a successful move, one can move at every step, but advance the clock by nτ I instead of τ I , generating n from the distribution (72). The resulting procedure is exactly equivalent to the original one. Strictly speaking, in this procedure the time step is variable, which in some cases is undesirable, as mentioned above; however, it is still always a multiple of the basic step τ I . It is also possible to always advance the clock by the same amount t Iw , equal to the average waiting time between successful steps, where η is either η L , η R , or η * . While not strictly equivalent to the algorithm given by Eqs. (58) and (59), such an algorithm should ultimately give the same results in the limit a, τ → 0. The advantage is that the need to draw random numbers from the distribution (72) is avoided. However, the time step will no longer be a multiple of τ . A potentially more serious problem is simulation artifacts that may appear when the waiting time is constant, as mentioned above. Regardless of the utility of this algorithm for actual simulations, it is useful in the analysis of the diffusion problem, as we discuss below.
In both variants of the algorithm (with random and deterministic clock advancement), at each step the moves to the left and to the right are equiprobable, except when the particle is right at the viscosity interface the probabilities are These probabilities are such that their ratio p ′0 I+ /p ′0 I− is the same as p 0 I+ /p 0 I− , but their sum is unity, so the rest probability is zero.
Note that in the algorithm presented here the particle moves to a neighboring site at every step, that is, the particle motion is simply an ordinary random walk on a 1D lattice, and it is unbiased everywhere, except perhaps at the interface. (We ignore for the moment the fact that each step takes a different amount of time in different regions.) In fact, in the special case of the Ito calculus the random walk is unbiased at the interface as well, since in that case Eqs. (74) and (75) give p ′0 I− = p ′0 I+ . In general, the trajectory of the particle can be thought of as consisting of segments of an unbiased random walk terminating at the interface; every time the interface is reached, the particle chooses to continue on the left side with the probability p ′0 I− or on the right side with the probability p ′0 I+ . Each such choice is independent of all of the previous choices and of whether the particle has reached the interface from the left or from the right. We will use this consideration in Sec. VIB and in Appendix B to derive some properties of the system.

B. MC II: Fixed ps, variable jump length
In the second MC approach, the probability of remaining still, p IIs , is independent of the viscosity and instead the length of a successful jump is dependent on the viscosity. To keep the particle on-lattice, we introduce a mesh with a different step in the two regions (Fig. 2). From Eqs. (55) and (8), when p IIs and τ II are fixed, the mesh step a ∝ √ D ∝ 1/ √ η, so the mesh steps on the left and right sides are, respectively, where a II0 and η II 0 are constants and While it is possible to put the interface at one of the lattice sites, as in MC I, we now place it between two sites to illustrate the possibility of this choice. The distance between the two sites closest to the interface (one on each side) is chosen to be the average of the mesh steps on the two sides, a * II = (a L II + a R II )/2, although the exact value does not matter in the limit a L II , a R II → 0.

FIG. 2:
A schematic of the MC II algorithm. Different shades of gray depict regions of different viscosity, with the darker region being more viscous. Vertical lines correspond to the locations of 1D lattice sites. The distance between sites is smaller in the region of higher viscosity. Arrows show transitions between sites: the thicker the arrow, the higher the transition probability. Notation for the probabilities corresponds to that used in the text. The top part illustrates the behavior of a particle away from the interface, which is independent of the calculus, while the rest of the figure is for particles on sites adjacent to the interface on both sides in the cases of the isothermal, Ito, and Stratonovich calculi.
For all sites not adjacent to the interface, on both sides, the probabilities of the moves are the same, with and the probability of remaining still, p IIs , is the same at all these sites, as well. The choice of p IIs is in principle arbitrary. While p IIs = 0 seems the most efficient choice at first glance, since the particle would move at every step of the algorithm, this, as mentioned, can introduce artifacts and, in fact, we have recently shown [50] that p IIs = 2/3 is optimal from the point of view of accuracy of the algorithm and may even be more efficient than p IIs = 0, since a coarser mesh can be used to achieve the same accuracy. For this reason, we choose p IIs = 2/3 in this work. As in MC I, the probabilities of the moves need to be modified for sites adjacent to the interface. A reasonable choice is to allow the probabilities of the two moves crossing the interface to change but keep the probabilities of the two other moves leaving these sites that do not cross the interface the same as given by Eq. (79). We use the following notation for the probabilities of the moves that are allowed to change: p 0 II+ denotes the probability of the move to the right from the site left of the boundary and p 0 II− is the probability of the move to the left from the site right of the boundary. This notation is similar to that we used for MC I, but it has a somewhat different meaning: in MC I these were the probabilities of the moves from the same interface site, but here these moves cross the interface starting from different sites.
The detailed balance condition is where n L (n R ) is the average number of particles per site to the left (right) of the interface in the stationary state.
To relate this condition to the stationary distribution of Eq. (37), we need to keep in mind that the density of sites (their number per unit length) is different in the two regions. The particle concentration ρ is proportional to both the number of particles per site and the site density; the latter is inversely proportional to the mesh step. Therefore, using Eqs. (76) and (77) for the mesh steps, or p 0 Again, as in MC I [cf. Eq. (65)], we have one equation for two probabilities. As an additional condition, we choose which, together with Eq. (82), gives This is a simple choice that guarantees that both p 0 II+ and p 0 II− are always between zero and one for any value of α and any choice of p IIs . Note that always either p 0 II+ = p L,R II± or p 0 II− = p L,R II± , so only one of these probabilities (at most) gets modified compared to the bulk value. The corresponding probability of staying put is changed accordingly. In the special case α = 1/2 (Stratonovich calculus), p 0 II+ = p 0 II− = p L,R II± , so all probabilities of all moves are the same (Fig. 2). While this means that for α = 1/2 there is no flux across the interface when n R = n L , this is not the same as ρ R = ρ L (a uniform particle distribution).
We note that a similar, but off-lattice approach was successfully used recently by one of the authors [51] to treat a problem corresponding to the α = 1 case.

C. Comparison between MC II and BD
It is instructive to compare the MC II algorithm described here in the Ito case (α = 0) to the BD algorithm, a single step of which is given by Eq. (15) with α = 0, specifically its variant with the random number r = ±1. Similarities are obvious: in both cases the time step is fixed and the jump size depends on what region the particle is in, being proportional to √ D (or inversely proportional to √ η). However, what happens at the interface is different for these two algorithms. In MC II, in order to stay on-lattice, the lengths of the jumps across the interface are chosen to be the same in both directions, but in BD the jump length depends on what region the jump originates in and thus is different for the left-toright jumps and the right-to-left jumps. This difference is crucial and makes it possible to achieve the correct particle concentration jump across the interface in the stationary state while keeping the probabilities of all jumps the same (effectively equal to 1/2, since there is no bias in BD and no staying at rest), whereas in MC II obtaining this concentration jump requires setting p 0 + = p 0 − . This is actually a rather subtle issue. Indeed, in the stationary state, the left-to-right and the right-to-left fluxes across the interface should be equal. But then one could argue for the BD algorithm that since the flux is equal to the product of the particle concentration and the jump size and the latter is proportional to √ D, then the stationary concentration ratio ρ R /ρ L should be equal to D L /D R , instead of the required D L /D R . The reason this argument is wrong is that instead of a perfectly sharp jump in the concentration, it actually changes from ρ L to ρ R over a length scale on the order of the length of the particle jump. Therefore the particle concentrations at the points where the jumps across the interface originate are not equal to the values far from the interface (ρ L and ρ R , respectively). Note that, while this is true even when r is continuously distributed, in the case r = ±1 there is an additional complication that affects the practical use of the BD algorithm: the particle concentration keeps varying wildly even far away from the interface, with the typical length scale of the oscillations again on the order of the particle jump size. As a consequence, the results depend strongly on how the data are binned. For this reason, using r from a discrete set is best avoided in practical applications of the BD algorithm.

V. SIMULATION PROTOCOL
The system is defined such that the left side extends from a wall at x = −d/2 to the interface at x = 0 and the right side extends from the interface to a wall at d/2, where d is the total length of the system. Thus the interface is placed equidistant between the walls. For all simulations described in this paper, the particle is initially placed at the interface. In the case of the MC II algorithm, when there is no site exactly at the interface, it starts at one of the sites immediately to the left and to the right of the interface with equal probabilities. Two types of simulations are performed. First, simulations with reflecting walls at x = −d/2 and x = d/2 are conducted. To obtain the stationary particle distributions, relatively long simulations must be performed and thus we focus on a single case where η R /η L = 4. Additional simulations with absorbing walls are performed to investigate whether particles preferentially end up at one wall (herein referred to as an ultimate preferential direction) and to measure the conditional mean first passage time (MFPT) to each wall. In most of these simulations, the viscosity on the left-hand side, η L , is held fixed and the viscosity on the right-hand side is varied from η R = 0.1η L to η R = 10η L ; we also use η R /η L = 50 or even, in Appendix A, η R /η L = 1000.
For static quantities, such as the stationary particle distributions and the ultimate preferential direction, only the ratio of the viscosities across the interface influences the results -the actual values of the viscosities and the proportionality factor κ in Eq.  (71)], if η L and η R are both changed by some factor and at the same time the adjustable constant η I 0 is changed by the same factor, then the probabilities of the moves remain the same. What does change is the time scale, since changing η I 0 changes the time step τ I , according to Eq. (60), which corresponds to speeding up or slowing down the dynamics. Since all approaches are expected to also give dynamically correct results such that the MFPTs and the decay of the transient behavior towards the stationary state should be consistent for any particular α value (with BD corresponding to α = 0 and LD to α = 1), a comparison of the dynamical quantities between the approaches is also meaningful. To make such comparisons, it is obviously necessary to know how the time steps of different algorithms are related.
It is convenient to choose the time step of the MC I algorithm as the unit of time and express all other times and time steps in terms of this unit. That is, numerically, the time step of MC I will be chosen to be Likewise, the mesh step of MC I is chosen as the unit of length (a I = 1). In these units, we According to Eq. (60), in these units In MC II, we choose arbitrarily a II0 = √ 2a I = √ 2 in Eqs. (76) and (77). The lattice is placed in such a way that the interface is at distance a L II /2 from the site adjacent to it on the left and at distance a R II /2 from the site on the right. For the leftmost site, −d/2 ≤ x < −d/2 + a L II , and for the rightmost site, d/2 − a R II < x ≤ d/2. It is convenient to express the viscosity in units of η II 0 . However, for direct comparison with MC I the viscosity unit should be the same in both cases; thus we need η II 0 = η I 0 (= 1). Using Eqs. (60) and (78), so or, using the chosen values for a II0 and a I and keeping in mind that we have chosen the rest probability in MC II p IIs = 2/3 (as explained in Sec. IV B), While choosing the values of the viscosities, we need to keep in mind that in MC I the viscosities cannot be smaller than η I 0 = 1. For the reflecting wall simulations, η L is set to 2, since, as we have mentioned, using a value too close to unity can introduce artifacts in the distribution function. Correspondingly, η R = 8 for η R /η L = 4. For the absorbing wall case, for the simulations where we vary η R /η L from 0.1 to 10, η L is set to 10 so that η R can be varied from 1 to 100. Using η R = 1 is not a problem in this case, since we are not interested in the details of the particle distribution or the FPT distribution, but only the MFPTs, in which case all oscillations and other artifacts average out. For η R /η L = 50 and 1000, we use η L = 2, as in the reflecting case, so η R is respectively 100 and 2000.
For both MC algorithms, when the particle attempts to move left from the leftmost site or right from the rightmost site, in the case of absorbing boundaries it disappears and the simulation is stopped, while in the case of reflecting boundaries the move is rejected.
As for the MD simulations (BD and LD), we likewise use d = 50 and treat the particle as pointlike in the sense that it can approach the walls infinitely closely. The random terms in the equations of motion are normally distributed. When attempting to cross a wall, the particle disappears in the absorbing case and is "reflected" in the reflecting case. The latter means that, for example, when the particle crosses the right wall at d/2 and would end up at x = d/2+∆, it is placed at x = d/2−∆ instead and in the LD case the sign of the particle velocity is changed; similarly for crossing the left wall. The equations of motion that are solved numerically [Eqs. (10) and (52), respectively] contain the diffusivity D. For comparison to MC, it needs to be expressed in terms of the viscosity given in the same units as in MC simulations, that is, in units of η I 0 . Using Eqs. (8) and (88), in these units For the BD time step δt, it is reasonable to choose a value similar to the time resolution of the MC schemes.
We have chosen δt = 2 for simulations with reflecting boundaries and δt = 20 for simulations with absorbing boundaries. While the latter may appear too large, we have checked that using a ten times smaller step produces equivalent results. Ultimately, this is not surprising: What matters is that the time step is much smaller than the typical time that it takes the particle to reach the boundary starting at the interface, and this condition is satisfied. The situation with LD is somewhat more complicated. First of all, the initial velocity of the particle needs to be specified; it is chosen equal to either + kT /m or − kT /m at random. Second, there is one other relevant combination of parameters, m/kT , and the associated time scale, mD/kT = m/2ηkT . It determines the inertia in the system, in particular, the correlation time of the particle velocity (or the time of relaxation to the terminal velocity). Since the intent is to model an overdamped system in which inertia is negligible, this time scale should be small, in any case smaller than the time it takes the particle to reach the wall from the interface. We chose m/kT = 4 for reflecting boundaries and m/kT = 400 for absorbing boundaries. In the least favorable case, η = 1 for absorbing boundaries, the relaxation time is t r = m/2ηkT = 200. The typical particle displacement over this time is √ 2Dt r = t r /η ≈ 14 which is less than a factor of 2 smaller than d/2 = 25. While some discrepancies between LD and the MC and analytic results are indeed found at the lowest viscosities studied, good agreement is obtained for the majority of η R values. The time step in LD should be much smaller than the relaxation time; we choose δt = 0.02 for reflecting boundaries and δt = 0.2 for absorbing boundaries. The least favorable case is the largest η in our series for absorbing boundaries, η = 100 (in the separate runs with η R =2000 we do not use LD, only MC I). The relaxation time in this case is t r = 2, so this condition is satisfied. Note that the time step for LD is much smaller than for BD and in fact, LD is the most time-consuming algorithm of all those considered in this paper. This is expected, since in LD the time step should be much smaller than the relaxation time, which should itself be small; we note, however, that the efficiency of the LD algorithm could be improved without sacrificing its accuracy by using different m/kT and δt for different viscosity ratios.

A. Reflecting Walls: Distributions
In this section, we look at the particle concentration distributions for the case when the particles start at the interface. We first present the distributions for the Ito and isothermal calculi using the MC I and MC II algorithms, comparing them with the results obtained using BD (for Ito) and LD (for isothermal). We then compare the distributions for several different values of α using one of the MC algorithms (MC I). In all cases, we show the distributions at a time shortly after the start of the simulation and at a long time corresponding to the stationary state. In addition, for the Ito and isothermal calculi we also include one intermediate time, such that the distribution on the low-viscosity side has already reached the quasistationary state, but that on the high-viscosity side still resembles the short-time distribution. For each curve, the total number of runs is 100000. In each run, a single particle is simulated and its position is determined at the times of interest. In MD simulations, the positions are binned with a unit bin size and the centers of the bins at half-integer x. The plotted values ρ i are proportional to the numbers of particles n i found at the given time at the site (or in the bin) i and inversely proportional to the distance between the sites (width of the bin). The distributions are normalized so the area under the curve is unity (thus, in effect, these are probability density distributions). In particular, for MC I, BD and LD, since the distance between the sites (or the width of the bins) is unity, the normalization condition is where the sum runs over all sites (bins), and thus For MC II, where the first sum is over the sites in the left region and the second sum is over the sites in the right region, so We also plot the ratio of the numbers of particles on the right and left sides as a function of time. When calculating this ratio, particles exactly at the interface (which is possible in the MC I case) are ignored.

Ito calculus
The probability distributions ρ as functions of position at three different times are shown in Fig. 3(a)-(c), for all algorithms yielding Ito calculus (MC I, MC II and BD). In Fig. 3(a), the distribution shortly after starting the simulations, at time t = 100, is shown. By this time, only very few particles reach the walls, so the distribution is essentially the same as it would be in infinite space, without the walls. In this case, theoretically the distribution is expected to follow (97) This satisfies the diffusion equation (45) in both regions, as well as the jump condition (44) and the flux matching condition (46), the latter because ∂ρ/∂x is zero both at x = −0 and x = +0. Thus the parts of the distribution on the left side (L) and the right side (R) should both have a Gaussian form, but there should be a significant discontinuity in the distribution with a higher concentration of particles on the high viscosity side. These features are indeed reproduced in Fig. 3(a). The discontinuity develops instantaneously (in a simulation, within a few steps from the beginning) and results from the trapping effect of high viscosity regions in the Ito formulation. Figure 3(c) shows the distribution after a long time (t = 30000) thus corresponding to the stationary state. While the transient piecewise-Gaussian shape has disappeared, the discontinuity persists yielding a stepfunction form. Recall the Ito stationary state condition [Eq. (39)] which, using D ∼ 1/η, can be written as ρ(x)/η(x) = const. Hence, as expected, the probability density on L, where the viscosity is low, is less than on R, where the viscosity is high, and the ratio is Given the normalization condition (ρ L + ρ R )d/2 = 1, For the case presented in Fig. 3 where d = 50, η L = 2, η R = 8, we have ρ L = 0.008, ρ R = 0.032. These values are indicated on Fig. 3(c) where ρ R /ρ L = 4 as expected for η R /η L = 4. Note that at short times [as in Fig. 3(a)], the ratio of the probability densities immediately to the right of the interface and immediately to the left of the interface is also the same. In fact, this is true at all times, as follows from Eq. (44). The sharp drops in the distribution near the walls in Fig. 3(c) are discretization and binning artifacts. In Fig. 3(d), the ratio of the number of particles on R to L, denoted N R /N L , is plotted as a function of time. This ratio is equal to the ratio of the areas under the parts of the distribution curves on R and L. Note that N R /N L increases from a value around 2 and saturates at a value around 4. Hence, there are always more particles on the high viscosity side: Under Ito conditions, the high viscosity side acts as a trap. While the saturation value agrees with the stationary state condition of Ito calculus, it is interesting that, beyond extremely short time behavior which depends more on the simulation algorithm than the physics, the initial point for the increase during the transient portion is 2. Indeed, from Eq. (97) valid at short times, the ratio should be η R /η L = √ 4 = 2. This is the ratio of the heights of the half-Gaussians (which is equal to D L /D R = η R /η L ) multiplied by the ratio of their widths (which is equal to D R /D L = η L /η R ). At longer times, the particles begin to find the edge of the system and bounce back. Correspondingly, the system begins to move towards equilibrium where four times as many particles are found on the high viscosity side. The transition starts when the particles on the lower-viscosity side reach the boundary, which happens around time t 1 such that 2D max t 1 = (d/2) 2 , thus t 1 = d 2 /(8D max ), or, using Eq. (92), where D max is the higher of the two values of the diffusivity and η min is the corresponding (lower of the two) value of the viscosity. The transition ends when the particles on the higher-viscosity site also reach the boundary, which occurs around For our values of the parameters (d = 50, η min = η L = 2, η max = η R = 8), t 1 = 1250 and t 2 =5000. These two times are marked on Fig. 3(d). The time at which the distributions in Fig. 3(a) are obtained is much smaller than t 1 , thus corresponding to the situation before the transition starts, and the value of the ratio N R /N L at that time is indeed close to two, according to Fig. 3(d).
On the other hand, the time at which the distributions in Fig. 3(c) are obtained is much larger than t 2 and thus corresponds to the stationary situation after the transition has ended; the value of N R /N L is close to four. In Fig. 3(b), we show the distribution at time t = 1500, which is between t 1 and t 2 . By time t such that particles on the low-viscosity side would have traversed that side many times, so the distribution on that side should be nearly flat resembling the stationary distribution, but on the other hand, very few particles would have reached the wall on the high-viscosity side, so the distribution on that side should still be close to the short-time Gaussian distribution. Even though in our case the ratio t 2 /t 1 = η R /η L is only four, so, strictly speaking, satisfying Eq. (103) is impossible, Fig. 3(b) can still serve as a reasonable illustration of that regime. Note that the distribution on the low-viscosity side is not truly stationary, but rather quasi-stationary, since there is a net outflow of particles to the high-viscosity side, as indicated by the gradual increase of N R /N L between t 1 and t 2 [see Fig. 3(d)]. This outflow is slow enough for the distribution to remain nearly flat, but the height of the distribution gradually decreases. In essence, the height of the plateau on the low-viscosity side adiabatically follows the height of the peak on the high-viscosity side (which decreases as this peak broadens), so that the boundary condition (44) is satisfied at all times. On the other hand, this outflow is large enough that the slope of the distribution function is significantly nonzero immediately to the right of the interface in Fig. 3(c). However, the viscosity ratio is rather small here and also t is close enough to t 2 that particles have already started to reach the right wall (note a nonzero particle concentration next to it). As we argue in Appendix A, for a very high viscosity ratio and when the condition (103) is truly satisfied, the derivative at x = +0 is negligible and the distribution on the high-viscosity side is a Gaussian. Taking into account the general character of the distribution at intermediate times between t 1 and t 2 , we can also obtain approximately the time dependence of the ratio N R /N L in this time interval. The number of particles on R, N R , is roughly proportional to ρ(+0, t) times the width of the distribution on that side, which is ∝ √ t. The number of particles on L, N L , is proportional to ρ(−0, t) times the length of that part of the system, d/2, which is time-independent. Thus √ t, and, given that, according to Eq. (44), the ratio in the square brackets is timeindependent, This agreement between the theory at α = 0 and the Brownian dynamics results confirms that the BD simulations correspond to Ito calculus. Further, note the overall good agreement between the BD simulations and both MC approaches for all the results shown in Fig. 3. Considering factors such as discretization, binning, and noise, the fact that the differences are small between the data sets is a convincing validation of the two MC approaches we have developed to simulate Ito calculus.

Isothermal calculus
The equivalent plots for all algorithms yielding the isothermal calculus (MC I, MC II and LD) are shown in Fig. 4. Note that the distributions are now continuous across the interface, in agreement with Eq. (44) for α = 1. For the short time behavior shown in Fig. 4(a), where, again, t = 100, the Gaussian-like distribution on the left side is wider than that on the right side, as in the Ito case. Theoretically, the distribution at short times is given by which is the same as Eq. (97) for the Ito case, except for the prefactors that ensure the continuity across the interface. At long times, in the stationary state, there is an equal probability to find a particle anywhere in the system, as is seen in Fig. 4(c) (t = 30000). This is consistent with the isothermal stationary state condition [Eq. (42)]. The stationary value of the density is obviously given by The time dependence of the ratio of the number of particles on R over L is shown in Fig. 4(d). Here, the transient behavior begins from the ratio of approximately 0.5. That is, given the bias at the interface to end up on the lower viscosity side, very rapidly there is a greater number of particles on L, and, contrary to the Ito results where there are twice as many particles on the high viscosity side near the beginning, for the isothermal results there are half as many particles on the high viscosity side. This ratio can again be obtained by considering the peak height and width. Here, the peak height is the same on both sides while the ratio of the widths is again η L /η R . Correspondingly, the ratio is N R /N L = √ η L / √ η R = 1/2. As the system moves towards equilibrium, N R /N L again climbs and plateaus, but now it saturates at 1.0 in accord with the isothermal stationary state condition. The transition between the initial and the final values again starts around t 1 given by Eq. (101) and ends around t 2 given by Eq. (102). At times between t 1 and t 2 , the distribution is nearly flat (quasi-stationary) on the low-viscosity side and Gaussian-like on the high-viscosity side, again, as in the Ito case [see Fig. 4(b) for t = 1500]. However, unlike in the Ito case, we show in Appendix A that even for a large viscosity ratio and t 1 ≪ t ≪ t 2 , the slope at +0 is significantly nonzero, and, in fact, the interface lies at the inflection point of the Gaussian-like curve rather than at its maximum. Despite this, the crude arguments leading to Eq. (104) still apply, so we expect that time dependence to still be approximately valid here. We compare the time dependences of N R /N L for different calculi in the next section.
In all four plots, again, there is good agreement between all three algorithms.

Comparison between calculi
In Fig. 5, we compare the results for the particle distributions and the ratio N R /N L for several different calculi (values of α), including the already considered Ito (α = 0) and isothermal (α = 1) limit cases, using the MC I algorithm. Essentially identical results are obtained using MC II (not shown). Figure 5(a) shows the distributions at t = 100 and Fig. 5(b) presents the distributions at t = 30000. At short times [ Fig. 5(a) is always piecewise-Gaussian with (107) There is a jump at the interface in all cases, except α = 1, in agreement with Eq. (44). At long times [ Fig. 5(b)], the distribution is piecewise-uniform, with which follows from Eq. (37) for the stationary distribution. Again, there is a jump for α = 1. In Fig. 5(c), the ratio N R /N L is plotted. The short-time value of this ratio is which is the product of the ratio of the peak widths [(η R /η L ) −1/2 ) and the ratio of the peak heights For 0 ≤ α < 1/2, N R /N L is always larger than 1, so there are always more particles on the higher-viscosity R side than on L. For α = 1/2, the ratio is one initially, but then grows above unity. For α = 1, there are always more particles on the lower-viscosity site, but the ratio approaches unity from below as t → ∞. Finally, for 1/2 < α < 1, there are more particles on the low-viscosity side initially, but the ratio crosses unity at a finite time and eventually there are more particles on the higher-viscosity side. The transition between the short-time and the long-time values of N R /N L is roughly between t 1 and t 2 given by Eqs. (101) and (102), and the time dependence should still be roughly given by Eq. (104). Note that at both small and large times, if the particle distribution on L is multiplied by (η R /η L ) 1−α , then the resulting function is the same for all α: This rescaling eliminates the jump at zero, and the shapes of the distributions in a given region and at any given time are calculus-independent, both at short times (when they are Gaussian with the calculus-independent width) and at long times (when they are flat). A reasonable question then is if this is also true at intermediate times and thus at all times. The answer, perhaps somewhat surprisingly, is no. This can be seen by comparing the shapes of the intermediate-time distributions for the Ito [ Fig. 3(b)] and isothermal [ Fig. 4(b)] cases that are indeed somewhat different, and, as mentioned, the difference is significant in the high-viscosity-ratio limit (see Appendix A). This difference is also easy to explain: If we take, for example, the solution for the Ito case and multiply its left half by η R /η L , the resulting function will satisfy the diffusion equation (45) and the jump condition (44) for α = 1, but not the flux continuity condition (46), since only the lefthand side of the latter will be multiplied by η R /η L , but not the right-hand side. The cases t → 0 and t → ∞ are exceptions, because in these limits ∂ρ/∂x is zero at both x = −0 and x = +0. However, for intermediate times this cannot be the case, since there must be some net flux between the two sides, or otherwise N R /N L would remain constant. As a consequence, while the curves in Fig. 5(c) for different α do look similar, and, in particular, the transition between the low-time and the high-time values of N R /N L occurs over roughly the same time interval and according to approximately the same ∼ t 1/2 dependence in all cases, these curves are not perfect rescaled copies of each other: Multiplying these functions by (η R /η L ) α so they match at t → 0 and t → ∞ does not make them coincide everywhere, as illustrated in Fig. 5(d), where small differences between the rescaled curves around t 2 are apparent. That plot also compares the curves with the t 1/2 dependence predicted for large viscosity ratios; the agreement is reasonable given that in our case the ratio η R /η L = 4 is not large.

B. Absorbing Walls: Preferential Direction
We now describe the results of simulations with absorbing walls, that is, each single-particle simulation run is stopped once the particle reaches one of the two walls. For the Ito and isothermal calculi, we first look at whether the particles prefer to end up on one of the walls (i.e., whether there is an ultimate preferential direction) and find the conditional MFPTs, separately for particles reaching each wall. The ultimate preference is then compared to the preference of the particle to spend more time on one or the other side, which is analyzed in two different ways. Unlike the case of reflecting boundaries, where a single set of viscosities was considered, here we vary the ratio of the viscosities, with the viscosity on L kept constant. For each set of viscosities, we carry out 10000 single-particle runs, again starting at the interface, as described in Sec. V. For the Ito calculus we also obtain the full conditional FPT distributions for a single representative viscosity ratio using 10 6 runs.

Ito calculus
Results for the absorbing wall simulations in the case of Ito calculus are shown in Fig. 6. Figure 6(a) displays the fraction of events absorbed on R, π R (dashed lines) and L, π L (solid lines) for all three algorithms. For each algorithm at any viscosity ratio, the number absorbed at L equals the number absorbed at R within error. Thus, under Ito conditions, there is no ultimate preferential direction in the system. This result may seem paradoxical at first, given that for large η R /η L there are much fewer particles on L than on R at any given moment of time, which we have seen for reflecting boundaries [ Fig. 3(d)] and which turns out to be true for absorbing boundaries as well. The primary reason for the equality in Fig. 6 (a) is that particles diffuse faster on L, so even a smaller concentration creates a significant flux towards the left wall. The conditional MFPTs to each wall (which we denote τ L and τ R ) are plotted in Fig. 6(b). As would be expected, the time to the lower viscosity side is always less than the time to the high viscosity side with an intersection corresponding to equal times at η R /η L = 1. A less obvious fact is that the two times are always of the same order of magnitude even when the viscosities are very different. In fact, both dependences are linear as demonstrated by straight-line fits (black lines) to the averages of all three algorithms, and the slope for the R data is, within the uncertainty, two (2.04 ± 0.07) times that of the slope for the L data, so the ratio of the two times approaches two as η R /η L → ∞. This may seem surprising, given that the time t 1 [Eq. (101)] around which a significant (not exponentially small) number of particles first reaches the left wall is much smaller in this limit than the corresponding time t 2 [Eq. (102)] for the right wall. The explanation is that, while it is indeed true that the particle flow to the wall on the low-viscosity side starts at time ∼ t 1 , it continues until time ∼ t 2 , as the particles keep "leaking" from the high-viscosity side until those that have not been absorbed by the wall on the low-viscosity side are absorbed by the other wall. Even if the means are of the same order of magnitude, the conditional first-passage time distributions P L,R (t) are very different when the viscosities are very different. For the wall on the low-viscosity side, the flux reaches its maximum very fast, by time ∼ t 1 and then decreases very gradually (as a power law) until time ∼ t 2 , after which it decays exponentially. But for the wall on the highviscosity side, the flux is essentially zero until time ∼ t 2 and then rises and decays relatively quickly. This behavior of the distributions is illustrated in Fig. 7 obtained for a large viscosity ratio (η L = 2, η R = 100, so η R /η L = 50) using the MC I method with the waiting time selected from a distribution as given by Eq. (72). The best fit to the power-law region has the exponent ≈ −0.8; we show in Appendix A that in the limit of a very large viscosity ratio this exponent is −0.5. To derive the properties revealed in Fig. 6, consider the variant of the MC I algorithm with all steps successful (no staying put) and the time increment at each step t w given by Eq. (73) (see Sec. IV A 2). Recall that in the Ito formulation, there is no bias at the interface [from Eqs. (74) and (75), p ′0 I− = p ′0 I+ = 1/2], or anywhere else in the system. Hence, the trajectory of a particle is just an unbiased random walk, except for the fact that the time increment at each step changes depending on what side of the interface the particle is at. This exception, of course, does not change the fact that this random walk is equally likely to first reach the left or the right wall, since it starts at equal distances from them, which explains the result in Fig. 6(a). To get insight into the results for the MFPTs in Fig. 6(b), we use the fact that if a particle is released equidistant between two absorbing walls that are at distance d from each other and this particle does a random walk with step size unity eventually reaching the right-hand wall, then on average it spends steps on the right-hand side and steps on the left-hand side (see Appendix B). Then and similarly, where t L Iw and t R Iw are the time steps [given by Eq. (73)] on L and R, respectively, or, using the conventions (86) and (87), For η L = const, the dependence of both times on η R /η L is indeed linear, and the ratio of the slopes is two, in agreement with the data in Fig. 6(b). In our case (d = 50, η L = 10), This should be compared with the fits to the data (black lines in Fig. 6), which are τ fit R = (2.35 ± 0.03) × 10 3 +(4.26 ± 0.09) × 10 3 (η R /η L ), (120) τ fit L = (4.58 ± 0.04) × 10 3 +(2.09 ± 0.06) × 10 3 (η R /η L ).
The slopes are indeed in agreement; the small disagreement in the intercepts (about 10%) is likely due to discretization effects. In the limit η R /η L → ∞, Eqs. (116) and (117) become where t 2 is given by Eq. (102). This shows that in this limit both times become independent of the lower viscosity and in particular, both are ∼ t 2 . We discuss the relation between this result and the FPT distribution in Appendix A.

Isothermal calculus
The same data for the case of isothermal calculus are shown in Fig. 8. There are significant quantitative and qualitative differences compared to the case of Ito calculus. First, examining Fig. 8(a), there is a clear ultimate preferential direction in the system. When η R /η L = 0.1, approximately 90% of the particles are absorbed at the R wall. Conversely, when η R /η L = 10, 90% of the particles are absorbed at the L wall. Plotting the data with a log scale for the viscosity ratio [inset to Fig. 8(a)] demonstrates that this symmetry holds for all simulated viscosity ratios: where π L,R (z) are the probabilities of absorption at the left (right) wall when the ratio of viscosities on R and on L is z. This symmetry is a necessary outcome considering that in the algorithms, it is only the ratio of viscosities, and not the absolute value, which influences the dynamics (beyond simply making all events take longer). The ultimate preferential direction under isothermal conditions arises from the bias at the interface favoring the low viscosity side. Since this bias increases with a larger disparity in viscosities, the strength of the preference in direction also grows. It must, however, saturate at 100% resulting in the behavior shown in Fig. 8(a). The MFPTs, shown in Fig. 8(b), are also quite different from those corresponding to the Ito condition. Again, as expected, the time to the lower viscosity side is lower than the time to the high viscosity side with an intersection at η R /η L = 1. However, the dependence across the entire viscosity range is clearly not linear. In fact, considering data for which η R /η L > 1, while τ R does increase linearly with increasing η R /η L (albeit only when this ratio is large), τ L appears to saturate.
To understand the results in Fig. 8, we turn again to the variant of the MC I algorithm with all steps successful. For generality, we consider an arbitrary α and then specialize to α = 1. For α = 0, the segments of the particle trajectory between visits to the interface site are still random walks, but when the particle is at the interface, it chooses to continue on L or on R with the probabilities given by Eqs. (74) and (75), respectively. Each such trajectory segment is contained entirely either within the left half or within the right half and for brevity, we refer to these groups of segments as left and right segments. Since the probabilities to reach the respective wall during the left and right segments are equal, the total probability to reach the left wall, π L , is equal to the fraction of left segments, which in its turn is equal to the probability to start a left segment while at the interface, which is p ′0 I− . Correspondingly, the probability to reach the right wall, π R , is p ′0 I+ . So Note that π L → 1 for η R /η L → ∞ whenever α > 0, whereas for α = 0, π L ≡ 1/2, as we have seen in the previous section. In the isothermal case (α = 1) These functions are plotted in Fig. 8(a), and agreement with the simulation data is excellent. Equations (127) and (128) satisfy the symmetry condition (124). At η R /η L = 0.1, π R = 10/11 ≈ 91%, while at η R /η L = 10, π L = 10/11 ≈ 91%, again, in agreement with the simulations.
To calculate the MFPTs, we use the result obtained in Appendix B that for a random walk that is unbiased, except at the interface (where it chooses to go left with probability p ′0 I− or right with probability p ′0 I+ ), the mean number of steps that the walk absorbed by the right wall spends on R is while on L it spends on average steps. Then the MFPT for the walk ending on the right wall is where we have used Eq. (73) for t L Iw and t R Iw and conventions (86) and (87). For walks ending on the left wall, analogously there are on average steps on L and steps on R, and then In particular, in the isothermal case (α = 1), Equations (135) and (136) with η L = 10 and d = 50 are plotted in Fig. 8(b). The agreement with the simulations is very good; discrepancies can be attributed to discretization effects as well as inertial effects in the LD case. Equations (116) and (117)  for α = 0, for 0 < α < 1, For η L = const, this always gives a linear dependence on (η R /η L ) for τ R , regardless of α, and, unsurprisingly, τ R ∼ t 2 , again for any α. For τ L , the dependence is ∝ (η R /η L ) 1−α , which is linear for α = 0, reproducing the result of the previous section, sublinear for 0 < α < 1, and a constant for α = 1, in agreement with Fig. 8(b). Note also that τ L is now of order t α 1 t 1−α 2 . This is determined by the power-law region between t 1 and t 2 in the FPT distribution that still exists for α > 0, as discussed in Appendix A.

Preferences during the diffusion process
It is interesting to compare the preferences of the particles to be absorbed by a particular wall (the ultimate preferential direction at the end of the trajectory) to their preference for a particular side at earlier times. The problem can be formulated in different ways.
First, we can take a sample of particles, calculate the total combined time spent by them on a particular side, and divide by the total combined duration of all trajectories. Obviously, this is the same as the average time that a trajectory spends on a particular side divided by the total average duration of the trajectory. Using again the MC I algorithm with all moves successful, we note that the average number of steps in each segment is the same for left and right segments. Denoting the average total number of steps n and given that the fraction of left segments is p ′0 I− , as discussed in the previous section, the average time a trajectory spends on L is similarly the average time on R is and then the fractions of the time on L and R are These functions for α = 0, 0.5, and 1 are plotted in Fig. 9(a) together with the results of MC I simulations in these cases, and agreement is excellent. Note that these probabilities are different from the probabilities to reach a particular wall, π L,R , given by Eqs. (125) and (126). In particular, while for Ito π L = π R , f L < f R when η L < η R . Conversely, for isothermal f L = f R , but π L > π R for η L < η R . In fact, π L for Ito equals f R for isothermal, and vice versa, or, in general, Second, we can also determine the fraction of the time spent on a particular side for each trajectory in the sample and then average these fractions. We denote these averaged fractions φ L,R . This definition gives more weight to short trajectories than the previous one, since a short trajectory contributes little to both the total duration of trajectories and to the duration on a particular side, the quantities that enter the definitions for f L,R , but it contributes as much as long trajectories to φ L,R . The results for φ L,R obtained in MC I simulations for α = 0, 1/2, and 1 are shown in Fig. 9(b). These results can be fitted using the approximate expressions, If the equalities (145) and (146) were exact, then, (1) φ L for Ito would be equal to φ R for isothermal, and (2) φ L and φ R for Stratonovich (α = 1/2) would be exactly 1/2 for all η R /η L . However, we have checked using several different lattice sizes and extrapolating to the infinite size (results not shown) that these relations are indeed approximate in the continuum limit and both statements (1) and (2) are approximate as well.
Comparing Eqs. (125), (126), (141), (142), (145), and (146), we note that the average fractions of time spent on the low-viscosity side calculated in both ways are always smaller than the probability to be absorbed by the wall on the low-viscosity side, but the second way of calculating the average gives a result intermediate between the other two. That is, for η R /η L > 1, The differences between these quantities are particularly striking when the viscosity ratio is large, in which case both "<" signs in Eq. (147) turn into "≪", so there are ranges of α where one or two of the quantities f L and φ L approach zero while π L approaches unity. These relations are easy to explain by recalling that many particles that are eventually absorbed by the wall on the low-viscosity side linger until time ≫ t 1 (up to ∼ t 2 ) necessarily spending much of that time on the high-viscosity side, which explains why f L , φ L < π L . In addition, such long trajectories have a relatively larger weight when calculating f L than when calculating φ L , which explains f L < φ L . Note that nearly all trajectories reaching the wall on the high-viscosity side spend most of the time on that side.

VII. DISCUSSION
In this paper, we have developed two dynamical Monte Carlo algorithms for simulating the dynamics of a particle in the presence of sharp viscosity (or, in general, diffusivity) changes applicable to cases where different interpretations of the stochastic term in the corresponding overdamped Langevin equation (or calculi, such as Ito, Stratonovich, and isothermal) are appropriate. The validity of these algorithms was demonstrated by the good agreement between the approaches and with Brownian dynamics for the Ito calculus and Langevin dynamics for the isothermal calculus, as well as with theoretical predictions. We have argued that MC approaches, such as those developed here, can be particularly advantageous when a sharp interface is present in the system, since solving the overdamped Langevin equation directly is problematic in this case, except for the Ito calculus (corresponding to BD), while the LD approach is only applicable to the isothermal case and generally more computationally intensive.
While our simulations of particles on a 1D interval with the interface in the middle were intended primarily to validate our approaches, the system itself turned out to be rather interesting, with sometimes non-obvious properties that depend significantly on the choice of the calculus characterized by the parameter α that can vary between zero (Ito) and one (isothermal). The differences between these extreme cases are summarized in Table I for a number of features of the system. Not only do the features vary greatly between calculi, but there are unexpected symmetries between certain features (indicated by arrows in the table). Clearly, when modeling a physical system, the choice of calculus in deciding how to treat the interface must be done with great care.
The goal for the future is extending the approaches presented here to more complex and physically interesting situations, such as two-and three-dimensional systems and biased diffusion under the action of a force. We currently work on applying these algorithms to studying the translocation of a polymer through a hole in a membrane, with fluids of different viscosities on the two sides of the membrane. Of course, extension to the case of several interfaces is straightforward and a gradual change in viscosity can most likely be modeled as many interfaces one lattice site apart. Besides this, even in the simplest case considered here, the algorithms can be further improved: There are several free parameters that we have chosen arbitrarily, but that can be optimized to speed up the convergence to continuum, so even coarse discretization produces accurate results.
It is useful to note that the physical scenario examined in this paper actually belongs to a more general class of problems. Of course, the discrimination between different calculi arises for any stochastically driven process that can be written as an overdamped Langevin equation in which the prefactor of the noise term is statedependent [3,5,12,23,52,53]. Hence, while we have focused on diffusion in an inhomogenous medium, similar concepts can be applied to many problems including financial modeling [54,55], NMR spectroscopy [57], population growth models [56], and climate modeling [58][59][60]. Another interesting example is the work by Smythe et al. [61] examining noise-induced phase transitions where it was found that digital simulations produced results in agreement with the Ito formulation while apparently equivalent analog simulations produced results that agreed with the Stratonovich formulation. Given these diverse scenarios, the simulation approaches introduced herein may be useful outside of the problem for which they were originally intended.

VIII. ACKNOWLEDGEMENTS
This research was supported by grants from NSERC and the University of Ottawa. Simulations were performed on the SHARCNET computer system.  Fig. 7, consider the particle distribution function ρ(x, t) in this regime for the problem with absorbing boundaries. Similar to the case of reflecting boundaries discussed in Sec. VI A, ρ(x, t) consists of a quasistationary distribution in the low-viscosity region and a "Gaussian-like" distribution in the high-viscosity region. (The exact meaning of "Gaussian-like" will be revealed later.) The difference is that because of the absorbing boundary conditions [ρ(x, t) = 0 at the walls], the quasistationary part is no longer a constant, but instead a linear function, Here and in what follows we assume that the left half of the system is the low-viscosity region (η L ≪ η R ). Equation (A.1) provides a relation between the value of ρ(x, t) at −0 and its derivative: The slope, [∂ρ(x, t)/∂x]| x=−d/2 ≈ [∂ρ(x, t)/∂x]| x=−0 , is related to the flux to the left wall and thus to the conditional FPT distribution for the left wall, From Eq. (A.2), using the jump condition (44) and the flux continuity condition (46), we obtain a similar relation on the high-viscosity side of the interface: For the diffusion problem in the high-viscosity region, this serves as a boundary condition (BC). BCs of this type, where the derivative is proportional to the value of the function, are known as Robin boundary conditions [62]. At times t ≪ t 2 the right boundary can be ignored, and we can treat this problem as being defined on a semiinfinite interval,   [63]. We will opt instead for a more intuitive approach, and from now on, we will omit all numerical constants of order unity from our considerations. Since the BC (A.4) is intermediate between the absorbing BC ρ(+0, t) = 0 and the reflecting BC [∂ρ(x, t)/∂x]| x=+0 = 0, we can expect that in some limits the former can be approximated by one of the latter two. For the reflecting BC, with particles starting near +0 at time t = 0, the solution is where B is a constant. At a given time, the maximum value of this function is at x = 0 and equals B; the x derivative at that point is zero. On the other hand, the maximum absolute value of the x derivative is at The ratio between the maximum absolute value of the derivative and the maximum value of the function is η R /t. If this value is much larger than the ratio of the values of the derivative and the function at +0 given by Eq. (A.4) [this latter ratio is ∼ (1/d)(η R /η L ) α ], then just a small shift in space of the solution of the reflecting problem, Eq. (A.6), will be sufficient to satisfy the BC (A.4), so in this situation, which is observed when the solution of the reflecting problem, Eq. (A.6), is an adequate approximation of the solution of the original problem with the BC (A.4). On the other hand, considering the solution of the absorbing problem, where C is a constant, the maximum value of the function is , the maximum value of the derivative is C/t 3/2 , and the ratio of the latter and the former is η R /t, as in the reflecting case. But now, in order for the solution with the BC (A.4) to be obtainable by a small perturbation of that of the absorbing problem, we need this ratio to be much smaller than that specified by Eq. (A.4), which happens for t ≫ t c , with t c still given by Eq. (A.7). Thus, the solution of the Robin boundary-value problem (A.5)-(A.4) crosses over from being "reflecting-like" at t ≪ t c to being "absorbinglike" at t ≫ t c . Note that the "reflecting-like" solution corresponds to the particle number on R staying nearly constant, whereas the "absorbing-like" solution corresponds to this particle number decreasing rapidly (as t −1/2 ) as the particles are "absorbed" by the interface, "transferred" to the left side, and "flow" to the wall.
However, it should be remembered that our consideration is only valid in the interval t 1 ≪ t ≪ t 2 . When α = 0, t c ∼ t 2 , so the "reflecting-like" solution should be observed in the whole range of validity, although some "traces" of the "absorbing-like" solution can appear as one approaches t 2 . On the other hand, when α > 1/2, t c < t 1 , so the solution is "absorbing-like" in the whole interval. Of course, since for t < t 1 the particles do not yet feel the walls, the distribution should be the same as in the reflecting case, that is, "reflecting-like"; the transition should occur around t 1 . It is only for 0 < α < 1/2 that the crossover should actually be observed in the interval t 1 ≪ t ≪ t 2 . If we choose t ∼ √ t 1 t 2 , then for α = 0 the particle distribution in the high-viscosity region should resemble Eq. (A.6); for α = 0.5 and α = 1 it should be similar to Eq. (A.8); finally, for α = 0.25 this time equals t c and the distribution should be a "hybrid" of these two cases. Figure A.1 showing the distributions for η L = 2, η R = 2000, t = 10 5 indeed matches these expectations. Even though t = 10 5 is about 2.5 times larger than √ t 1 t 2 ≈ 40000 (this larger time is chosen to have broader distributions, so their features are more visible), the crossover is broad enough that this time can still be considered to correspond to the crossover for α = 0.25. Note that the distributions in the low-viscosity region are linear, as given by Eq. (A.1); note also that a different scale had to be used for the low-viscosity region, since the values of the particle density in that region are much lower than in the high-viscosity region. Thus, for any α, the fraction of particles in the low-viscosity region is very small in this time range; contrast this with the reflecting case, where, for example, in the isothermal case there are more particles in the low-viscosity region than in the high-viscosity region at all times (except when t → ∞, these numbers become equal). To obtain the plot in Fig. A.1, the variant of the MC I algorithm with the waiting time drawn from the distribution (72) is used providing significant computational time savings for this large viscosity ratio.
We are now in a position to obtain the FPT distribution for the left wall, P L (t). In the "reflecting-like" situation, considering the solution with the reflecting BC, Eq. (A.6), by itself, the x derivative at +0 is, of course, zero. The next iteration is obtained by inserting the value of the function at this point, B/ √ t, into the BC (A.4), which gives ∼ B(η R /η L ) α /(d √ t). Using Eqs. (46) and (A.3), Here we have used the fact that π L is always ∼ 1. In the "absorbing-like" case, the x derivative at +0 is C/t 3/2 , and in the same way we obtain Thus, the P (t) dependence should always be a power law; the exponent can be −1/2 in the "reflecting-like" situation or −3/2 in the "absorbing-like" situation. In Fig. A.2, we plot P L (t) for η L = 2 and η R = 2000. We see that indeed, for 0.5 ≤ α ≤ 1 (which corresponds to the "absorbing-like" situation at all times), the dependence is ∼ t −3/2 . For α = 0 (the "reflecting-like" case), the t −1/2 fit is adequate, although the slope does get steeper close to t 2 , as the crossover is starting to be "felt." For a smaller viscosity ratio, as in Fig. 7, the crossover is felt over the whole intermediate range of times, which is probably what causes the apparent t −0.8 dependence in that case. Finally, for α = 1/4 the curve cannot be fitted by a simple power-law function; apparently, the crossover is so broad that the already large viscosity ratio needs to be increased by at least another order of magnitude before the distinct t −1/2 and t −3/2 regions can be observed.
To obtain the prefactors (or the coefficients B and C) in P L (t), we assume that the number of particles on R, N R , at the beginning of the intermediate time interval (t ∼ t 1 ) is still the same as at t → 0, which, in its turn, is the same as in the reflecting case, since at t → 0 the particles have not reached the boundaries yet. Note that this does not mean that the ratio N R /N L at t ∼ t 1 is the same as at t → 0, because particles on L do start to be absorbed around t 1 , so N L does change by ∼ t 1 . This means that while Eq. (110) can be used to get N R at t ∼ t 1 , N L needs to be replaced with N 0 − N R , where N 0 is the total initial number of particles on both sides. Two cases need to be considered. For α < 1/2, Eq. (110) (modified as above) gives N R ≫ N 0 − N R , so On the other hand, for α ≥ 1/2, The first case is "reflecting-like" at t 1 and until t c ; normalizing the particle distribution function so the area under the whole curve (L and R) is unity initially [and therefore, given Eq. (A.11), the area under R is still approximately unity at ∼ t 1 ], we obtain (B/ √ t 1 ) √ D R t 1 ∼ 1, thus B ∼ 1/ √ D R ∼ √ η R , and using Eq. (A.9), After the crossover, P L (t) ∝ t −3/2 ; by matching to Eq. (A.13) at the crossover point t c , (A.14) The case α ≥ 1/2 is "absorbing-like" for all intermediate t; the area under the curve given by Eq. (A.8) at , so, according to Eq. (A.12), the normalization condition is C/(η R t 1/2 , and, using Eq. (A.10), which is the same expression as Eq. (A.14).
Note that Eqs. (A.13) and (A.15) are not necessarily valid, even approximately, at t ≈ t 1 , mostly because around t 1 , and until time t * = βt 1 , where β is a sufficiently large factor, the particle distribution on L may still be very different from that given by Eq. (A.1) .
[Strictly speaking, this means that using t ∼ t 1 in considerations leading up to Eqs. (A.13)-(A.15) is not entirely correct, but since for all practical purposes β can be chosen as a large, but finite and constant numerical factor, this does not affect the scalings, apart from at most logarithmic corrections.] For α > 1/2, the number of particles on R relative to the total initial number of particles in the system N 0 is always ≪ 1. On the other hand, using Eq. (A.3), the number of particles on L (again, relative to N 0 ) at t = t * (at which the above considerations are already valid but still as close to t 1 as possible) is (A.16) Thus by time t * most particles have already disappeared from the system. This absorption of particles must have happened around time ∼ t 1 over the time interval ∼ t 1 . Thus, P L (t 1 ) ∼ 1/t 1 . However, Eq. (A.15) predicts P L (t 1 ) ∼ (1/t 1 )(t 1 /t 2 ) α−1/2 ≪ 1/t 1 . So, on top of the power law predicted by Eq. (A.15) there should be an additional peak (or "bulge") at ∼ t 1 of height ∼ 1/t 1 : This "bulge", in fact, contains nearly all particles, although the rest of the distribution is still important when calculating the MFPT, as we will see. On the other hand, for α < 1/2 only a small fraction of all particles, ∼ (η R /η L ) α−1/2 , are on L initially. Even if all of them are absorbed over time ∼ t 1 , the resulting FPT density is , which is the same as what Eq. (A.13) at t = t 1 gives, so there is no additional "bulge". In fact, in Fig. A.2 there appear to be "bulges" for all α. However, it should be kept in mind that our discretization is not entirely adequate for this problem with a very high viscosity ratio: at t = t 1 = 1250, the width of the distribution on R is less than one, that is, less than one mesh step. In the same figure, we also show the results obtained for d = 100, with the resulting distributions appropriately rescaled. The "bulges" on the corresponding curves with α < 1/2 are smaller, so it is reasonable to assume that they largely disappear in the continuum limit. On the other hand, the "bulges" on curves with α > 1/2 are of the same size for d = 50 and d = 100 and therefore are "real", rather than discretization artifacts.
Using our results for FPT, we can now analyze the expression for the MFPT we have obtained [Eq. (138)], by decomposing the MFPT into parts associated with different regions of the FPT distribution. For α > 1/2, the weight of the bulge is nearly one, while the weight of the power-law region is On the other hand, the average time over the "bulge" is ∼ t 1 , but the average time over the power-law region is  which is much larger. The total MFPT is (A.20) where the first term comes from the "bulge" and the second one comes from the power-law region at intermediate times. Note that the second term is always larger than the first one (except for α = 1, when they are equal), and therefore it is the second term that determines the character of the MFPT dependence, which indeed coincides with Eq. (138). For α < 1/2, there are two power-law regions. The weight of the first, "reflecting-like" region is The weight of the second, "absorbing-like" region is So both regions carry equal weight. In fact, the largest contribution comes from the region around t c -the same result is obtained if the lower limit of the first integral and the upper limit of the second integral are replaced by arbitrary values, as long as the first one is much smaller than t c and the second one is much larger than t c . To put it differently, of all intervals [t a , t b ] with equal t b /t a , the one containing t c will contribute the most. However, other regions cannot be ignored when calculating the MFPT. The average time over the "reflecting-like" region is The resulting MFPT is The contribution of the second term dominates (except for α = 0, when these contributions are equal, even though the second region "shrinks into a point"). Again, this second term coincides with Eq. (138). Note that for both α < 1/2 and α > 1/2, contributions of the parts of the distribution with t > t 1 are significant. Particles spending such a long time before being absorbed necessarily spend most of that time in the high-viscosity region. Note in this regard that the term in Eq. (134) associated with the time spent on R dominates for 0 ≤ α < 1 when η R /η L → ∞, and even at α = 1 the terms corresponding to the time spent on L and on R contribute equally (by order of magnitude). It is interesting to compare the particle distributions in Fig. A.1 with those obtained for the same η L , η R , d, and t for reflecting boundaries. The crucial difference is that, while in the absorbing case the flow of particles is from the high-viscosity to the low-viscosity region, the opposite is true in the reflecting case, which follows from the fact that for η R /η L > 1 and for all α, the ratio N R /N L grows with time [see Fig. 5(c)], yet N L +N R stays constant. As mentioned, this growth is proportional to √ t; given that it occurs between t 1 and t 2 and using the short-time [Eq. For α < 1/2 this ratio is much larger than unity at all times, and the fraction of particles on R is The derivative of this is equal to the flux of particles into the right region, which in its turn is proportional to the space derivative of the particle distribution function at x = +0. Thus, Assuming that the distribution on R can be approximated by one half of a Gaussian curve with area ∼ 1 [Eq. (A.6) with B ∼ 1/ √ D R ∼ √ η R ], the maximum absolute value of the x derivative is ∼ η R /t. This is much larger than the absolute value of the derivative given by Eq. (A.29) for any t ≫ t 1 . Thus, arguing as in the absorbing case, the Gaussian with the peak at zero [Eq. (A.6)] is indeed a good approximation in this case.
Consider now α > 1/2. In this case, two situations are possible: for (A.30) N R /N L is still large, and the considerations above for α < 1/2 mostly apply. In particular, Eq. (A.29) and the expression for the maximum absolute value of the x derivative of the trial Gaussian are still valid, and the absolute value of the former is much smaller than the latter for any t ≫ t cr , so the Gaussian with a peak at zero still approximates the distribution well. On the other hand, for t ≪ t cr , N R /N L ≪ 1. This means that N L stays nearly constant and thus ρ(x, t) on L, besides being nearly constant in space, is also nearly constant in time and Then, according to the jump condition (44), and is likewise time-independent. The solution of the diffusion equation satisfying the condition ρ(+0, t) = const for t > 0 is the complementary error function erfc [64], and we get Note that N R ∼ √ t, as expected, and also that the inflection point is now at x = +0 (indeed, ∂ρ/∂t = 0 at x = +0 implies ∂ 2 ρ/∂x 2 = 0).
Similar to the absorbing case, we can consider the particle distributions at t ∼ √ t 1 t 2 . We expect the Gaussian distribution for α between 0 and 0.5, the complementary error function with the inflection point at 0 for α = 1 and a "hybrid" situation for α = 0.75. In Fig. A.3, we show the results obtained using MC I simulations with t Iw given by Eq. (73). These results confirm our expectations. The α = 0.75 curve resembles the α = 1 curve, but in the former case the inflection point is clearly shifted to the right of the interface. Comparing the results in Figs. A.1 and A.3, it is interesting to note a kind of "mirror symmetry" between them. For example, in the absorbing case the distributions with α between 0.5 and 1 are very similar in the high-viscosity region, and in the reflecting case the distributions with α between 0 and 0.5 have this property. In the absorbing case, α = 0.25 is "hybrid," while in the reflecting case, α = 0.75 is "hybrid." Finally, in the absorbing case the distribution for α = 1 is given by Eq. (A.8), which is the derivative of the Gaussian distribution, and in the reflecting case the distribution for α = 1 is given by the complementary error function [Eq. (A.33)], which is the integral of the Gaussian distribution.
Appendix B: Proportion of the number of steps on left and right for a random walk absorbed at a particular wall Using the MC I algorithm with all move attempts successful (Sec. IV A 2) reduces the problem of particle dif- All particles start at the interface in the middle of the system. fusion in a system with a diffusivity interface to a lattice random walk that is unbiased everywhere, except perhaps at the interface site. Depending on what side of the interface the particle is on, a single step of this random walk takes a different amount of time, but this can be taken care of separately, by working in terms of the numbers of steps on each side and then multiplying by the respective time steps. In this way, our problem is reduced to that of an ordinary lattice random walk with a constant diffusion rate that starts at the site in the middle between two absorbing walls and can be biased at that site, but nowhere else. We first consider the case with no bias even at the middle site (corresponding to the Ito calculus) and then generalize. First, consider the diffusion on a lattice with mesh step a = 1 of an ensemble of particles released half-way between walls at x = ±d/2, with d/2 integer and the initial position of each particle given by x 0 = 0. By symmetry, the average mean first passage times (MFPT) to either wall (the conditional MFPTs) are equal and thus equal to the MFPT to either wall (the unconditional MFPT). If the time is measured in terms of the number of steps, then the MFPT is the average number of steps, which is [50] n = d 2 4 . (B.1) To proceed, we can decompose the full trajectory of a particle into segments separated by its visits to the middle site. Each such segment is entirely either in the left half or in the right half and for brevity, we refer to these groups of segments as left and right segments. Consider the ensemble of trajectories ending on the right wall. Taking one such trajectory and replacing all of the right segments (except for the last one, of course) with its mirror-image left segment and vice-versa produces another valid trajectory belonging to the ensemble. Therefore, up to the last segment, there should be an equal number of left and right segments on average and they should have the same average duration. Thus, we can divide the trajectory into two parts. The first part contains all but the last segment and on average, the time spent in both halves of the interval is the same, but the second part, consisting of only the last segment, is entirely in the right half.
The average time corresponding to the second part of the trajectory can be obtained by recasting this problem as an equivalent one. Consider the interval consisting of just the right half of the original interval. If both boundaries of this interval are treated as absorbing, the particle is released right next to the boundary corresponding to the middle site and the conditional MFPT, given that the particle is absorbed by the other boundary, is calculated, then all trajectories which return to the middle site will be discarded. While the probability of being absorbed by a wall vanishes in the limit when the distance from the point where the particle is released to the opposite wall approaches zero (on a lattice this probability is ∼ a/d), the MFPT to the far wall is well-defined and is, for the continuum problem [65], where l = d/2 is the distance between the walls (the length of the half-interval under consideration). Using Eq. (55) with p s = 0, a = 1 and τ = 1, we convert this result into that for the average number of steps in the corresponding lattice problem,