I'm currently studying this article because I'm interested in the technique for finding the Green's function $G$ of the following PDE: $$\frac{\partial B}{\partial t}+\frac{\sigma^2}{2}\frac{\partial^2 B}{\partial x^2}+(a-bx)\frac{\partial B}{\partial x}-xB=0 \label{1} \tag{1}$$ This equation comes from the bond princing of the Vasicek model, which is a stocastic differential process, like the CIR model of the article, following $dX_t=(a-bX_t)dt+\sigma dW_t$ with $W_t$ Brownian motion and $a,b,\sigma>0$. The solution to this equation can be defined also as the expected value $B(x,t)=\mathbb{E}(e^{\int_s^tX_\tau d\tau}\mid \mathcal{F_s})$ under the natural filtration for Brownian motion and using Feynman-Kac theorem we get our equation. Since $t\in[0,T]$ I have the classic final condition for ZCB: $B(x,T)=1$.
Anyway I'm not understanding some steps used in the article: at the beginning they define the probability density using the expected value as $f(x,y,\tau)=\mathbb{E}(\delta(X_t-y)\mid X_s=x)$, with $\tau=t-s$, because they start from a stochastic differential equation of type $dX_t=(a-bX_t)dt+\sigma\sqrt(X_t)dW_t$, with $W_t$ Brownian motion and $a,b,\sigma>0$ (the CIR model). Then they use Feynman-Kac theorem to obtain the corrisponding partial differential equation $$\partial_tf(x,y,t-s)=-Hf(x,y,t-s) \tag{*}$$ in their case it's the Kolmogorov backward equation so their differential operator is $H=(a-bx)\partial_x+\frac{\sigma^2}{2}x\partial^2_x$. In my case I have a different differential operator $\tilde{H}=-x+(a-bx)\partial_x+\frac{\sigma^2}{2}\partial^2_x$. They define the operator $\hat{x}$ to be the multiplication-by-$x$ operator, that is a linear operator acting on some space of functions of the variable $x$, that simply multiplies any function by $x$. In the CIR model one restricts the analysis to the $x\geq0$ half-line due to the fact that the process cannot attain negative values. This means that the spectrum (set of eigenvalues $\lambda$) of $\hat{x}$ is continuous and consists of all $\lambda\geq0$. The corresponding eigenfunctions are normalized delta functions $$\hat{x}\delta(x-\lambda)=\lambda\delta(x-\lambda)$$ Going forward they extend this operator into $\hat{x}(\tau)=e^{H\tau}xe^{-H\tau}$ and define also $\hat{p}(\tau)=e^{H\tau}\partial_x e^{-H\tau}$ which are connected by a pair of ODEs. Solving them leads to a solution for $\hat{x}(\tau)=(a_2+b_2x)\partial^2_x+(a_1+b_1x)\partial_x+(a_0+b_0x)$. Then they denote the eigenfunction $Q(x,y,\tau)$ in general obeying $$(a_2+b_2x)\partial^2_xQ+(a_1+b_1x)\partial_xQ+(a_0+b_0x)Q=yQ \tag{2}$$ This equation determines the correct $x$-dependence of the density function but due to its linearity any solution could be multiplied by any function of $y$ and $\tau$ and would still be a solution. Hence this equation alone is not enough to get the correct dependence on the other two variables. In order to get the correct $\tau$-dependence one will have to substitute $Q$ back into the backward Kolmogorov equation and then make use of the initial condition $f(x,y,0)=\delta(x-y)$ to get the correct $y$-dependence.
Let's skip to page $5$ of the article section A. Vasicek Model: they get the following operator $$\hat{x}=\frac{\sigma^2}{2b}(e^{b\tau}-e^{-b\tau})\partial_x+\frac{a}{b}(1-e^{-b\tau})+e^{-b\tau}x$$ which should lead, according to (2), to the following ODE $$a_1\partial_xQ+(a_0+b_0x)Q=yQ \tag{3}$$ with $a_1=\frac{\sigma^2}{2b}(e^{b\tau}-e^{-b\tau})$, $a_0=\frac{a}{b}(1-e^{-b\tau})$ and $b_0=e^{-b\tau}$. They claim: " There is only one linearly independent solution for every $y \geq 0$ eigenvalue: $$Q(x,y,\tau)\sim exp\Biggl(-\frac{(y-a_0-b_0x)^2}{2e^{b\tau}a_1}\Biggr)" \tag{4}$$ The problem is that if I solve $(3)$, I'll get $Q=C\,exp(-\frac{b_{0}\,{x}^{2}}{2\,a_{1}}+\frac{y\,x}{a_{1}}-\frac{a_{0}\,x}{a_{1}})$ that doesn't have $y^2$ anywhere. Is $(4)$ wrong or what am I missing?
Also after that they use equation $(*)$ to get $$f(x,y,\tau)\sim \frac{1}{\sqrt{1-e^{-2b\tau}}}Q(x,y,\tau) \tag{5}$$ and claim " Putting everything together we see that the last piece missing from the density is a constant factor, $(2\pi\sigma^2/(2b))^{-1/2}$, in this case independent of $y$. This factor assures that the initial condition $f(x,y,0)=\delta(x-y)$ is met as well." This means that the limit for $\tau$ that tends to $0$ in $(5)$ should give $(2\pi\sigma^2/(2b))^{1/2}\delta(x-y)$, but I don't understand how the delta function should pop up in this limit.