[latexpage]
Basic notions
$df/dt$: 全微分，只随时间变化，没有空间坐标依赖关系 $\Rightarrow$ 描述的是随体运动坐标系中的物理量的性质（Lagrangian坐标系）。
转换到时间$t$-固定空间$x$坐标系（Eulerian坐标系）中描述，流体元$f(x,t)$的时空演化: $df/dt = \partial f/\partial t + \partial f/\partial x \cdot dx/dt = \partial f/\partial t + u_x \cdot \partial f/\partial x $. 其中，$\partial f / \partial t$是固定坐标系中的时间变化。
分布函数$f(x,v,t) \Rightarrow df/dt = \partial f/\partial t + u_x \cdot \partial f/\partial x + F/m \cdot \partial f/ \partial v $.
如果：1）$df/dt=0$,即无碰撞；2）$F=q/m(E+v\times B)$，即作用力为电磁力, 则$df/dt=0$这个方程叫做Vlasov Eq.
这里的记号$x,v$均应看作三维空间的矢量，每个矢量都有沿三个正交坐标方向的分量，如果采用笛卡尔坐标($e_x,e_y,e_z$)，则$x=(x,y,z)$和$v=(v_x,v_y,v_z)$. Vlasov codes就是在Eulerian坐标系上的($x,v$)6D相空间中用有限差分方法来数值求解Vlasov Eq.这个偏微分方程，需要计算和存储6D网格上每个格点的数据及其随时间的演化，计算量很大。
The meaning of distribution function
$f(x,v,t)$代表在时间$t$、位置$x$、速度为$v$的粒子数目。更严谨的说法是$fdxdv$代表$t$时刻在$x$和$x+dx$之间以及$v
$和$v+dv$之间的相空间体积元中的粒子数。它是“相空间”坐标($x,v$)和时间$t$的函数。如果默认考虑单位体积，即$f \equiv f/Vol$，$Vol$是整个体系的体积，则$f$代表粒子密度。如果$f$是归一化的，即$\int fdv=1$,则$f$代表几率。
对应通常密度$n=n(x,t)$来理解：如果对速度积分，即考虑在时间$t$和位置$x$包含了所有速度的(单位体积)粒子数目，则$\int_{-\infty}^{\infty}f(x,v,t)dv=n(x,t)$，这就是通常所说的密度$n(x,t)$，它是四个标量变量即位置$x$和时间$t$的函数，因为这个密度是通常坐标空间单位体积的粒子数，只随位置$x$和时间$t$发生变化。
本来完全描述一个N个粒子-体系的分布函数的形式应该是$F(x_1,v_1,\cdots,x_N,v_N,t)$，相空间是6N维。但是根据统计物理得出的结论，对于每个种类(电子或离子)的粒子来说，用一个单粒子分布函数就可以很好地描述这类粒子体系的性质，相空间只有6维，前提是粒子间的关联(碰撞)效应很弱，可以用微扰处理。
Liouville theorem
如前所述，全微分$df/dt$就解释为随体(相空间流体元)运动的坐标系中$f$随时间的变化率。那么为啥除了有碰撞外，相空间密度$f$不随时间变化，即有Liouville定理:$df/dt=0$ ?
我的理解是：因为在随体坐标系看来，$f$只是$t$的函数了，但是$f(t)$表示在时间$t$的粒子数目，如果没有碰撞或热运动，流体元中的所有粒子一起运动，流体元中的粒子数目不随时间变化；如果考虑单位体积，就是说相空间密度不随时间变化。如果有碰撞或热运动，这个随体的“体”就会逐渐瓦解，不能再代表一团流体元整体一起运动，流体元中的粒子数目$f$会随时间$t$发生变化，即此时$df/dt\neq 0$.
固定($x,v$)坐标系来看，相空间分布(即$f$)会随时间变化，这只是说$\partial f / \partial t \neq 0$. 而说无碰撞时分布函数不随时间变化，说的是$f$的全微分为零。
更严格深刻的理解，可以参考B. V. Somov的Plasma Astrophysics, Part I: Fundamentals and Particle一书，开篇就是Liouville定理，并给出了数学推导及守恒条件：
在相空间运动的粒子应遵守连续性方程，
\begin{equation}
\partial f / \partial t + \nabla \cdot (f\dot{X}) = 0,
\end{equation}
其中$X={x,v}$. 可以类比电荷连续性方程$\partial \rho / \partial t + \nabla \cdot j = 0$（其中$j=\rho v$）来理解：$f \rightarrow \rho$. 用${x,v}$将连续性方程的$X$具体写出：
\begin{equation}
\partial f / \partial t + \nabla_x \cdot (f \dot{x}) + \nabla_v \cdot (f \dot{v}) = 0.
\end{equation}
其中，左边第二项：
\begin{equation}
\nabla_x \cdot (f \dot{x}) = \nabla_x \cdot (f v) = f \nabla_x \cdot v + (v\cdot \nabla_x)f = 0 + (v\cdot \nabla_x)f
\end{equation}
上式最后那个0是因为在相空间中$x$和$v$是独立变量。再看上上式左边第三项：
\begin{equation}
\nabla_v \cdot (f \dot{v}) = f \nabla_v \cdot \dot{v} + \dot{v} \cdot \nabla_v f.
\end{equation}
所以连续性方程现在可写成
\begin{equation}
\partial f / \partial t + (v\cdot \nabla_x)f + f \nabla_v \cdot \dot{v} + \dot{v} \cdot \nabla_v f =0.
\end{equation}
到目前为止都还是严格的，没有任何假设和近似。
假如，现在开始我们只考虑上式左边第三项为零的情形，即：
\begin{equation}\label{fcondition}
f \nabla_v \cdot \dot{v} = 0.
\end{equation}
这样连续性方程变成
\begin{equation}
\partial f / \partial t + (v\cdot \nabla_x)f + \dot{v} \cdot \nabla_v f =0.
\end{equation}
这就是Liouville方程，而且$\partial / \partial t + v\cdot \nabla_x + \dot{v} \cdot \nabla_v $ 刚好就是全微分$d / dt$.
所以，相空间密度守恒即Liouville定理成立的关键就是满足条件(\ref{fcondition})，也就是说要求\textbf{某个方向的加速度或者力，与对应方向的速度无关}。可以验证，Lorentz力和重力满足这个条件，对相空间密度守恒来说是保守力；而摩擦力与碰撞力不满足，因为它们与速度有关，换言之，对相空间密度守恒来说是耗散力。\footnote{Mar. 27, 2017 Mon}
连续性方程是从物理出发得到的最根本性的方程，它不含任何假设。由它可以推导出Liouville和Vlasov方程。
Distribution function and particle-in-cell simulation
现在有了$df/dt=0$，即$f$在无碰撞时为常数，这意味着粒子在相空间中是沿等$f$线运动的，也就是说$f$的等值线就是粒子在相空间($x,v$)运动的轨道。如果考虑2D相空间($x,v_x$)，$f$的等值线是一条线，e.g.,所有粒子速度都为$v_0$，或产生一个沿$x$传播的波动，速度分布为一条正弦曲线$v_0\sin x$。考虑6D相空间，则$f$张成一个3D的超曲面，粒子轨道沿这个超曲面。随着系统演化，这个等值曲面可能发生拉伸、收缩、变形，但始终是3D的。用一定数量的Finite Phase-Fluid Elements (FPFE)，即一定数量的在相空间中的有限大小流体元（或称为宏粒子），来近似(approximate)或取样(sample)这个等$f$曲面，再跟踪它们的分布演化，就可得到整个粒子系统的演化信息，所以在相空间中的FPFE方法与通常说的PIC方法等价，不同的是前者中的宏粒子是相空间中的有限流体元，而后者中的宏粒子代表一大堆德拜球内的真实粒子。
采用FPFE方法的PIC的好处：
1）和Vlasov code比，大大减小了计算量，因为只用sample和follow有粒子分布或有感兴趣的物理过程发生的3D等$f$超曲面区域，而不用计算和存储在6D相空间中即使为零或不感兴趣的每个格点，从这个意义来说，这样的PIC codes可看作“压缩”或“Lagrangian”的Vlasov codes，而且用较少的宏粒子FPFE就能准确地取样超曲面;
2）很容易推广到one (macroparticle) to one (real particle) simulation的情形，只需把粒子间的碰撞力加上就行。
但如何保证这个方程是沿着3D等$f$超曲面演化的呢？
等$f$意味着没有加热(热运动造成相空间的diffusion)，因此保证系统总能量守恒就保证了等$f$的条件。系统总能量由粒子动能、静电能、电磁能、磁场能(Lorentz力自动保证能量守恒)组成。如果运动方程满足正则方程，则总能量即总Hamiltonian守恒。第$n$个FPFE(宏粒子)的正则运动方程为
\begin{equation}
dx_n/dt = p/m\gamma,
\end{equation}
\begin{equation}
dp_n/dt=F(+F_{collision}).
\end{equation}
再通过特定的电流分配方法，可保证电荷守恒，即保证满足连续性方程$\partial \rho / \partial t + \nabla \cdot j = 0$，那么电磁场变化的Maxwell方程组
\begin{equation}
\partial E / \partial t = c\nabla \times B ## 4\pi j,
\end{equation}
\begin{equation}
\partial B / \partial t = -c \nabla \times E
\end{equation}
也满足能量守恒。这样得到的系统总能量守恒。
current interpolation to the grid: 粒子有一定的中心位置，当粒子运动，沿运动路径积分得到电流密度，通过保证电荷守恒的方法，将电流密度分配到网格上。然后计算网格上的电场。（电荷密度是定义在网格中心的，通过形状因子将粒子位置分配到中心）
force interpolation to the actual particle position: 将网格上的电场分配到当前粒子位置，其分配方法和电流密度被分配到网格上的方法相同。通过正则运动方程由场推动粒子得到动量。\footnote{Mar. 21, 2017}
Why no collisions in PIC
描述N个粒子的分布函数$ F^N(\textbf{x}_1,\cdots,\textbf{x}_N,\textbf{p}_1,\cdots,\textbf{p}_N) $是6N维的。
\textbf{如果}满足粒子间碰撞(inter-particle correlations/collisions)效应很弱(可微扰处理)，则：
1）用单种粒子的分布函数$ f(\textbf{x},\textbf{p}) $来描述系统就足够了，只有6维；
2）粒子的运动满足Boltzmann-Vlasov Eq.: $ \partial f /\partial t + (P/\gamma m) \nabla f + (F/m) \nabla_p f = s.t. $，其中s.t.是粒子间碰撞项。
而sample the distribution function by a set of Finite Phase-Fluid Elements (FPFEs): $ f(\textbf{x},\textbf{p}) = \sum_n WEIGHT_{n}^{ph} SHAPE^{ph}(\textbf{x}-\textbf{x}n,\textbf{p}-\textbf{p}_n); \; d \textbf{x}_n / d t = \textbf{p}/ m \gamma; \; d\textbf{p}_n /d t = \textbf{F}{ext} + \textbf{F}_{St} $，这是PIC中采用的相对论粒子运动方程，它比Boltzmann-Vlasov Eq.还来得基本，因为Boltzmann-Vlasov Eq.的前提是粒子间关联或碰撞足够小，而FPFEs可以推广到处理包含了碰撞在内的one(numerical particle)-to-one(real particle) simulation！
为啥通常PIC里面没有包含碰撞呢？因为只考虑了外加电磁场的力$ \textbf{F}{external} $这一项，而忽略了$ \textbf{F}{St} $这一项，它是带电的粒子与粒子之间的静电库仑力，即inter-particle correlations/collisions.