## Abstract

Consider a particle diffusing in a confined volume which is divided into two equal regions. In one region, the diffusion coefficient is twice the value of the diffusion coefficient in the other region. Will the particle spend equal proportions of time in the two regions in the long term? Statistical mechanics would suggest yes, since the number of accessible states in each region is presumably the same. However, another line of reasoning suggests that the particle should spend less time in the region with faster diffusion, since it will exit that region more quickly. We demonstrate with a simple microscopic model system that both predictions are consistent with the information given. Thus, specifying the diffusion rate as a function of position is not enough to characterize the behaviour of a system, even assuming the absence of external forces. We propose an alternative framework for modelling diffusive dynamics in which both the diffusion rate and equilibrium probability density for the position of the particle are specified by the modeller. We introduce a numerical method for simulating dynamics in our framework that samples from the equilibrium probability density exactly and is suitable for discontinuous diffusion coefficients.

## 1. Introduction

Consider a particle diffusing in a two-dimensional box with reflecting boundary conditions. We show a portion of a simulated trajectory of such a system in figure 1. Suppose that in the left half of the box the particle diffuses with coefficient *D*_{1}, and that in the right half of the box the particle diffuses with coefficient *D*_{2}=2*D*_{1}. We assume that there are no external forces acting on the particle. Our question is: does the particle spend an equal fraction of time on each side of the box in the long run?

One answer is based on statistical mechanics.

*Statistical mechanics prediction:* the particle will spend an equal proportion of time on each side of the box.

The justification for this prediction is the principle of statistical mechanics which states that ‘An isolated system in equilibrium is equally likely to be in any of its accessible states’ (Reif 1965, p. 54). Since all states in the box are accessible, and there are an equal number of states on each side of the box, the particle should spend an equal proportion of its time on each side of the box.

Another answer is based on the idea of rescaling time in one side of the box.

*Time-change prediction:* the particle will spend less time on the side of the box where the diffusion coefficient is greater.

The justification for this prediction is that, in the absence of any drift, faster diffusion is equivalent to time passing more quickly. This means that the periods of time the particle spends on the right-hand side of the box will be shorter than those spent on the left-hand side. Hence the total time the particle spends on the right-hand side of the box will be less. We show how this prediction is a straightforward consequence of interpreting the particle's motion as drift-free diffusion, where we interpret the state-dependent diffusion coefficient using the Itô convention (Gardiner 2004). We can write the equation for the particle's motion as
1.1where *x*=(*x*_{1},*x*_{2}) and *B* is standard two-dimensional Brownian motion. (Equivalently, we may write this equation as d*x*/d*t*=*b*(*x*)*η*(*t*), where *η* is two-dimensional Gaussian white noise.) We specify *b*(*x*)=*b*_{1} for *x*_{1}<0 and *b*(*x*)=*b*_{2} for *x*_{1}>0. Here , *i*=1,2 where *D*_{i} is the corresponding diffusion coefficient. We enforce reflecting boundary conditions at the four walls of the box. The (Itô-)Fokker–Planck equation for *ρ*(*x*,*t*), the probability density of the particle's position at time *t*, is (Gardiner 2004, p. 118)
The equilibrium density *ρ*_{eq}(*x*) satisfies
Reflecting boundary conditions for the diffusion correspond to Neumann boundary conditions (zero-flux) for the Fokker–Planck equation. With these boundary conditions, the unique equilibrium density is (*ρ*_{eq}(*x*)=*C*/*D*(*x*)) for some constant *C*. Thus, since *D*_{2}=2*D*_{1}, the particle spends half as much time on the right-hand side of the box as on the left. (We discuss the relation of our question to other interpretations of (1.1) in §4.)

Neither the statistical mechanics prediction nor the time-change prediction are definitive. The statistical mechanics prediction relies on the principle of equal probability of all accessible states, which needs to be independently justified for the mesoscopic level of description we are considering here. The time-change prediction is not definitive since Itô stochastic differential equations (SDEs) are themselves mesoscopic models whose use can only be rigorously justified by showing how they arise as the coarse-scale limit of microscopic dynamics. Since the two predictions contradict each other, at least one of them must be wrong for any given physical system. The approach by which we resolve this apparent contradiction is to study a microscopic model of the box system we have described above, and then see what proportion of the time the particle spends on each side in simulations of that system. If one prediction turned out to always be true for our model system, we could use our result as a baseline for investigating under which more general situations the prediction was still true. However, we will see that even for our simple system, there are parameter choices that make the statistical mechanics prediction correct and parameter choices that make the time-change prediction correct. This demonstrates that there is no *a priori* reason to conclude that one prediction or another is correct, given only the diffusion coefficient *D*(*x*) for the system.

In §2, we describe our model system, perform numerical experiments on it, and show that the result of the numerical experiments can be determined analytically from the properties of the model system. The system we consider is a deterministic Hamiltonian billiard system: the random Lorentz gas (Dettmann 2000). The system has two free parameters: disc radius and free volume fraction. There is a one-parameter family of values of these parameters that can generate the same effective diffusion coefficient. We will show how to choose these parameters to obtain arbitrarily good approximations to diffusive motion for the particle on each side of the box. Using the degree of freedom in the choice of parameters, we show that the equilibrium density of the particle is underdetermined by the diffusion coefficient on each side of the region. For a given *D*_{1} and *D*_{2}, exploiting the flexibility in the parameter choice allows us to create systems in which either the statistical mechanics prediction or the time-change prediction is correct.

Our justification for studying a particular microscopic system is threefold. Firstly, there is a solid analytical understanding of the random Lorentz gas on which we can base our simulations. Secondly, our purpose is not to conclude that a particular style of mesoscopic modelling is always the correct one, but to show that postulating a state-dependent diffusion rate is not enough to fully specify mesoscopic behaviour. For this objective, it is enough to show that multiple mesoscopic behaviours are possible for a single class of simple models, as we do in this paper. Finally, we have resorted to a purely deterministic microscopic model, rather than a stochastic microscopic model (such as Langevin dynamics or a random walk) to avoid concerns that the method used to introduce randomness at the microscopic level somehow biases the results at the mesoscopic level. This last point distinguishes our work from similar discussions of van Kampen (2007, §XI.3) and Othmer & Stevens (1997, §2). These authors show that spatially inhomogeneous random walks, with natural choices of parameters, can lead to either prediction holding true at the mesoscopic level.

Similarly, the conclusions of our §2 closely parallel those of Korabel & Barkai (2011), in which the same question is answered using a random-walk model on a one-dimensional lattice. There, following the earlier work of Ovaskainen & Cornell (2003), they show that choosing how the random walk behaves at the interface of two regions of differing diffusion coefficient leads to different equilibrium probability densities for the system. Korabel & Barkai (2011) explain how to determine the correct interface behaviour of the random walk model using experimentally measurable quantities. Our approach differs in that, in our model, the interface behaviour is determined indirectly through the microscopic dynamics that we describe.

The results in §2 show that choosing a diffusion coefficient *D*(*x*) is not enough to fully specify diffusive dynamics in the absence of other assumptions. In particular, our results show that fixing *D*(*x*) is not enough to specify *ρ*_{eq}(*x*), the equilibrium probability density for the position of the particle. In §3, we propose a framework for modelling state-dependent diffusion that makes this fact explicit. Rather than simply specifying a state-dependent diffusion coefficient, we specify diffusion coefficient *D*(*x*) and an equilibrium density *ρ*_{eq}(*x*) which together with a detailed balance assumption (no-flux in equilibrium) completely determine the dynamics. We then introduce a new method for the numerical simulation of diffusive dynamics which makes use of our framework: the numerical method is expressed in terms of *D*(*x*) and *ρ*_{eq}(*x*) and makes no reference to a drift term. The method consists of Euler–Maruyama steps for a purely diffusive Itô SDE together with Metropolis rejections. The method is similar to the Metropolis-adjusted Langevin algorithm (MALA) (Roberts & Tweedie 1996; Bou-Rabee & Vanden-Eijnden 2010), except that there is no drift term in the Euler–Maruyama step, and it is the Metropolis rejections that induce any drift in the trajectories. The advantages of our method over the Euler–Maruyama method are that it samples space with the correct equilibrium density *ρ*_{eq}(*x*) and performs well even when *D*(*x*) and *ρ*_{eq}(*x*) have discontinuities.

One area in which the proposed framework in §3 could be used is cellular biology. Although earlier models of chemistry in the cytoplasm of cells assumed that chemical species were ‘well-mixed’ and thus ignored diffusion, more recent models have taken the geometry of the cell and the diffusion coefficient of various molecules into account (Turner *et al.* 2004). Effective diffusion coefficients of a molecule in a cell differ from that in water owing to the crowding effect of other molecules. One approach is to model the motion of a molecule as diffusion with constant coefficient, but then have the effective diffusion be modified by interaction with other particles which are also included in the model (Ridgway *et al.* 2008). Another is to not include the crowding particles in the model, but to model their effect with a modified diffusion coefficient (Hall & Hoshino 2010). Given the inhomogeneity of the cytoplasm, we can expect that this effective diffusion coefficient of the molecule (and its equilibrium probability density) will vary with location within the cell. Our framework and numerical method are developed with this latter situation in mind.

In §4, we conclude by explaining the relation between our results and three different interpretations of state-dependent diffusion: Itô, Stratonovich and the isothermal convention of Lançon *et al.* (2001).

## 2. A model system for state-dependent diffusion

In §1, we described a two-dimensional system consisting of a single particle diffusing inside a rectangular box and reflecting off the boundaries. The diffusion coefficient is twice as large on the right-hand side of the box as the left. In this section, we demonstrate how to construct a family of deterministic Hamiltonian systems that approximates this behaviour on a coarse scale.

In §2*a*, we describe the random Lorentz gas (Dettmann 2000), a deterministic system that when given a random initial condition yields constant-coefficient diffusion at a coarse scale. In §2*b*, we show how to approximate the model system of the introduction by creating two adjacent domains of the random Lorentz gas within a bounding box. In §2*c*, we describe numerical experiments with the box system demonstrating that the fraction of time a particle spends on each side of the box cannot be determined solely from the values of the diffusion coefficient. In §2*d*, we explain how the result of the numerical experiments in §2*c* can be predicted from properties of the dynamics of the microscopic system.

### (a) Random Lorentz gas

Consider infinitely many discs with positions fixed in . The centres of the discs are distributed randomly with uniform density subject to the constraint that the discs do not overlap. We consider a single point particle interacting with the discs in the following manner. Given an initial velocity and an initial position not on a disc, the particle moves with constant velocity until it meets a disc. Then it undergoes an instantaneous elastic reflection with the boundary of the disc, with the angle of incidence equalling the angle of reflection. This model is the two-dimensional *random Lorentz gas* (Dettmann 2000), a mathematical formulation of a model originally due to Lorentz (Lorentz 1905). We will always consider the case when the particle has initial velocity of magnitude 1. Figure 2 shows trajectories of the random Lorentz gas at two different scales.

The random Lorentz gas has two parameters: the radius of the discs *r*, and the number of discs per unit area *λ*. The radius *r* can take any positive value. The density *λ* has a maximum value corresponding to the close-packed hexagonal pattern of discs. We define *ϕ* to be the free volume fraction, the proportion of the area not occupied by discs. We have that *ϕ*=1−*πr*^{2}*λ*. The free volume can take any value in [*ϕ*_{min},1) where , regardless of the value of *r*. The pair (*r*,*ϕ*) provides an alternative parametrization of the random Lorentz gas, with the advantage that *ϕ* is dimensionless.

For *ϕ*=*ϕ*_{min}, adjacent discs are touching and the particle remains in a small region of the plane for its entire trajectory. For *ϕ*∈(*ϕ*_{min},1), the probability that two discs are touching anywhere in the plane is zero (Dettmann 2000, p. 325) and motion of the particle is conjectured to be diffusive (Dettmann 2000). (This contrasts with the situation of overlapping discs for which there is a phase transition to anomalous transport for dense enough placement of discs (Höfling *et al.* 2008).) Specifically, if we start the particle at initial position *x*(0) not coincident with a disc and give it initial speed 1 and uniformly distributed direction in [0,2*π*), then *x*(*t*), the position of the particle at time *t*, is approximately distributed as a Gaussian random vector with mean 0 and variance matrix 2*DtI*. Here *I* is the 2×2 identity matrix and *D* denotes the diffusion coefficient. This conjecture is supported by analytical calculations (Ernst & Weyland 1971; Dettmann 2000) and numerical simulations (Bruin 1972; Dettmann & Cohen 2000).

A stronger and more formal statement of the conjecture is stated in the language of weak convergence (Billingsley 1999). Specifically, let *x*(*t*) denote the position of the moving particle at time *t*. Fix *x*(0) to be some point not coincident with a disc and let the initial velocity be chosen as above. It is conjectured that
as , where *B*(*t*) is standard two-dimensional Brownian motion (meaning 〈*B*(*t*)〉=0, 〈*B*(*t*)*B*(*t*)^{T}〉=*tI*), *D*_{r,ϕ} is the diffusion coefficient which depends on *r* and *ϕ*, and ⇒ denotes weak convergence in the space of continuous functions (Billingsley 1999). Such a result holds for the certain periodic Lorentz gasses (Bunimovich & Sinaĭ 1980/1981; Klages & Dellago 2000), but remains open for the random Lorentz gas. (We chose for our study not to use the standard periodic Lorentz gas with discs centred on a hexagonal lattice since for large enough *ϕ* the particle undergoes superdiffusive motion in this case.)

A scaling argument shows that, for fixed *ϕ*, *D*_{r,ϕ} is proportional to *r*. To see this, letting |·| denote the Euclidean norm, observe that the coefficient *D* can be obtained as , assuming *x*(0)=0. If we rescale space by a factor *R*, we increase both the size of the discs and the distance the particle travels by a factor of *R*, without changing *ϕ*. So 〈|*x*(*t*)^{2}|〉 increases by a factor of *R*^{2}. To maintain the speed of the particle as 1, we also have to rescale time, increasing *t* by a factor *R*. The net effect on the ratio 〈|*x*(*t*)|^{2}〉/4*t* is to increase it by a factor *R* (Sanders 2005, §3.1.6). So
2.1for some function *f* of *ϕ*.

Figure 3 shows the relation between *f*(*ϕ*) and *ϕ* that was computed using the techniques similar to those described in Dettmann & Cohen (2000). It appears that the function *f*(*ϕ*) is continuous on its domain, it is monotonically increasing, as and *f*(*ϕ*) goes to infinity as . Indeed, calculations from kinetic theory show that *f*(*ϕ*)∼3*π*/[16(1−*ϕ*)] in the limit (van Leeuwen & Weijland 1967; Bruin 1972).

Given any fixed diffusion coefficient *D*>0 there is a one-parameter family of choices of *r*,*ϕ* such that *D*=*D*_{r,ϕ}: for any *ϕ*∈(*ϕ*_{min},1), just choose *r*=*D*/*f*(*ϕ*).

### (b) Box with two domains

In order to investigate the main question of this paper, we take a rectangular box and divide it into two equal regions. Each region is filled with randomly placed discs, with different *r* and *ϕ* on each side. As in the Lorentz gas, the centres of the discs are placed uniformly at random with the condition that they not overlap with each other. Discs can intersect with the walls of the box but not the dividing line between the two sides of the box. Figures 4 and 5 show two examples of such configurations of discs. The dynamics of the point particle are the same as in the infinite random Lorentz gas, with the added condition that the particle reflects off the walls of the box.

Suppose we are given diffusion coefficients *D*_{1},*D*_{2}>0. We can choose *r*_{1},*ϕ*_{1} and *r*_{2},*ϕ*_{2} such that *D*_{1}=*D*_{r1,ϕ1} and *D*_{2}=*D*_{r2,ϕ2}. For small enough *r*_{1},*r*_{2}, the dynamics of the particle in this system will be well approximated by a particle that diffuses with coefficient *D*_{1} on the left-hand side of the box and diffuses with coefficient *D*_{2} on the right-hand side of the box. With appropriate choices of the parameters, we can investigate the question of the proportion of time the particle spends on each side of the box.

### (c) Numerical experiments

We start the particle off at some position in the box not on a disc with velocity of magnitude 1 and randomly chosen direction. The motion of the particle is simulated with an event-driven simulation, computing a trajectory that is accurate up to the errors of floating point arithmetic (Dettmann & Cohen 2000). Periodically the position of the particle is recorded. At the end of a long trajectory, we compute the number of times the particle is on the right-hand side of the box divided by the number of times the particle is on the left-hand side of the box.

We consider two choices of parameters on each side of the box, or set-ups, each of which leads to the effective diffusion coefficients satisfying *D*_{2}=2*D*_{1}. We have contrived the set-ups so that in set-up 1, the statistical mechanics prediction is correct, and that in set-up 2, the time-change prediction is correct. The parameters for each set-up are summarized in table 1.

In the first set-up *ϕ*_{1}=*ϕ*_{2}=0.5 and 2*r*_{1}=*r*_{2}=0.6. We show the position of the discs in the box in figure 4. Over a trajectory of length 5×10^{5} time units the ratio between the occupation times is approximately 1, as we show in table 1. This result agrees with the statistical mechanics prediction.

In the second set-up 2*ϕ*_{2}=*ϕ*_{1}=0.60 and 8.3*r*_{1}≈*r*_{2}=0.75. We show the position of the discs in the box in figure 5. Over a trajectory of length 7×10^{6} time units the ratio between the occupation times is approximately , as we show in table 1. This result agrees with the time-change prediction.

Thus, by fixing the parameters appropriately, both the statistical mechanics prediction and the time-change prediction can be seen to be correct for the given mesoscopic behaviour.

### (d) Analysis and discussion

We first explain that, given the properties of random Lorentz gas, the results of the simulation in the previous subsection are predictable. We first note that the system is ergodic. The ergodicity of periodic Lorentz gas is shown by Sinaĭ (1970), since it is equivalent to dispersing billiards on the torus. Our system is not dispersing because the walls of the box are not convex, but Sinaĭ (1970, §9) shows how ergodicity still holds in this case by unfolding the box to obtain a periodic domain.

Ergodicity implies that the duration of time that the particle spends in a region of phase space is proportional to the volume of the region. For our system, this implies that the amount of time spent by the particle on each side is proportional to the free volume fraction *ϕ* on each side. For fixed *ϕ*, the parameters *r* and *D* are irrelevant to the proportion of time the particle spends on each side of the box. In set-up 1 above, *ϕ*_{1}=*ϕ*_{2} so the ratio of times spent on each side of the box is equal and the statistical mechanics prediction is correct. In set-up 2, *ϕ*_{1}=2*ϕ*_{2}, and so the particle spends twice as much time on the left-hand side of the box, and the time-change prediction is correct.

Despite appearances, the principles of statistical mechanics are not violated in set-up 2. There are two ways to reconcile the apparent contradiction. Firstly, one can say that at the microscopic scale statistical mechanics is not violated because there are not an equal number of states on each side. The number of states is proportional to the free volume fraction on both sides, and so the system spends more time where the free volume fraction is higher. The other way of reconciling the disagreement is at the mesoscale. Suppose we divide the box into many rectangular cells of equal area, each much smaller than the whole box, but much larger than the size of the discs. Each cell corresponds to a mesoscopic state which the particle may be in. Since the system is in the microcanonical ensemble (no exchange of energy with the outside) the probability of the system in equilibrium being in a particular mesoscopic state is determined by the entropy of that state, where the entropy of a mesoscopic state is proportional to the logarithm of the amount of microscopic states it contains. The greater free volume fraction on the left in set-up 2, implies that mesoscopic states have greater entropy there, and so the system will spend more time on the left-hand side of the box.

On the other hand, in set-up 1, the time-change prediction proves to be wrong. This means that the motion of the particle in the box is not well-described by the drift-free Itô SDE (1.1). Since, by the properties of the uniform random Lorentz gas, (1.1) *is* a good model for the dynamics of the particle within each of the two regions of constant disc radius *r*, it must be that the equation is no longer a good model at the boundary of the two regions.

We see that there are choices of *ϕ*_{1},*r*_{1} and *ϕ*_{2},*r*_{2} such that either the time-change prediction or the statistical mechanics prediction are correct, while still *D*_{2}=2*D*_{1}. Indeed, for any *D*_{1} and *D*_{2}, parameters can be chosen to induce arbitrary ratios between the times spent on the left-hand side and the right-hand side.

Although for generic values of the parameters in our Lorentz gas model neither prediction will be valid, we point out that the statistical mechanics prediction holds for a natural set of parameter settings, whereas the same is not true of the time-change prediction. For the time-change prediction to be correct, it is necessary for the ratio between *ϕ*_{1} and *ϕ*_{2} to match the ratio between *D*_{2} and *D*_{1}. We see no natural way that the parameters in our model may be set for this matching to occur. On the other hand, for the statistical mechanics prediction to be correct it is necessary that *ϕ*_{1} and *ϕ*_{2} be equal. There is at least one case where this condition approximately holds for a naturally occurring system. Consider a situation in which both *ϕ*_{1},*ϕ*_{2}≈1. This requires no fine tuning, only that the discs take up a small fraction of the total volume. Choosing *r*_{1} and *r*_{2} to be unequal leads to different diffusion coefficients on each side, but the particle still spends approximately equal proportions of time in each region. Likewise, in any physical system where a particle diffuses by interacting with small, sparsely placed scatterers, we expect the statistical mechanics prediction to be correct.

## 3. A proposal for modelling with state-dependent diffusion

The numerical experiments of §2*c* demonstrate that the equilibrium density *ρ*_{eq}(*x*) is not determined solely by the local diffusion rate *D*(*x*), even in situations with no external forces acting on the particle. This raises the practical issue of how to model systems with state-dependent diffusion. Ideally, mesoscopic diffusive models would be derived from microscopic models via an asymptotic technique, such as the van Kampen system-size expansion (van Kampen 2007, ch. XI.3). But in many circumstances, deriving a realistic microscopic model may be impractical.

Instead, in §3*a*, we describe a more phenomenological approach. We assume that the modeller posits an isotropic state-dependent diffusion rate *D*(*x*) and an equilibrium density *ρ*_{eq}(*x*). We then derive a drift coefficient *a*(*x*) for an Itô SDE with diffusion coefficient *D*(*x*) that gives the desired *ρ*_{eq}(*x*). In §3*b*, we show how the Lorentz gas system of §2 can be modelled at a mesoscopic level in this way. In §3*c*, we then describe an algorithm for simulating such SDEs that does not refer to *a*(*x*) but is expressed in terms only of *D*(*x*) and *ρ*_{eq}(*x*).

### (a) A modelling framework for state-dependent diffusion

We model the diffusion in *k* spatial dimensions by the Itô SDE
3.1where we will determine *a*(*x*) in terms of *ρ*_{eq}(*x*) and *D*(*x*). Here *B*(*t*) is standard *k*-dimensional Brownian motion. (In alternate notation, we write this equation as where *η* is *k*-dimensional Gaussian white noise.) The Fokker–Planck equation for this system is
where *J* is the probability flux. Since the problem is underdetermined as stated, we will stipulate that the probability flux vanishes in equilibrium, which is the same as detailed balance holding; see van Kampen (2007, ch. XI.4) for criteria on systems under which this condition holds. We choose *a*(*x*) so that *J*(*x*) is zero in equilibrium:
Solving for *a*(*x*) gives
3.2Thus given an equilibrium density *ρ*_{eq}(*x*) and diffusion coefficient *D*(*x*) the appropriate Itô SDE is
3.3The Fokker–Planck equation of (3.3) is
or simplifying,
3.4which is the special isotropic case of eqn XI.4.14 of van Kampen (2007). Note that if *ρ*_{eq}(*x*) is constant with respect to *x*, which corresponds to the statistical mechanics prediction being true, then *a*(*x*)=∇*D*(*x*) and (3.3) reduces to
On the other hand, to obtain a drift-free Itô SDE ((*a*(*x*)=0) in this framework) requires *D*(*x*)*ρ*_{eq}(*x*) to be constant in *x*, which means that *ρ*_{eq}(*x*) is determined completely by *D*(*x*).

### (b) Connection to the Lorentz gas model

We explain the connection between the framework of §3*a* and the microscopic Lorentz gas model of §2. Given an instance of our random Lorentz gas model with two domains, what are the corresponding functions *D*(*x*), *ρ*_{eq}(*x*) in (3.4), the mesoscopic equation for *ρ*?

Suppose on the left-hand side of the box the Lorentz gas model has parameters *r*_{1} and *ϕ*_{1} and on the right it has parameters *r*_{2} and *ϕ*_{2}. The diffusion coefficients on each side, *D*_{1} and *D*_{2} are determined by the relation (2.1). To determine the equilibrium probability density of the particle on each side, let 2*A* be the total area of the box, so that each side has area *A*. We know that *ρ*_{eq,i}, the probability density on side *i*, is proportional to *ϕ*_{i}, and thus *ρ*_{eq,1}/*ρ*_{eq,2}=*ϕ*_{1}/*ϕ*_{2}. We also know that since the total probability must be 1, *ρ*_{eq,1}*A*+*ρ*_{eq,2}*A*=1. Solving for *ρ*_{eq,i} gives
for *i*=1,2. We define *D*(*x*) and *ρ*_{eq}(*x*) for *x*=(*x*_{1},*x*_{2}) in the box by
These functions determine *a*(*x*), the drift coefficient in (3.1), via (3.2). The drift *a*(*x*) is zero everywhere except along the line *x*_{1}=0 where it is not defined.

We determine the appropriate boundary conditions for *ρ* along the boundary line *x*_{1}=0. In order for the right-hand side of (3.4) to be well defined, we require *ρ*(*x*)/*ρ*_{eq}(*x*) to be continuous. So for any *x* along the boundary line we must have
where *x*^{−} denotes taking the limit from the left, and *x*^{+} denotes taking the limit from the right. For our particular choice of *ρ*_{eq}, this gives
The second boundary condition comes from assuming continuous flux across the boundary line:
For our particular choices of *D* and *ρ*_{eq}, since *ρ*_{eq} is constant away from the line *x*_{1}=0, this gives the boundary conditions

A possible direction for further investigation is to consider Lorentz gas models where disc radius *r* and free volume fraction *ϕ* vary smoothly with *x*, and to determine what *D*(*x*) and *ρ*_{eq}(*x*), and hence *a*(*x*) are in this case.

### (c) Numerical simulation

If both *ρ*_{eq} and *D* are smooth then the Euler–Maruyama scheme or Milstein's method is effective for simulating (3.3) (see Higham 2001). If either *D*(*x*) or *ρ*_{eq}(*x*) are discontinuous with respect to *x*, then the drift term will have a singularity that is not resolved by standard time-integration schemes, and numerical simulations may not correctly approximate the equilibrium density. We propose using a variant of MALA as introduced by Roberts & Tweedie (1996) and analysed by Bou-Rabee & Vanden-Eijnden (2010). MALA is obtained by taking a convergent numerical method for the SDE (in this case, Euler–Maruyama) and introducing Metropolis step-rejections in order to have a discrete-time process with the correct equilibrium density. Our approach is different in that we start with a convergent method for only the diffusive part of the Itô SDE, and then we introduce step rejections to induce the correct drift and equilibrium density.

Let *h* be the step length of our numerical discretization and let *X*_{n} be the numerical approximation to *X*(*nh*), where *X* is the solution to the Itô SDE (3.3) for *n*≥0. Given a numerical value *X*_{n}, we let the trial value for the next step be defined by
3.5where *B* is standard *k*-dimensional Brownian motion. Then *X*_{n+1} is given by
3.6where
and *ξ*_{k},*k*≥1 is an independent, identically distributed sequence of random variables, uniform on [0,1] and independent of *B*. Here
is the transition probability density for *X**_{n+1} being at *y* given that *X*_{n} is at *x*. The definition of *X**_{n+1} in (3.5) is just the Euler–Maruyama method for the Itô SDE . The Metropolis rejection procedure (3.6) accepts the Euler–Maruyama step with probability *α*_{h}(*X*_{n},*X**_{n+1}) and rejects it otherwise.

The Metropolis rejection procedure guarantees that the process *X*_{n},*n*≥1 has *ρ*_{eq}(*x*) as its equilibrium density. On the other hand, in any region of the state space where *D*(*x*) and *ρ*_{eq}(*x*) are constant with respect to *x*, *α*_{h} is 1 and so the method reduces to the Euler–Maruyama method for the constant coefficient diffusion without drift. Future work will investigate rigorously the convergence of the above method to the solutions of (3.3).

Here we numerically demonstrate the convergence of the method for the simple case where
and
Equation (3.3) with this choice of *D*(*x*) and *ρ*_{eq}(*x*) models a particle diffusing on the interval [−1,1] with reflecting boundary conditions at ±1 and a piecewise constant diffusion coefficient. The reflecting boundary conditions are conveniently implemented by our choice of *ρ*_{eq}(*x*).

Figure 6 shows the results of simulating (3.3) with these choices of *D*(*x*) and *ρ*_{eq}(*x*) using the method we have described. To show results, we divide the interval [−1,1] into 20 equal subintervals, and plot the density for the amount of time the particle spends in each subinterval. We also plot the effective diffusion coefficient for each bin, which we define to the average observed value of (*X*_{n+1}−*X*_{n})^{2}/2*h* over the trajectory, for all *n* such that *X*_{n} lies in the given bin. Simulations were conducted with *h*=0.01,0.001,0.0001 and for trajectories long enough so that standard statistical errors in the plots are smaller than the symbols used. We see that for all values of *h* the equilibrium density *ρ*_{eq}(*x*) is correctly reproduced. The effective diffusion coefficient converges to *D* as *h* goes to zero.

## 4. Discussion: Itô, Stratonovich or isothermal

We conclude by discussing our results in the context of the apparent ambiguity between Itô, Stratonovich and isothermal interpretations of stochastic integrals (Kupferman *et al.* 2004; Lau & Lubensky 2007; Volpe *et al.* 2010). For the purposes of discussion, we consider a particle moving in one spatial dimension whose position at time *t* is *X*(*t*). We assume that *X*(*t*) is a Markov stochastic process with continuous sample paths. We model the motion of the particle with the SDE
4.1where *B*(*t*) is standard Brownian motion. We call *a*(*x*) the drift and *b*(*x*) the diffusion of the SDE.

As is well known (Gardiner 2004; van Kampen 2007), unless we specify a particular interpretation, the SDE (4.1) does not unambiguously define the stochastic process *X*(*t*). To see this, we integrate (4.1) over [0,*T*] to get
The first term on the right has a unique interpretation as a Riemann integral, but the second term cannot be simply viewed as a Riemann–Stieljes integral, since *B* is not of bounded variation. If we compute the second term as the limit of Riemann sums, the answer depends on where in each subinterval the argument *b*(*X*(*t*)) is evaluated. For example, choosing *h*=*T*/*N*, *t*_{n}=*hn* and letting *B*_{n}=*B*(*t*_{n}) and *X*_{n}=*X*(*t*_{n}), suppose we take the integral with respect to *B* to be
where
Famously, unless *b*(*x*) is a constant, the limit depends on the choice of *α* (Volpe *et al.* 2010). If we choose *α*=0, we obtain the Itô interpretation of the integral, which yields a stochastic process *X*(*t*) with Fokker–Planck equation
If we choose *α*=1/2, we obtain the Stratonovich interpretation of the integral, which yields a process *X*(*t*) with Fokker–Planck equation
Note that the Fokker–Planck equation shows that the Stratonovich interpretation of the SDE with drift *a*(*x*) and diffusion *b*(*x*) yields the same stochastic process as the Itô interpretation of the SDE with drift *a*(*x*)+*b*(*x*)*b*′(*x*)/2 and diffusion *b*(*x*) (Gardiner 2004, p. 99). Finally, if we choose *α*=1 we obtain the anti-Itô or isothermal interpretation (Lau & Lubensky 2007; Volpe *et al.* 2010), which yields a process with Fokker–Planck equation
In this case, the isothermal interpretation of the SDE with drift *a*(*x*) and diffusion *b*(*x*) gives the same stochastic process as the Itô interpretation of the SDE with drift *a*(*x*)+*b*(*x*)*b*′(*x*) and diffusion *b*(*x*) (Lau & Lubensky 2007).

As we can see in the various Fokker–Planck equations above, if we fix *a*(*x*) and *b*(*x*), varying the parameter *α* gives different stochastic processes for the motion of the particle. This fact may make it seem like there should be a physically correct choice of the parameter *α*. We argue that this is false. For any fixed *α*, the range of stochastic processes that can be captured by an appropriate choice of *a*(*x*) and *b*(*x*) is the same. For example, suppose we fix a choice of *a*(*x*) and *b*(*x*) and choose to interpret (4.1) with a given *α*∈[0,1]. The process defined is identical to what we would obtain with the Itô interpretation (*α*=0) of the SDE with drift *a*(*x*)+*αb*′(*x*)*b*(*x*) and diffusion *b*(*x*).

Though the families of stochastic processes described using each convention are the same, it may still be the case that some choice of *α* is more natural or convenient for some purposes than others. Frequently, the rationale is based on the idea that if the drift is zero, then *X*(*t*) should have certain properties. For example, if one wants 〈*X*(*t*)−*X*(0)〉=0 for all *t* when *a*(*x*)≡0, regardless of *b*(*x*), then the Itô convention with *α*=0 guarantees this. If one wants that when *a*(*x*)≡0 the equilibrium density is constant, then the isothermal convention with *α*=1 guarantees this. We summarize the properties of the various interpretations of the SDE (4.1) when *a*(*x*)≡0 in table 2.

The question we posed in §1 may be rephrased as follows: in the absence of external forces, and given a diffusion , what is the correct choice of parameter *α* and drift *a*(*x*) to model the motion of the particle? A natural way to approach the problem is to interpret the absence of external forces as meaning that *a*(*x*)≡0. Then the problem boils down to the choice of *α*: the statistical-mechanics prediction follows from taking *α*=1, and the time-change prediction follows from taking *α*=0. The results in §2 showed that neither answer is justified universally.

In §3, we recommended a different approach. We fix *α* and then choose *a*(*x*) to generate the desired equilibrium distribution. As we have explained here, the choice of *α* is not crucial once we allow a non-zero *a*(*x*). Accordingly, we have chosen *α*=0, corresponding to Itô calculus. This is the main choice in the mathematics literature, and numerical methods such as the Euler–Maruyama method take a particularly simple form with it. Once we have made this choice of *α*, we are free to choose *a*(*x*) appropriately. In §3, we chose *a*(*x*) to ensure a given equilibrium density.

Beyond the particular needs of the present work, we believe the framework we describe in §3 provides a natural and flexible way to model diffusive systems. In situations where a researcher is confident for physical reasons that the equilibrium probability is constant, then our framework takes a simple form. The numerical method we present functions even when the diffusion and equilibrium density are discontinuous. Future work will study the convergence properties of the method, as well exploring applications of our general framework.

## Acknowledgements

The authors thank John Bechhoefer for suggesting the problem, and Florian Theil, Michael Plischke and Martin Zuckermann for helpful discussions. We also thank Nilima Nigam, John Bechhoefer and Nancy Forde for extensive comments on an earlier draft of this manuscript.

- Received April 27, 2012.
- Accepted August 6, 2012.

- This journal is © 2012 The Royal Society