從離散功率譜到星系相關函數：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

參考