从离散功率谱到星系相关函数：Hankel变换

在宇宙学中，星系数密度功率谱的计算方法一般是直接对三维坐标 $(x,y,z)$ 密度涨落进行FFT，这样可以得到三维 $(k_x, k_y, k_z)$ 傅里叶变换，但是为了更好的观察数据，一般会把三维FFT结果平均投影到以 $k = \sqrt{k_x^2 + k_y^2 + k_z^2}$ 为半径的球壳上，从而得到一维功率谱 $P(k)$ 即monopole 功率谱 $P_0(k)$，是角向平均的功率谱。

另外为了探测RSD(Redshift space distortion)效应，还会计算多级功率谱(但是这里我们只讨论阶数为0的情况)，这样计算效率会变得非常慢，并且各种观测效应的处理非常复杂，其中一种较为简便的方式利用到Hankel transform。

Hankel transform

假设密度涨落是 $f(r)$，相当于以 $r$ 为半径的球壳里的密度涨落，$r$ 可以看作共动距离，$f(r)$ 的hankel transform：
\begin{align}
F_\nu(k) = \int_0^{+\infty} f(r) J_\nu(k r)r \mathrm{d}r
\end{align}
这里的 $J_\nu$ 是阶数为 $\nu$ 的第一类贝塞尔函数，与球贝塞尔函数 $j_\nu$ 的关系是：
\begin{align}
j_\nu (x) = \sqrt{\frac{\pi}{2x}}J_{\nu+1/2}(x)
\end{align}

逆变换：
\begin{align}
f(r) = \int_0^{+\infty} F_\nu(k) J_\nu(k r)k \mathrm{d}k
\end{align}

当 $k$、$k^\prime$大于0时，满足正交性：
\begin{align}
\int_0^{+\infty} J_\nu(kr)J_\nu(k^\prime r)r \mathrm{d}r = \frac{\delta(k^\prime-k)}{k}
\end{align}

如果 $f(r)$ 满足球对称(实际上是 $f(r, \theta)$，但是满足球对称忽略了 $\theta$，本质上反映了角向平均的特性)，那么FFT可以直接利用Hankel transform表示：
\begin{align}
F_0(k) = 2\pi \int_0^{+\infty} f(r) J_0(k r)r \mathrm{d}r
\end{align}

星系功率谱与Hankel transform

既然上面提到Hankel transform可以计算FFT，你或许想到对三维星系数密度涨落进行Hankel transform然后乘以 $2\pi$ 计算FFT，但是这样做并不正确，两个原因：

Hankel transform的计算需要在极坐标里，而原始数据是三维 $(x,y,z)$ 的，因此不适用
星系功率谱是两个FFT的乘积，但是Hankel transform主要用于一维问题，而星系数密度的功率谱需要通过多维FFT来计算，因此不适用

那么Hankel transform有什么用呢？答案是可以实现功率谱与相关函数之间的转换

利用Hankel transform实现功率谱与相关函数之间的转换

假设功率谱是 $P(k)$，相关函数是 $\xi(r)$ (对应 $f(r)$)：

三维笛卡尔坐标系 $(x,y,z)$ 下，功率谱 $P(k_x,k_y,k_z)$ 与相关函数 $\xi(x,y,z)$ 之间的转换很简单，功率谱的逆FFT就是相关函数；反之，相关函数的FFT就是功率谱；

而球壳坐标下(极坐标系)，功率谱 $P(k)$ 与相关函数 $\xi(r)$ 之间的转换就是通过Hankel transform计算：
\begin{align}
P_\ell(k) = 4\pi (-i)^\ell \int_0^{+\infty} \xi_\ell(r) j_\ell(k r)r^2 \mathrm{d}r
\end{align}

对应逆变换：
\begin{align}
\xi_\ell(r) = 4\pi i^\ell \int_0^{+\infty} P_\ell(k) j_\ell(k r)k^2 \mathrm{d}k
\end{align}

这两个公式里 $\ell$ 代表阶数，等同于 $\nu$，$j_\ell$是球贝塞尔函数，注意区别，多级功率谱的意义是计算角向的各向异性。一般取0，即只计算角向平均功率谱即可。

利用Hankel transform计算相关函数的好处是，在计算多级功率谱时可以方便得到相关函数，从而考虑多级功率谱的window function效应。

python实现Hankel transform

程序上实现的方法一个是使用 hankel 库计算，另一方面是写函数手算，这里建议手写函数计算，因为计算星系功率谱与相关函数的转换必须要手动计算。

导入库并定义函数：

import numpy as np
import matplotlib.pyplot as plt

import scipy.special as sp
from scipy.integrate import quad,simps
from scipy.interpolate import interp1d

import hankel
from hankel import HankelTransform

# 定义手动计算hankel transform和逆变换函数
def hankel_transform(fr_dis, r, k, nu):
    
    hankeled = np.zeros_like(k)
                
    for ii,i in enumerate(k):
        coff_bessel = fr_dis*sp.jv(nu, i*r)*r
        
        # hankeled[ii] = quad(func, 0, np.inf)
        hankeled[ii] = simps(coff_bessel, r)
        
    return hankeled

def ihankel_transform(fk_dis, r, k, nu):

    ihankeled = np.zeros_like(r)
                
    for ii,i in enumerate(r):
        coff_bessel = fk_dis*sp.jv(nu, k*i)*k
        
        # hankeled[ii] = quad(func, 0, np.inf)
        ihankeled[ii] = simps(coff_bessel, k)
        
    return ihankeled

进行测试：

# 定义k比较符合物理
k = np.linspace(0.02, 1.0, 20)
r = 2*np.pi/k   # 单调递减
r = r[::-1]   # 为了插值准确需要变成单点递增

# 离散点数太少，为了使得计算更精确，进行插值扩充点
nr = 200  # 扩充到100个点
r_dense = np.linspace(r.min(), r.max(), nr)

k_dense = 2*np.pi/r_dense   # 最好是这样计算，保证k-r一一对应即可，比如log分bin必须一致
k_dense = k_dense[::-1]


f_r_func = lambda r: 1/(2*r)  # 连续函数形式
f_r = f_r_func(r)


# 一定注意，插值保证自变量x单调递增，否则插值可能不准确
f_r_inte = interp1d(r, f_r, kind='linear', bounds_error=False, fill_value=0)
f_r_dense = f_r_inte(r_dense) # 离散插值数据


# 手动计算HF，离散数据
hfr_dis = hankel_transform(f_r_dense, r_dense, k_dense, 0)
ihfr_dis = ihankel_transform(hfr_dis, r_dense, k_dense, 0)

# hankel程序
ht = HankelTransform()
hfr = ht.transform(f_r_inte, k_dense)

对比正向变换：

plt.plot(k_dense, hfr_dis, label='Discrete manual')
plt.plot(k_dense, hfr[0], label='hankel software')
plt.legend()
plt.xlabel('k')
plt.ylabel('f(k)')

对比逆向变换：

plt.plot(r, f_r, label='origin')
plt.plot(r_dense, f_r_dense, label='interp1d')
plt.plot(r_dense, ihfr_dis, label='Reconstructed ')
plt.legend()
plt.xlabel('r')
plt.ylabel('f(r)')

从离散功率谱到星系相关函数：Hankel变换

Hankel transform

星系功率谱与Hankel transform

利用Hankel transform实现功率谱与相关函数之间的转换

python实现Hankel transform

参考