is non-zero only when the comoving separation between $(x,\theta,\phi)$ and $(x',\theta',\phi')$ is comparable to, or less than $R$. Practically, $R$ is much less than the Hubble scale, and the linear operator in the equation (\[eqa3-16-1\]) is commutable with the smoothing integration of equation (\[eq3-21\]), since the functions $D(z)$, $b(z)$, $f(z)$, $\alpha(z)$, and ${S_K}(x)$ are practically constants over individual cells, and the Laplacian is an Hermitian operator. Therefore, the equation (\[eq3-21\]) reduces to $$\delta^{\rm (s)}_{R}(x(z),\theta,\phi) =
  D(z)
  \left\{
    b(z) + f(z)
      \left[\frac{\partial}{\partial x} + \alpha(z)\right]
        \frac{\partial}{\partial x} (\triangle + 3K)^{-1}
  \right\}
  \delta_R(x(z),\theta,\phi),
\label{eq3-22}$$ where $$\delta_R(x,\theta,\phi) =
  \int {S_K}^2(x') dx' \sin\theta' d\theta' d\phi'
  W_{R}(x,\theta,\phi;x',\theta',\phi')
  \delta(x',\theta',\phi'),
\label{eq3-23}$$ is the smoothed density contrast in real space. A complex shape of the smoothing kernel complicates further analytic treatment, and therefore slows down the numerical evaluation of the correlation matrix based on the analytic formula. It is desirable to choose the smoothing kernel as simple as possible. In this repect, the simplest choice of the smoothing function $W_R$ is the spherically symmetric function in comoving space, in which case, the kernel is expanded by complete set of orthonormal modes $X_0$ (c.f., eq.\[\[eqa2-6b\]\]), $$W_R(x_{\rm rel}) =
  \int \frac{k^2dk}{2\pi^2} W(kR) X_0(k,x_{\rm rel}),
\label{eq3-24}$$ where $x_{\rm rel}$ is the 