\documentclass[journal=jpcbfk,manuscript=article,layout=singlecolumn,articletitle=true]{achemso}
\setkeys{acs}{articletitle = true}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Place any additional packages needed here.  Only include packages
%% which are essential, to avoid problems later.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{graphicx,bm,lineno,mathpazo}  
\usepackage{amssymb,amsfonts,amsmath,mathtools}
\usepackage{algorithm,algpseudocode,booktabs,dcolumn}
\captionsetup{labelfont=bf}
\usepackage[input-decimal-markers={.},separate-uncertainty=true,table-auto-round,scientific-notation = false,zero-decimal-to-integer=true]{siunitx}
\DeclareMathOperator*{\tr}{tr}
\DeclareMathOperator*{\argmax}{arg\,max}
\usepackage[caption=false,font=normalsize,position=top]{subfig}
\usepackage[kerning,spacing,tracking]{microtype}
\microtypecontext{spacing=nonfrench}
\usepackage{indentfirst}
\overfullrule=5pt
\usepackage[normalem]{ulem} 
\usepackage{color}
\usepackage{empheq}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% New command
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\boxedeq}[2]{\begin{empheq}[box={\fboxsep=6pt\fbox}]{align}\label{#1}#2\end{empheq}}

\title{The conservation of probability in Fokker-Planck Equation without dark operator}
\author{Yi-Tsao Chen}
\date{December 2019}

\begin{document}
 
\maketitle

\section{Fokker-Planck Equation(FPE)}
\begin{equation}
\label{eq:FPE}
\frac{\partial p(x,t)}{\partial t} = -D\frac{\partial}{\partial x}\left [ F(x)p(x,t) \right] + D \frac{\partial^2 p(x,t)}{\partial x^2}
\end{equation}
where $D$ is the diffusion coefficient, which is now assumed to be constant, and 
\begin{equation}
\label{eq:force}
F(x)=-\frac{d V(x)}{d x}
\end{equation}
for a given potential $V(x)$. Then, we imposed natural boundary conditions,
\begin{align}
\label{eq:BC1}
\lim_{x\to\pm\infty} p(x,t) &= 0
\end{align}
\begin{align}
\label{eq:BC2}
\lim_{x\to\pm\infty} \frac{\partial p(x,t)}{\partial x}&= 0
\end{align}

\section{Continuous equation and probability current}
According to eq (\ref{eq:FPE}),
\begin{equation}
\begin{split}
\label{eq:FPE2}
\frac{\partial p(x,t)}{\partial t}  &= -D\frac{d}{d x}\left [ F(x) p(x,t)\right ] +  D \frac{d^2 p(x,t)}{d x^2}\\
  &= D\frac{d}{d x}\left [ F(x) p(x,t)\right ] -  D \frac{d^2 p(x,t)}{d x^2}\\
  &= -D\frac{d}{d x}\left [ \frac{d V(x)}{d x}p(x,t) + \frac{d p(x,t)}{d x}   \right ] 
\end{split}
\end{equation}
Then,
\begin{equation}
\frac{\partial p(x,t)}{\partial t} + D\frac{d}{d x}\left [ \frac{d V(x)}{d x}p(x,t) + \frac{d p(x,t)}{d x}\right] = 0
\end{equation}
As a general continuity equation in a differential form:
\begin{equation}
\frac{\partial p}{\partial t} + \nabla \cdot \mathbf{j} = \sigma
\end{equation}
Here, $\sigma=0$, so
\begin{equation}
\label{eq:continutyeq}
\frac{\partial p}{\partial t} + \nabla \cdot \mathbf{j} = 0
\end{equation}
And $\mathbf{j}(x,t)$ is the current associated with the probability density $p$:
\begin{equation}
\label{eq:currentdefine}
\mathbf{j}(x,t) \equiv \frac{d V(x)}{d x}p(x,t) + \frac{d p(x,t)}{d x}
\end{equation}

\section{The conservation of probability}
Define $\mathcal{N}(t)$ as the integral of the probability density throughout space:
\begin{equation}
\mathcal{N}(t) \equiv \int p(x,t) dx
\end{equation}
We want to show the conservation of probability by the following equation
\begin{equation}
\label{eq:conserpb}
\frac{d \mathcal{N}(t)}{dt} = 0
\end{equation}
We begin by computing the time derivative of $\mathcal{N}(t)$ by using FPE
\begin{equation}
\begin{split}
\frac{d\mathcal{N}(t)}{dt} &= \int_{-\infty}^{\infty}\frac{\partial p(x,t)}{\partial t}dx \\
&= \int_{-\infty}^{\infty}-D\frac{\partial}{\partial x}\left [ \frac{d V(x)}{d x}p(x,t) + \frac{\partial p(x,t)}{\partial x}\right]dx \\
&= -D\int_{-\infty}^{\infty}\frac{\partial \mathbf{j}(x,t)}{\partial x}
\end{split}
\end{equation}
Here, we used eq (\ref{eq:FPE2}) and (\ref{eq:currentdefine}). Then, by natrual boundary conditions eq (\ref{eq:BC1}) and eq (\ref{eq:BC2}), we can get
\begin{equation}
\label{eq:bcconserve}
\begin{split}
\frac{d\mathcal{N}(t)}{dt} &= -D\left [ \mathbf{j}(\infty,t) - \mathbf{j}(-\infty,t)\right ] \\
&= -D\{[\frac{d V(\infty)}{d x}p(\infty,t) + \frac{\partial p(\infty,t)}{\partial x}]-[\frac{d V(-\infty)}{d x}p(-\infty,t) + \frac{\partial p(-\infty,t)}{\partial x}]\} \\
&= -D\{0-0\} = 0
\end{split}
\end{equation}
As you can see, the natural boundary conditions gurantee us the conservation of probability.

\section{The equilibrium probability density}
The equilibrium probability density,~$p_{eq}(x)$, is the long-time limit of $p(x,t)$ at $t\to\infty$ 
\begin{equation}
\label{eq:dpeqtconstant}
\frac{d p_{eq}(x)}{dt}=0
\end{equation}
Therefore, according to eq (\ref{eq:continutyeq})
\begin{equation}
\nabla \cdot \mathbf{j}_{eq} = 0
\end{equation}
That is,
\begin{equation}
\frac{d}{d x}\left [ \frac{d V(x)}{d x}p_{eq}(x) + \frac{d p_{eq}(x)}{d x}\right] = 0
\end{equation}
So
\begin{equation}
\frac{d V(x)}{d x}p_{eq}(x) + \frac{d p_{eq}(x)}{d x} = C
\end{equation}
Here, $C$ is a constant. And we have already known the value in two boundary as shown in eq. (\ref{eq:bcconserve}). That tells us the constant $C=0$. Therefore,
\begin{equation}
\label{eq:jzero}
\frac{d V(x)}{d x}p_{eq}(x) + \frac{d p_{eq}(x)}{d x} = 0
\end{equation}
Then,
\begin{align}
\frac{d p_{eq}(x)}{d x} &= -\frac{d V(x)}{d x}p_{eq}(x)\\
\frac{d p_{eq}(x)}{p_{eq}(x)} &= -d V(x) \\
\ln{(p_{eq}(x))} &= -V(x)
\end{align}
Finally, we have
\begin{equation}
\label{eq:defpeq}
p_{eq}(x) & = N^{-1} \exp{(-V(x))}= N^{-1}\mathcal{Y}(x)
\end{equation}
where we define
\begin{equation}
\label{eq:Yofx}
\mathcal{Y}(x)\equiv \exp{(-V(x))}
\end{equation}
and $N$ is a constant for normalization,
\begin{equation}
N = \int_{-\infty}^{\infty}\exp{(-V(x))} dx 
\end{equation}

\begin{figure}
\includegraphics{figures/peq_dotproduct}
\caption[eigenfunctions]{The left panel: For triple well potential, $\mathcal{Y}(x)= \exp{(-V(x))}$, and after normalizing, we get $p_{eq}(x)$, its quadrature is 1 with respect to weights $w(x)$. The right panel: the eigenvectors $|\psi_i(x)\rangle$ of $\bm{H}^0$ satisfies the completeness relationship.}
\end{figure}

\section{Transform FPE into an equation of Schr\"{o}dinger type}
Starting with the transformation $p(x,t) \rightarrow \rho(x,t)$
\begin{equation}
\label{eq:transform}
p(x,t) = \mathcal{Y}^{1/2}(x)\rho(x,t) \equiv \exp{(\frac{-V(x)}{2})}\rho(x,t)
\end{equation}
This transformation removes the first derivative in the original Fokker-Planck equation and generates the equation of Schr\"{o}dinger type(detail in Appendix A) for the function $\rho(x,t)$, i.e.,
\begin{align}
\label{eq:FPH0}
\frac{\partial \rho(x,t)}{\partial t} &= -\bm{H}^0 \rho(x,t) \\
\bm{H}^0 &= -D\nabla^2  + V_{\rm eff}(x) \\
V_{\rm eff}  &= \frac{DF'(x)}{2} + \frac{DF^2(x)}{4}
\end{align}

\section{Stationary states}
A stationary state is a quantum state with all observables independent of time. We can build stationary states if the potential energy function $V_{\rm eff}(x)$ is time-independent, and that is our case.
Recast our Schr\"{o}dinger-like FPE:
\begin{equation}
\label{eq:recastFPE}
\frac{\partial \Psi_i(x,t)}{\partial t} = -(-D\frac{\partial^2}{\partial x^2}+\frac{DF'(x)}{2} + \frac{DF^2(x)}{4})\Psi(x,t) = -\bm{H}^0 \Psi_i(x,t)
\end{equation} 
The signature property of a stationary state is that the position and the time dependence of the wavefunction factorize, which means that
\begin{equation}
\label{eq:factorize}
\Psi_i(x,t) = g_i(t)\psi_i(x)
\end{equation} 
Substitute eq(\ref{eq:factorize}) into eq(\ref{eq:recastFPE})
\begin{equation}
(\frac{d g_i(t)}{dt})\psi_i(x)=g_i(t) (-\bm{H}^0 \psi_i(x))
\end{equation} 
Then,
\begin{equation}
-\frac{1}{g_i(t)}\frac{d g_i(t)}{dt}=\frac{1}{\psi_i(x)}(\bm{H}^0 \psi_i(x))
\end{equation} 
The left-hand side is a function of $t$ only, while the right-hand side is a function of $x$ only because $\bm{H}^0$ contains no time dependent potential. The only way the two sides can equal each other for all values of $t$ and $x$ is for both sides to be equal to a constant $\lambda_i$ with units of $s^{-1}$
\begin{equation}
-\frac{1}{g_i(t)}\frac{d g_i(t)}{dt}=\frac{1}{\psi_i(x)}(\bm{H}^0 \psi_i(x)) = \lambda_i
\end{equation} 
Therefore 
\begin{equation}
g_i(t) = \exp{(-\lambda_i t)}
\end{equation}
and 
\begin{equation}
\label{eq:eigenvalue}
\bm{H}^0\psi_i(x) = \lambda_i \psi_i(x)
\end{equation}
This equation is the eigenvalue equation of $\bm{H}^0$.
Finally, we have
\begin{equation}
\label{Psi_last}
\Psi_i(x,t) = \exp{(-\lambda_i t)}\psi_i(x)
\end{equation} 

Furthermore, because $\bm{H}^0$ is a Hermitian operator(detail in Appendix B), it gives the orthonormal eigenvectors and form a complete set. According to spectral theorem, any solution $\rho(x,t)$ of eq(\ref{eq:recastFPE}) can be expressed as:
\begin{equation}
\begin{split}
\label{eq:rhoexpress}
\rho(x,t) &= c_1 \Psi_1(x,t) + c_2 \Psi_2(x,t) + \cdots + c_n \Psi_n(x,t) \\
&= c_1 \exp{(-\lambda_1 t)}\psi_1(x) + c_2 \exp{(-\lambda_2 t)}\psi_2(x) + \cdots + c_n \exp{(-\lambda_n t)}\psi_n(x) \\
&= \sum_{i=1}^{n} c_i \exp{(-\lambda_i t)}\psi_i(x)
\end{split}
\end{equation} 

\section{p(x,t)}
Given an initial probability distribution $p(x,0)$, and accroding to eq(\ref{eq:transform}) and eq(\ref{eq:rhoexpress}), we can determine the coefficients $c_i$
\begin{equation}
\begin{split}
p(x,t=0) &= \mathcal{Y}^{1/2}(x)\rho(x,t=0)\\
&= \mathcal{Y}^{1/2}(x)[c_1\psi_1(x)+c_2\psi_2(x)+\cdots+c_n \psi_n(x)]  
\end{split}
\end{equation}
After we have $c_i$, we can calculate the probability density at time $t$
\begin{equation}
\label{eq:pevolute}
\begin{split}
p(x,t) &= \mathcal{Y}^{1/2}(x)\rho(x,t)\\
&= \mathcal{Y}^{1/2}(x)[c_1 \exp{(-\lambda_1 t)}\psi_1(x)+c_2 \exp{(-\lambda_2 t)}\psi_2(x)+\cdots+c_n \exp{(-\lambda_n t)}\psi_n(x)]  
\end{split}
\end{equation}
Notice $c_1$ will be a constant(Appendix D and Appendix E) for arbitrary initial distribution if we fix the potential $V(x)$ (also $p_eq(x)$). To check whether the probability density is normalizable, we calculate
\begin{equation}
\mathcal{N}(t) \equiv \int p(x,t) dx
\end{equation}

\begin{figure}
\includegraphics{figures/onewell}
\caption[onewell]{The initial probability density corresponds to a one-well harmonic potential. This is the evolution of probability followed the equation \ref{eq:pevolute}.}
\end{figure}

\begin{figure}
\includegraphics{figures/twowell}
\caption[twowell]{The initial probability density corresponds to a double-well potential.}
\end{figure}


\clearpage
\section{Appendix A: Transform FPE into an equation of Schr\"{o}dinger type}
Insert eq (\ref{eq:transform}) into eq (\ref{eq:FPE})
\begin{equation}
\begin{split}
\frac{\partial p(x,t)}{\partial t} = -D\frac{d F(x)}{dx}p(x,t) - DF(x)\frac{\partial p(x,t)}{\partial x} + D \frac{\partial^2 p(x,t)}{\partial x^2}
\end{split}
\end{equation}
Then
\begin{equation}
\begin{split}
\exp{(\frac{-V(x)}{2})}\frac{\partial \rho(x,t)}{\partial t} &= D\frac{d^2 V(x)}{dx^2}\exp{(\frac{-V(x)}{2})}\rho(x,t) -\frac{1}{2}D(\frac{d V(x)}{dx})^{2}\exp{(\frac{-V(x)}{2})}\rho(x,t)\\
&+D\frac{dV(x)}{dx}\exp{(\frac{-V(x)}{2})} \frac{\partial \rho(x,t)}{\partial x} -\frac{1}{2}D\frac{d^2 V(x)}{dx^2}\exp{(\frac{-V(x)}{2})}\rho(x,t)\\
&+\frac{1}{4}D(\frac{d V(x)}{dx})^2\exp{(\frac{-V(x)}{2})}\rho(x,t) -D\frac{d V(x)}{dx}\exp{(\frac{-V(x)}{2})}\frac{\partial \rho(x,t)}{\partial x}\\ 
&+D\exp{(\frac{-V(x)}{2})} \frac{\partial^2 \rho(x,t)}{\partial x^2} \\ 
&= \frac{1}{2}D\frac{d^2 V(x)}{dx^2}\exp{(\frac{-V(x)}{2})}\rho(x,t) -\frac{1}{4}D(\frac{d V(x)}{dx})^{2}\exp{(\frac{-V(x)}{2})}\rho(x,t)\\
&+D\exp{(\frac{-V(x)}{2})} \frac{\partial^2 \rho(x,t)}{\partial x^2} 
\end{split}
\end{equation}
Multiply $\exp{(\frac{V(x)}{2})}$ on both sides,
\begin{equation}
\frac{\partial \rho(x,t)}{\partial t} = \frac{1}{2}D\frac{d^2 V(x)}{dx^2}\rho(x,t)-\frac{1}{4}D(\frac{d V(x)}{dx})^{2}\rho(x,t) +D\frac{\partial^2 \rho(x,t)}{\partial x^2}
\end{equation}
According to eq(\ref{eq:force}), 
\begin{align}
\frac{d^2 V(x)}{dx^2} &= -\frac{dF(x)}{dx}\\ 
(\frac{d V(x)}{dx})^{2} &= F^2(x)
\end{align}
So, we can get the last form of FPE
\begin{equation}
\frac{\partial \rho(x,t)}{\partial t} = -\frac{1}{2}D\frac{d F(x)}{dx}\rho(x,t)-\frac{1}{4} D F^{2}(x) \rho(x,t) +D\frac{\partial^2 \rho(x,t)}{\partial x^2}
\end{equation}

In the derivation, we used the following equations:
\begin{align}
\frac{\partial p(x,t)}{\partial t} &= \exp{(\frac{-V(x)}{2})}\frac{\partial \rho(x,t)}{\partial t}\\
\frac{\partial(\exp{(\frac{-V(x)}{2})})}{\partial x}&=-\frac{1}{2}\frac{d V(x)}{dx}\exp{(\frac{-V(x)}{2})}
\end{align}
\begin{equation}
\begin{split}
\frac{\partial p(x,t)}{\partial x} &= \frac{\partial \exp{(\frac{-V(x)}{2})}}{\partial x} \rho(x,t)+ \exp{(\frac{-V(x)}{2})} \frac{\partial \rho(x,t)}{\partial x}\\
&= -\frac{1}{2}\frac{d V(x)}{dx}\exp{(\frac{-V(x)}{2})}\rho(x,t)+\exp{(\frac{-V(x)}{2})} \frac{\partial \rho(x,t)}{\partial x}
\end{split}
\end{equation}
\begin{equation}
\begin{split}
\frac{\partial^2 p(x,t)}{\partial x^2} = &-\frac{1}{2}\frac{d^2 V(x)}{dx^2}\exp{(\frac{-V(x)}{2})}\rho(x,t) +  \frac{1}{4}(\frac{d V(x)}{dx})^2\exp{(\frac{-V(x)}{2})}\rho(x,t) \\ 
&-\frac{1}{2}\frac{d V(x)}{dx}\exp{(\frac{-V(x)}{2})}\frac{\partial \rho(x,t)}{\partial x} -\frac{1}{2}\frac{d V(x)}{dx}\exp{(\frac{-V(x)}{2})}\frac{\partial \rho(x,t)}{\partial x}\\
&+ \exp{(\frac{-V(x)}{2})} \frac{\partial^2 \rho(x,t)}{\partial x^2} \\
= &-\frac{1}{2}\frac{d^2 V(x)}{dx^2}\exp{(\frac{-V(x)}{2})}\rho(x,t) + \frac{1}{4}(\frac{d V(x)}{dx})^2\exp{(\frac{-V(x)}{2})}\rho(x,t) \\ 
&-\frac{d V(x)}{dx}\exp{(\frac{-V(x)}{2})}\frac{\partial \rho(x,t)}{\partial x} + \exp{(\frac{-V(x)}{2})} \frac{\partial^2 \rho(x,t)}{\partial x^2} 
\end{split}
\end{equation}

\clearpage
\section{Appendix B: $\bm{H}^0$ is a Hermitian operator}
If $\bm{H}^0$ is a Hermitian operator for functions $\psi(x)$ which satisfies boundary conditions like eq(\ref{eq:BC1}) and eq(\ref{eq:BC2}), the following equation should hold:
\begin{equation}
\langle \psi_i | \bm{H}^0 \psi_j \rangle = \langle \bm{H}^0 \psi_i | \psi_j \rangle
\end{equation}
or
\begin{equation}
\int_{-\infty}^{\infty}dx \psi^{*}_i \bm{H}^0 \psi_j = \int_{-\infty}^{\infty}dx (\bm{H}^0\psi_i)^* \psi_j
\end{equation}
Starting from the left hand side:
\begin{equation}
\begin{split}
\int_{-\infty}^{\infty}dx \psi^{*}_i \bm{H}^0 \psi_j &= \int_{-\infty}^{\infty}dx \psi^{*}_i [-D \frac{\partial^2 \psi_j}{\partial x^2} + \frac{1}{2}D F'(x)\psi_j + \frac{1}{4}DF^2(x)\psi_j] \\
&= -D\int_{-\infty}^{\infty}\psi^{*}_i \frac{\partial^2 \psi_j}{\partial x^2} dx +\frac{1}{2}D \int_{-\infty}^{\infty} F'(x) \psi^{*}_i \psi_j dx + \frac{1}{4}D\int_{-\infty}^{\infty} \psi^{*}_i F^2(x)\psi_j dx \\
& = D\int_{-\infty}^{\infty} \frac{\partial \psi^*_i}{\partial x} \frac{\partial \psi_j}{\partial x}dx +\frac{1}{2}D \int_{-\infty}^{\infty} F'(x) \psi^{*}_i \psi_j dx + \frac{1}{4}D\int_{-\infty}^{\infty} \psi^{*}_i F^2(x)\psi_j dx 
\end{split}
\end{equation}
The right hand side:
\begin{equation}
\begin{split}
\int_{-\infty}^{\infty}dx (\bm{H}^0\psi_i)^* \psi_j &= \int_{-\infty}^{\infty}dx [-D \frac{\partial^2 \psi_i}{\partial x^2} + \frac{1}{2}D F'(x)\psi_i + \frac{1}{4}DF^2(x)\psi_i]^{*}\psi_j \\
&= -D\int_{-\infty}^{\infty}\frac{\partial^2 \psi^*_i}{\partial x^2}\psi_j dx + \frac{1}{2}D\int_{-\infty}^{\infty}F'(x)\psi^*_i \psi_j dx + \frac{1}{4}D\int_{-\infty}^{\infty} F^2(x)\psi^*_i \psi_j dx \\
&= -D\int_{-\infty}^{\infty}\frac{\partial \psi^*_i}{\partial x}\frac{\partial \psi_j}{\partial x} dx + \frac{1}{2}D\int_{-\infty}^{\infty}F'(x)\psi^*_i \psi_j dx + \frac{1}{4}D\int_{-\infty}^{\infty} F^2(x)\psi^*_i \psi_j dx \\
\end{split}
\end{equation}
As you can see, the two sides result in the two functional form, so $\bm{H}^0$ is a Hermitian operator.
Here, we used integration by parts and boundary conditions:
\begin{equation}
\begin{split}
\int_{-\infty}^{\infty}\psi^{*}_i \frac{\partial^2 \psi_j}{\partial x^2} dx &= \int_{-\infty}^{\infty}\psi^{*}_i \frac{\partial}{\partial x}(\frac{\partial \psi_j}{\partial x}) dx \\
& = [\psi^{*}_i \frac{\partial \psi_j}{\partial x}]_{-\infty}^{\infty} - \int_{-\infty}^{\infty} \frac{\partial \psi^*_i}{\partial x} \frac{\partial \psi_j}{\partial x}dx\\
& = - \int_{-\infty}^{\infty} \frac{\partial \psi^*_i}{\partial x} \frac{\partial \psi_j}{\partial x}dx 
\end{split}
\end{equation}
\clearpage
\section{Appendix C: When the initial distribution is $\delta(x-x_0)$}
According to eq (\ref{eq:transform}),
\begin{equation}
\label{eq:C1}
p(x,t=0) = \mathcal{Y}^{1/2}(x)\rho(x,t=0) = \delta(x-x_0)
\end{equation}
and further substitute eq (\ref{eq:rhoexpress}) into eq (\ref{eq:C1}),
\begin{equation}
\label{eq:C2}
\mathcal{Y}^{-1/2}(x) \delta(x-x_0) = \sum_{i=1}^{n} c_i \psi_i(x)
\end{equation}
Then we multiply both sides of eq (\ref{eq:C2}) by $\psi_j$
\begin{equation}
\label{eq:C3}
\mathcal{Y}^{-1/2}(x) \psi_j(x) \delta(x-x_0) = \sum_{i=1}^{n} c_i \psi_j(x) \psi_i(x)
\end{equation}
And integrate eq (\ref{eq:C3}) over space
\begin{equation}
\label{eq:C4}
\int_{-\infty}^{\infty} \mathcal{Y}^{-1/2}(x) \psi_j(x) \delta(x-x_0)dx = \int_{-\infty}^{\infty} \sum_{i=1}^{n} c_i \psi_j(x) \psi_i(x) dx
\end{equation}
At last, we can get the coefficients $c_i$
\begin{equation}
\label{eq:C5}
c_i = \mathcal{Y}^{-1/2}(x_0) \psi_i(x_0)
\end{equation}
\clearpage
\section{Appendix D: Solve $\psi_1(x)$}
For $p_{eq}(x)$, there is a corresponding $\rho_{eq}(x)$. 
\begin{equation}
p_{eq}(x) = \exp{(\frac{-V(x)}{2})}\rho_{eq}(x)
\end{equation}
According to eq(\ref{eq:jzero})
\begin{equation}
\begin{split}
&\frac{d V(x)}{d x}p_{eq}(x) + \frac{d p_{eq}(x)}{d x}\\
&= \frac{d V(x)}{d x} \exp{(\frac{-V(x)}{2})}\rho_{eq}(x) + \frac{d}{dx}(\exp{(\frac{-V(x)}{2})}\rho_{eq}(x))\\
&= \frac{d V(x)}{d x} \exp{(\frac{-V(x)}{2})}\rho_{eq}(x) - \frac{1}{2}\frac{d V(x)}{d x} \exp{(\frac{-V(x)}{2})}\rho_{eq}(x) + \exp{(\frac{-V(x)}{2})} \frac{d\rho_{eq}(x)}{dx}= 0
\end{split}
\end{equation}
Then, multiply $\exp{(\frac{V(x)}{2})}$ on the both sides
\begin{equation}
\begin{split}
\frac{d\rho_{eq}(x)}{dx} = -\frac{1}{2}\frac{d V(x)}{d x}\rho_{eq}(x) 
\end{split}
\end{equation}
Therefore, we obtain
\begin{equation}
\rho_{eq}(x) = C \exp{(-\frac{V(x)}{2})}
\end{equation}
where C is an arbitrary constant.

On the other hand, the equilibrium probability is steady
\begin{equation}
\frac{\partial p_{eq}(x)}{\partial t} = \frac{\partial }{\partial t}(\exp{(\frac{-V(x)}{2})}\rho_{eq}(x)) = 0
\end{equation}
so
\begin{equation}
\frac{\partial \rho_{eq}(x)}{\partial t} = 0
\end{equation}
And according to eq(\ref{eq:recastFPE}) 
\begin{equation}
\label{eq:D1}
\frac{\partial \rho_{eq}(x)}{\partial t}  = -\bm{H}^0 \rho_{eq}(x) = 0
\end{equation} 
This means that there is a eigenvalue is zero, $\lambda_1=0$, and the corresponding eigenvector $\psi_1(x)$
\begin{equation}
\label{eq:D2}
\bm{H}^0 \psi_1(x) = \lambda_1 \psi_1(x) = 0 
\end{equation} 
At last, by combining eq(\ref{eq:D1}) and eq(\ref{eq:D2}), we get
\begin{equation}
\psi_1(x) = \rho_{eq}(x) = C \exp{(-\frac{V(x)}{2})}
\end{equation} 
Here, $C=\sqrt{\frac{1}{\int_{-\infty}^{\infty}\exp{(-V(x))}dx}}=N^{-1/2}$, because
\begin{equation}
\int_{-\infty}^{\infty}|\psi_1(x)|^2 dx = 1
\end{equation} 
As a result
\begin{equation}
\label{eq:D3}
\psi_1(x)=\rho_{eq}(x) = \sqrt{p_{eq}(x)} 
\end{equation} 
\clearpage
\section{Appendix E: $c_1$ is a constant for all initial distributions}
When $t >> 1$, according to eq(\ref{eq:pevolute}) and also $\lambda_1=0$
\begin{equation}
\label{eq:E0}
p(x,t) = p_{eq}(x) = c_1 \mathcal{Y}^{1/2}(x)\psi_1(x)
\end{equation}
Substitute eq(\ref{eq:D3}) into eq(\ref{eq:E0})
\begin{equation}
\label{eq:E1}
p_{eq}(x) = c_1 \exp{(-\frac{V(x)}{2})} N^{-1/2}\exp{(-\frac{V(x)}{2})} = c_1 N^{-1/2} \exp{(-V(x))}
\end{equation}
Then, eq(\ref{eq:defpeq}) gives us  
\begin{equation}
\label{eq:E2}
p_{eq}(x) = N^{-1}\mathcal{Y}(x) = N^{-1}\exp{(-V(x))}
\end{equation}
Compare eq(\ref{eq:E1}) and eq(\ref{eq:E2})
\begin{equation}
c_1 = N^{-1/2} = \sqrt{\frac{1}{\int_{-\infty}^{\infty}\exp{(-V(x))}dx}}
\end{equation}
As you can see, $c_1$ is totally determined by $V(x)$, so any initial distributions won't affect it.
\begin{figure}
\includegraphics{figures/c1_constant}
\caption[c1]{$c_1$ is a constant for any initial distribution.}
\end{figure}

\end{document}
