Brownian motion in a harmonic potential is described by

\(m\ddot x + c\dot x + kx = F_r(t)\)

where \(F_r(t)\) is random force, which we represent as white noise:

\(<F_r(t)> = 0, <F_r(t)F_r(t+\tau)> = 2ck_BT\delta(\tau)\)

where \(k_B\) is the Boltzmann constant, and \(T\) is the temperature.

Initialization

time = c(0.,0.001,20001.) #min,step,n
initial = c(100.,0.) #position, velocity
parameter = c(1.,10.,100.,100.) #m,c,k,kT
iter <- 5
max_tau <- 1000
set.seed(1234)

Functions

Time_Grid<-function(time){
     n = as.integer(time[3])
     min = time[1]
     step = time[2]
     x <- seq(from=min,by=step,length=n)
     x
}

Langevin<-function(time,initial,parameter,time_grid){
     n <- time[3]
     step <- time[2]
     c <- parameter[2]
     kT <- parameter[4]
     position <- vector(mode = "double",length = n)
     velocity <- vector(mode = "double",length = n)
     position[1] <- initial[1]
     velocity[1] <- initial[2]
     for(i in 2:n){
          force <- rnorm(n = 1,mean = 0,sd = (2*c*kT)^0.5)
          acc <- Solve(parameter,position[i-1],velocity[i-1],force)
          position[i] <- 0.5*acc*step*step + velocity[i-1]*step + position[i-1]
          velocity[i] <- velocity[i-1] + acc*step
     }
     output <- data.frame(position,velocity)
}

Solve<-function(parameter,position,velocity,force){
     m <- parameter[1]
     c <- parameter[2]
     k <- parameter[3]
     acc <- (force - k*position - c*velocity)/m
     acc
}

Probability Density Functions

First, we show the time-position and time-velocity graphs starting from the initial conditions and continuously progressing for specified iterations. At the first iteration, the system has not yet been in the stationary state. The stationary state starts since the second iteration. The probability density of the stationary state seems to be Gaussian with center at zero for both position and velocity.

time_grid <- Time_Grid(time)
for(i in 1:iter){
     output <- Langevin(time,initial,parameter,time_grid)
     position <- output$position
     velocity <- output$velocity
     par(mfrow=c(2,2))
     plot(time_grid,position)
     plot(time_grid,velocity)
     hist(position,breaks=12)
     hist(velocity,breaks=12)
     cat(sprintf("Iter=%i\n",i))
     cat(sprintf("Position: mean=%g, sd=%g\n",mean(position),sd(position)))
     cat(sprintf("Velocity: mean=%g, sd=%g\n",mean(velocity),sd(velocity)))
     initial <- c(position[time[3]],velocity[time[3]])
}

## Iter=1
## Position: mean=0.499725, sd=7.05913
## Velocity: mean=-4.99989, sd=49.8633

## Iter=2
## Position: mean=0.00307897, sd=0.030466
## Velocity: mean=0.00164598, sd=0.30539

## Iter=3
## Position: mean=-0.00226191, sd=0.0335634
## Velocity: mean=-0.00145065, sd=0.328978

## Iter=4
## Position: mean=0.000536344, sd=0.0307944
## Velocity: mean=-0.0010217, sd=0.305495

## Iter=5
## Position: mean=0.00526463, sd=0.0297039
## Velocity: mean=0.00305947, sd=0.303701

Theoretically, \(SD(x;T) = sqrt(k_BT/k), SD(v;T) = sqrt(k_BT/m)\) which equals

m <- parameter[1]
c <- parameter[2]
k <- parameter[3]
kT <- parameter[4]
sd_position <- (kT/k)^0.5
sd_velocity <- (kT/m)^0.5
cat(sprintf("SD(x;T) = %g, SD(v;T) = %g",sd_position,sd_velocity))
## SD(x;T) = 1, SD(v;T) = 10

Autocorrelation Functions

n <- length(position)
Rxx <- c()
for(i in 1:max_tau){
     sum <- 0.
     for(j in 1:(n-i)){
          sum <- sum + position[j]*position[j+i]
     }
     sum <- sum/j
     Rxx <- c(Rxx,sum)
}
tau <- seq(from=time[2],by=time[2],length=max_tau)
plot(x = tau,y = Rxx,type="l",main="Aucorrelation function of position: Rxx")

n <- length(position)
Rvv <- c()
for(i in 1:max_tau){
     sum <- 0.
     for(j in 1:(n-i)){
          sum <- sum + velocity[j]*velocity[j+i]
     }
     sum <- sum/j
     Rvv <- c(Rvv,sum)
}
tau <- seq(from=time[2],by=time[2],length=max_tau)
plot(x = tau,y = Rvv,type="l",main="Aucorrelation function of velocity: Rvv")

Theoretically,

\(R_{xx}(\tau) = \frac{ck_BT}{\pi m^2} \int d\omega \frac{e^{i\omega\tau}}{(\omega_0^2 - \omega^2)^2+(\omega\beta)^2}\)

where \(\omega_0 = sqrt(k/m)\), and \(\beta = c/m\). For underdamped oscillation (i.e. \(\omega_1^2 = \omega_0^2 - (\beta/2)^2 > 0\), and \(\tau >= 0\)), we get

\(R_{xx}(\tau) = \frac{k_BT}{k}e^{-\frac{c\tau}{2m}}(\cos(\omega_1\tau) + \frac{\beta}{2\omega_1}\sin(\omega_1\tau))\)

and

\(R_{vv}(\tau) = -\frac{d^2R_{xx}(\tau)}{d\tau^2}\).

A <- exp(-c*tau/(2*m))
B <- (k/m) - (c/(2*m))^2
C <- cos(B*tau) + (c/(2*m*B))*sin(B*tau)
Rxx_theory <- (kT/k)*A*C
plot(x = tau,y = Rxx_theory,type="l",main="Aucorrelation function of position: Rxx_theory")

n <- length(Rxx_theory)
step <- tau[2] -  tau[1]
Rxv_theory <- c()
for(i in 1:(n-1)){
     xx <- (Rxx_theory[i+1] - Rxx_theory[i])/step
     Rxv_theory <- c(Rxv_theory,xx)
}
n <- length(Rxv_theory)
Rvv_theory <- c()
for(i in 1:(n-1)){
     xx <- -1*(Rxv_theory[i+1] - Rxv_theory[i])/step
     Rvv_theory <- c(Rvv_theory,xx)
}
plot(x = tau[1:i],y = Rvv_theory,type="l",main="Aucorrelation function of velocity: Rvv_theory")