时间之矢

本文是《Statistical and Thermal Physics Using Python》¹一书的学习笔记。

一个和多个粒子

实心球的可逆弹跳

球撞击弹性板，其受到两个势的影响：重力势能\(U_G=mgz\)和模拟与板相互作用的线性弹簧\(U_F = -(1/2)kz^2\) (当\(z<0\))及0 (当\(Z>0\))。 该过程的曲线(\(t\)-\(z\))在时间上是完全对称的：令\(t\to -t\)，曲线看起来一样。你看不出我们是在后退还是在前进方向上打球。这是因为单个粒子的运动在时间上是可逆的。 该点可从牛顿第二定律看出： \[\begin{equation} m \frac{d^2 \mathbf{r}}{d t^2}=-\nabla U(\mathbf{r})=m \frac{d^2 \mathbf{r}}{d(-t)^2} \end{equation}\] 这也能从能量守恒得出： \[\begin{equation} \frac{1}{2} m v^2=U(\mathbf{r}) \Rightarrow v= \pm \sqrt{2 U(\mathbf{r}) / m} \end{equation}\] 其中，有两个可能的符号对应着两个可能的运动方向。

对于单个粒子(受保守势力影响)的运动在时间上是可逆的， 单个粒子的运动将永远是保守和可逆的。

弹性球的不可逆弹跳

如果将上述实心球换成弹性球，则球不会反弹到原来的高度，而是在较低的高度上反弹。

系统中的总能量没有发生变化。运动是保守的，能量的变化可以忽略不计。然而，我们可以把能量分为两部分，质心的能量和相对质心运动形式的内能以及弹簧内部压缩或伸长的内能。

质心的能量为\(E_{cm}\)，则内能\(\Delta E = E_{TOT}-E_{cm}\)，且是时间的函数。 我们看到，最初所有的能量都在质心运动中，但是经过短暂的时间后，能量就分布在质心运动和内能之间。

显然，最初状态是特殊点。系统似乎不太可能移动到这个特殊状态。注意，是不大可能，不是不可能。 时间之矢的核心是：在时间中向前运动并不是唯一的可能性，而是最可能的可能性。 而且，我们可以看出系统朝向一种状态发展 where 能量分布在系统中所有可能的自由度上。并且这种状态对应一种“平衡” (equilibrium)。 在多粒子系统中，能量有多种分配方式，较可能的状态对应于能量的均匀分布(本段理解需要结合原书Fig. 2.1)。

我们现在进入了许多粒子的新物理定律。它与可能性有关——与概率有关——因此它们不同于牛顿定律，牛顿定律是绝对定律。 我们将从一个简单的原子系统——原子气体出发，来讨论它们之间的差异。

我们将一维弹性变形球建模为一组质量\(m\)，由弹簧常数\(k\)的弹簧连接。弹簧从下往上依次枚举\(i = 1，..，N\)。最底层的质量也与地板相互作用。则颗粒\(i\)受到的作用力为 \[\begin{equation} m_i a_i=F_i=-m_i g-k\left(z_{i+1}-z_i-b\right)+k\left(z_i-z_{i-1} -b\right)-F_w \end{equation}\] 其中，当\(i=1\)且\(z_i < 0\)时，\(F_w = -k_w z_i\)。其余地方均为零。

import numpy as np
import math
from numba import jit
import matplotlib.pyplot as plt

N = 10 # Nr of elements 
k = 100000
kw = 1000000 # Spring constants 
d = 0.1 # Size in m 
m = 0.1 # Mass in kg 
t = 20 # Simulation time 
g = 9.81 # Acc of gravity 
z0 = 0.1 # Initial height 
dt = 0.00001 # Timestep 
nstep = math.ceil(t/dt) #向上取整

z = np.zeros((N,nstep),float)
vz = np.zeros((N,nstep),float)
timex = np.zeros(nstep)
z[:,0] = d*(np.array(range(N)))/N+z0 # Initial positions 

l = d/N # Distance between nodes 

@jit
def diff(arr):
    diff_arr = np.empty(len(arr)-1, dtype=arr.dtype)
    for i in range(len(arr) - 1):
        diff_arr[i] = arr[i+1] - arr[i]
    return diff_arr

#@jit 加速存在问题，因此该函数暂不加速
def calculate():
    for i in range(nstep-1): # Calculate motion 
        dz = diff(z[:,i])-l 
        F = -k*np.append(0.0,dz) + k*np.append(dz,0.0) - m*g # 内部力, 添加0对应最上和最下两个球的受力情况，
        F[0] = F[0] - kw*z[0,i]*(z[0,i]<0) # Bottom wall
        a = F/m
        vz[:,i+1] = vz[:,i]  + a*dt
        z[:,i+1] = z[:,i] + vz[:,i+1]*dt
        timex[i+1] = i*dt
        if i%50000==0:
            print(i)

calculate()
plt.plot(timex, z[9,:], '-')

更多的自由度

能量守恒并不是对宏观行为重要的唯一定律。

能量守恒不足以描述时间之箭。

达到平衡——分子动力学

原子假说本身与牛顿力学的结合以及对原子间相互作用的简化描述为宏观系统的行为提供了有价值的见解和良好的预测。

可以利用量子力学来寻找多原子体系的有效势能。此处，假定该能量仅取决于原子间的距离。

两个惰性气体原子间的相互作用主要有两个贡献：

• 由于偶极-偶极相互作用，存在吸引力。其正比于\((1/r)^6\)，其中\(r\)为原子间距。(范德瓦尔斯相互作用，范德瓦尔斯力);
• 当两个原子被推挤在一起时，由于电子轨道重叠的可能性，存在一个排斥力，这是一种QM效应。我们用一个形式为\((1/r)^n\)的幂律来表示这种相互作用。通常选择\(n = 12\)，这对氩的行为给出了很好的近似。

该组合模型称为 Lennard-Jones potential: \[\begin{equation} U(r)=4 \epsilon\left(\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^6\right) \end{equation}\]

为了求解运动方程，需要使用LAMMPS。

达到平衡

通过在平衡之外启动系统的方法来详细研究系统如何达到平衡。

将一个盒子分成两部分，用原子气体填充一部分，另一部分为空。原子以随机速度开始，但总动量为零。从\(t=0\)开始，我们允许系统在时间上发展。

lmp_serial < in.gastwosection01.lmp

在终端中输入上述命令以上述

# 2d Lennard-Jones gas
units lj
dimension 2
atom_style atomic
lattice hex 0.25
region mybox block 0 20 0 10 -0.5 0.5
create_box 1 mybox
region 2 block 0 10 0 10 -0.5 0.05
create_atoms 1 region 2
mass 1 1.0
velocity all create 0.5 87287
pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0 2.5
neighbor 0.3 bin
neigh_modify every 20 delay 0 check no
fix 1 all nve
dump 1 all custom 10 twosec01.lammpstrj id type x y z vx vy vz
restart 2000 mydata.restart01 mydata.restart011
thermo 100
run 2000

如果令时间倒流，从2000时间步开始运行模拟，会发生什么？

# 2d Lennard-Jones gas
read_restart mydata.restart01
# Reverse all velocities
variable vx atom -vx
variable vy atom -vy
velocity all set v_vx v_vy NULL
fix 1 all nve
dump 1 all custom 10 twosec11.lammpstrj id type x y z vx vy vz
restart 2000 mydata.restart11 mydata.restart111
thermo 100 # Output thermodyn variables every 100 timesteps
run 2000 # Number of timesteps

粒子几乎回到了开始的地方，该行为是不物理的：我们不会期望原子突然组织起来。相反，我们期望原子继续很好的混合。 事实上，如果我们将第一个模拟继续进行，粒子在随机移动中保持混合。从某种意义上说，这种混合状态似乎比开始的状态更有可能。

可以通过测量系统左半边和右半边粒子数。若用\(n(t)\)表示左边的原子数，由于原子守恒，则右边原子个数为\(N-n(t)\)。

原文是使用官方的pizza.py读取不同时间步下原子的分布，然而不幸的是pizza.py是使用python2写的，不兼容目前较为广泛使用的python3。 为了实现该目标，我们需要自己了解数据结构并建立读取的函数。

filename = './twosec21.lammpstrj'

## 给定计算原子个数的范围
ranges = [[0, 25.31],
        [0, 3.7224194364083985e+01], 
        [-1.0745699318235420e+00, 1.0745699318235420e+00]]

## 获得所有时间步
timesteps = getAllTimesteps(filename)
results = np.zeros(len(timesteps))

for t in range(len(timesteps)):
    results[t] = calculateTheAtomsNumbers(filename, timesteps[t], ranges)

plt.plot(timesteps, results)

上图显示了左边的粒子数收敛到近似恒定的值，然后围绕这个值上下波动。可以看出系统明显朝特定方向运动。

这一时间箭头对于我们理解宏观系统是至关重要的。而且它显然是一个新的定律- -除了动量守恒或能量守恒之外的定律。

达到平衡--算法模型

本节将建立仅包含过程最关键特征的算法模型。

左侧原子的模型

分析上述过程的核心特征：

• 封闭盒子；
• 最初一侧有\(N\)个原子，另一侧没有；
• 随着原子向四周移动，原子发生混合；
• 一段时间后原子的分布变得更加均匀。

通过在系统的中间引入实际的有孔的墙来使得系统更加理想化。 这会减慢平衡过程，因为原子需要一段时间才能命中孔。 我们应该简化这个过程，并引入一个算法模型，该模型只捕获过程中最重要的部分- -原子从一个盒子到另一个盒子的运动

简化系统--第一部分

需要做一些简化假设：

• 假设所有的原子都是独立的；
• 对于每一个原子，在单位时间内都有一定的概率穿过孔洞向另一侧移动。由于原子是独立同构的，所以这个概率对于所有的原子来说一定是相同的。

现在我们已经将问题转化为一个基于规则的模型，我们可以用它来发现系统的行为：

• 该系统由N个独立原子组成，这些原子被带有小孔的壁隔开。每个原子的位置由\(s_i\)给出，因此若原子原子\(i\)在左框中，则\(s_i = 1\)；若在右框，则\(s_i = 0\)。

• 原子在时间间隔\(\Delta t\)内从盒子一侧穿过孔到另一侧的概率为\(p=R\Delta t\)

• 左侧有\(n(t) = \Sigma_i s_i\)个原子，右侧有\(N-n(t)\)个

这意味着可以为系统随时间的发展建立简单的原子。在时间\(t\)系统的状态由所有N个原子的位置\(s_i(t)\)所决定。在\(t\)到\(t+\Delta t\)的\(\Delta t\)时间间隔内，原子\(i\)有\(p\)的概率从一侧到另一侧。也就是每一个原子都有\(p\)的可能改变它的位置状态。

简化系统--第二部分

上述模型仍然比必要的复杂得多，为了进一步简化，我们以原子通过孔的典型时间\(\Delta t\)来划分时间，则规则是：

• 在每个间隔时间，随机选择一个原子将其移动到另一侧；

这是一个可能的模型，但它仍然有点过于复杂。我们可以使模型更简单，更切中要害，因为我们对每个原子的位置不感兴趣 - 我们只对每边有多少原子感兴趣。 如果我们随机选择一个原子并将其移动到另一侧，则我们在左侧选择原子的概率为\(n/N\)。在右侧选择一个原子的概率是\((N − n)/N\)。 因此，我们将原子从左侧移动到右侧的概率为\(n/N\)。因此，我们得出了一个简单的算法:

• 生成一个介于 0 和 1 之间的随机数 \(r\)；

• 如果 \(r ≤ n/N\)，则从左向右移动一个原子 ，否则我们将一个原子从右向左移动；

• 将时间增加\(\Delta t\)

使用该算法来研究趋于平衡的方法：

import numpy as np
import matplotlib.pyplot as plt

N = 210 # 粒子数量
nstep = 4000 #步数
n = np.zeros(nstep)
n[0] = 210 # 左侧初始值

for i in range(1,nstep):
    r = np.random.random()
    if (r<n[i-1]/N):
        n[i] = n[i-1] - 1 # Move atom from left to right
    else: 
        n[i] = n[i-1] + 1 # Move atom from right to left

plt.plot(range(0,nstep),n)
plt.xlabel('t')
plt.ylabel('n')
plt.show()

涨落和原子数

这个简化的模型使得我们可以有效地研究达到平衡态和平衡态的涨落，例如通过解决对原子数N的依赖关系。 当我们增加原子数量时会发生什么？ 基于我们的讨论，我们期望系统变得更加不可逆。

import numpy as np
import matplotlib.pyplot as plt

N = np.array([20,40,60,80]) # 粒子数量
nstep = 4000 #步数

for j in range(0,np.size(N)):
    NN = N[j]
    n = np.zeros(nstep)
    n[0] = NN # 左侧初始值
    for i in range(1,nstep):
        r = np.random.random()
        if (r<n[i-1]/NN):
            n[i] = n[i-1] - 1 # Move atom from left to right
        else: 
            n[i] = n[i-1] + 1 # Move atom from right to left
    label = 'N = '+str(NN)
    plt.plot(range(0,nstep),n,label=label)
    plt.xlabel('t')
    plt.ylabel('n')

plt.xlim(0, 400)
plt.legend()
plt.show()

可以看到，随着N的增大，波动变小。这意味着系统在较长时间内保持在与平均n相对应的n的平衡值附近。当粒子数N变大时，发生大漂移的概率变小。这也意味着当粒子数很小时，完全有可能在所有粒子都在左侧(\(n = 1\))的地方产生波动，但是随着\(N\)变大，这变得越来越不可能。这说明了我们在这里探讨的物理规律的原理：一个远离平均数的行为，随着原子数的增加，变得越来越不可能--不是不可能，只是不大可能。

因此，宏观系统的物理规律是典型宏观系统中原子数目众多的结果。数量的庞大，保证了物理规律的规律性。 数量非常大的时候涨落变得非常非常小。

\(n\)的平均值

通过求\(n\)在一定时间步内的平均值：

print(np.average(n[1000:]))

结果为0.497384126984127。 非常接近我们的期望值\(1/2\)。

初始条件的依赖性

经过一定时间后系统的行为将不依赖初始状态，这意味着当系统达到平衡状态时，行为不再依赖于初始条件- -行为反而是一般的，取决于系统处于平衡状态时的性质。

达到平衡--理论模型

总结

• 在原子层次上，所有定律在时间上都是可逆的。
• 在宏观层次上，时间具有特定的方向。
• 能量守恒不是一个充分的原则。我们需要新的原理来描述能量的分布。
• 我们可以将时间的方向理解为朝向系统中最有可能出现的状态的发展。

练习

暂略，复习的时候再做。

---

1. Malthe-Sørenssen, A. & Dysthe, D. K. Statistical and Thermal Physics Using Python.↩︎