Random walk in finite VERSUS infinite space: Probability density functions and their interpretation

Question

I am studying the probability density function of a random walk in a confined geometry (2D-BOX). I am also comparing this probability density function to its equivalent in infinite two-dimensional plane.

1) RANDOM WALK IN CONFINED GEOMETRY (2D-BOX)

The PDF for a random walk in a TWO-DIMENSIONAL BOX is obtained by multiplying the PDF for a random walk in a ONE-DIMENSIONAL BOX for the x- and y-directions. The PDF for a random walk in a ONE-DIMENSIONAL BOX is:

$$P(x,t|x_0,0)=\frac{1}{L}+\frac{2}{L}\sum_{i=1}^{\infty}\exp[-(\frac{i\pi}{2})^2\frac{t}{\tau}]\cos(\frac{ix\pi}{L})\cos(\frac{ix_0\pi}{L})$$

where $P(x,t|x_0,0)$ represents the probability that a particle move from $x_0$ to $x$ in time $t$, and $\tau = (\frac{L}{2})^2\frac{1}{D}$ , which represents the characteristic time-scale; $L$ is the length of the box (equal for each direction); $D$ is the diffusion coefficient.

2) RANDOM WALK IN INFINITE PLANE

In the case of a two-dimensional random walk in an infinite plane we can write:

$$P(p, t | p_0,0) = \frac{1}{4\pi Dt} \exp[-(p-p_0)^2/4Dt]$$

where $p = (x,y)$ and $p_0 = (x_0, y_0)$.

Now, my objective is to make sense of these PDFs. I implemented these functions in R (see code) and measured the PDF for time values between 1 and 10000, with time interval 100. I set D = 0.0015, L = 11, $p_0 = (2, 5)$ and $p = (8, 5)$. Since the summation to infinity cannot be execute, I chose a value of n_lim, such that the PDF does not change, which is equivalent to sum up to infinity.

prob1Dbox<-function(invL, t, invtau, x0, x, n_lim) {
  
  c =  pi * (pi/4) * (t * invtau)
  res = 0
  for(n in 1:n_lim){
    res = res + (exp(-1 * (n * n) * c) * cos((n * pi * x) * invL) * cos((n * pi * x0) * invL))
  }
  return(invL + ((2 * invL) * res))
}

prob2Dinfinite<-function(t, r, sigma) {
  a1 = 1/(2 * pi * sigma * t)
  b1 = exp(-(r * r)/(2.0 * sigma * t))
  return(a1 * b1)
}


prob2DBox<-function(t,  x0,  xt_pos,  y0,  yt_pos,  invLx,  invLy, invtau_x,  invtau_y, n_lim) {
  
  pbx = prob1Dbox(invLx, t, invtau_x, x0, xt_pos, n_lim)
  pby = prob1Dbox(invLy, t, invtau_y, y0, yt_pos, n_lim)
  return(pbx*pby)
}

If I plot the PDF of each model with respect to time, I get:

Random walk in confined geometry (2D-BOX)

Random walk in infinite plane

My concern is related to the "tail" of each distribution. In fact, in the infinite plane, we observe that the PDF decreases over time. This is because in an infinite plane a particle can indefinitely move. In the random walk in a 2D box, we observe that the PDF reaches a plateau over time. This seems to suggest that the probability to move from $p_0$ to $p$ remains constant after a certain time point. I am not sure whether it makes sense. Do you think that the PDF should be different? May be I wrongly implemented the PDF in my code.

UPDATE (29/06/2022)

Following the method of images suggested by @lpz, I implemented the proposed formula for the method of images in one-dimension as:

  prob1DboxImages<-function(t, D, x, x0_series){
  res = 0
  for(xi in x0_series){
    res = res + ((1/sqrt(4 * pi * D * t)) * exp(-((x-xi)^2)/(4 * D * t)))
  }
  return(res)
}

I then compared the PDF obtained from prob1DboxImages with the PDF given by prob1Dbox. Using the same parameter values as in the example above, setting Zlim = 1000, Zx = c(seq(0,Zlim,1), seq(0,Zlim,1) * (-1)) and x0_series = union((xA0 + (2 * Lx * Zx)),(-xA0 + (2 * Lx * Zx))) I do get:

Interestingly, I do not find differences between the two PDFs, for shorter time scales, at least for the parameter values here considered.

Answer 1

You are correct, the boundary conditions (von Neumann along the boundary of the box), the initial conditions (Dirac delta) and the equation of motion (diffusion equation) are separable. Therefore, the solution is separable as well. Expressed in $[0,L]^2$, after expanding out the factorization: $$ P(x,y,t|x_0,y_0) = \frac{1}{L^2}\\ +\frac{4}{L^2}\sum_{n_x,n_x\in\mathbb N^*{}^2}\cos(\pi n_xx/L) \cos(\pi n_xx_0/L)\\ \cos(\pi n_yy/L) \cos(\pi n_yy_0/L)\\ e^{-(n_x^2+n_y^2)t/\tau}\\ +\frac{2}{L^2}\sum_{n_x,n_x\in\mathbb N^*\times \{0\}\cup \{0\}\times N^*}\cos(\pi n_xx/L) \cos(\pi n_xx_0/L)\\ \cos(\pi n_yy/L) \cos(\pi n_yy_0/L)\\ e^{-(n_x^2+n_y^2)t/\tau}\\ $$ with $\tau =\frac{L^2}{D\pi^2}$, which is solution to the diffusion equation $\partial_t P=D\Delta P$ with von Neuman boundary conditions. This has the advantage of showing the decay rates of the various modes.

The difference between the two is that the box is compact while the plane is not. In the first case, the walk converges exponentially to the uniform distribution $f(x)=\frac{1}{L^2}$ hence the plateau effect. You can interpret this as ergodicity, the distribution will converge to the invariant measure. You can also view this spectrally. In this case, the spectrum is discrete. At long times, only the first constant term remains, and at next to leading contribution only the two degenerate modes $n=(0,1)$ and $n=(1,0)$ will contribute and give the rate of convergence $\tau$.

In the continuous case, there is no uniform distribution so it converges to $0$. Intuitively, this is because the walk is not recurrent and random walkers run off to infinity. Furthermore, the spectrum is not gapped, so the rate of convergence is slower (power law).

Your code looks correct, and log scales will help you better test out the calculations.

Hope this helps and tell me if you need more details.

--EDIT:

Sorry, my bad. You were right about the separability, all along, I had forgotten some modes. This time, I checked it numerically so it's correct. By the way, a good sanity check is to compare different formulas. While the above formula is good at long times, at short times, the one derived from the method of images is more accurate: $$ P(x,t|x_0) = \sum_{x'\in X}\frac{1}{\sqrt{4\pi Dt}}e^{-\frac{(x-x')^2}{4Dt}} $$ with $X = (x_0+2L\mathbb Z)\cup (-x_0+2L\mathbb Z)$.