2025年度

25/06/26 進捗報告

Grackleを用いて cooling timeを計算して、以下の(Correa et al. 2018)の図を作成しようとしました。

が。。。。計算時間がえぐく一旦中止

20 halosに対して3000万粒子ぐらいあるので、その 0.1% を計算すると18.3 sほどかかった。なので 20 haloの計算で、2日ほどかかることが分かりました。

そこでより簡易化を行うことにしました。

CloudyData_UVB=HM2012.h5 のテーブルを用いて線形補間を行って導出しようとしました。

CloudyData_UVB=HM2012.h5 の中身

h5pyを用いて構造を確認すると以下の通りでした。

Group: CoolingRates
Group: CoolingRates/Metals
Dataset: CoolingRates/Metals/Cooling
Dataset: CoolingRates/Metals/Heating
Group: CoolingRates/Primordial
Dataset: CoolingRates/Primordial/Cooling
Dataset: CoolingRates/Primordial/Heating
Dataset: CoolingRates/Primordial/MMW
Group: UVBRates
Group: UVBRates/Chemistry
Dataset: UVBRates/Chemistry/k24
Dataset: UVBRates/Chemistry/k25
Dataset: UVBRates/Chemistry/k26
Dataset: UVBRates/Chemistry/k27
Dataset: UVBRates/Chemistry/k28
Dataset: UVBRates/Chemistry/k29
Dataset: UVBRates/Chemistry/k30
Dataset: UVBRates/Chemistry/k31
Dataset: UVBRates/Info
Group: UVBRates/Photoheating
Dataset: UVBRates/Photoheating/piHI
Dataset: UVBRates/Photoheating/piHeI
Dataset: UVBRates/Photoheating/piHeII
Dataset: UVBRates/z

きっと以下のとおりでしょう：

• Cooling/Heating: おそらくガスやプラズマの冷却/加熱率データ。
• Primordial/Metals: 原始的元素（H, Heなど）・金属成分（C, O, Si等）ごとの情報でしょう。
• UVBRates: 紫外線背景（UVB）に関するデータセット。
• Chemistry/kXX: 化学反応に関するレート係数と思われます。
• Photoheating/piHI など: フォトン（紫外線など）による水素やヘリウムの加熱率。

このcooling tableを使用します。

cooling table

構造は以下の通りでした。

データセット名	shape	dtype	属性 (一部)
/CoolingRates/Metals/Cooling	(29, 26, 161)	float64	Dimension, Parameter1, Parameter2, Parameter3 など
/CoolingRates/Metals/Heating	(29, 26, 161)	float64	同上
/CoolingRates/Primordial/Cooling	(29, 26, 161)	float64	同上
/CoolingRates/Primordial/Heating	(29, 26, 161)	float64	同上
/CoolingRates/Primordial/MMW	(29, 26, 161)	float64	同上

また Metals と Primordial があります。

【Primordial（原始成分）】

• 「Primordial」はビッグバン直後に存在した主な元素、つまり水素（H）とヘリウム（He）だけを考慮しています。
• 宇宙の初期、まだ恒星が誕生する前の「原始ガス」の冷却・加熱率です。
• 分子や金属（C, O, Feなど）は含みません。
• 典型的には、温度や密度だけがパラメータになります。

【Metals（金属成分）】

• 「Metals」は水素・ヘリウム以外の重元素（天文学ではこれらを全部まとめて「金属」と呼びます）も含めた冷却・加熱率です。
• 金属は星の内部で生成され、超新星爆発などで宇宙空間にまき散らされます。
• 金属成分があると、原始成分に比べて低温域での冷却率が大きく上昇します（特にT < 10⁶ Kの領域）。
• このデータセットは、金属量（metallicity, Z）もパラメータになっていることが多いです。

テーブルの線形補間

今回、テーブルは CoolingRates/Metals/Cooling を使用します。

これは、水素数密度, 赤方偏移, 温度を変数として冷却率を返すテーブルになっています。そこで3次元グリッド（密度・赤方偏移・温度）上のデータなので、「三次元線形補間（trilinear interpolation）」を使うことにします。

# 水素数密度（hden）は0で固定
hden_val = 0.0
# 赤方偏移（z）は0で固定
z_val = 0.0
# 温度グリッド
temps = temp_grid

# 各温度に対する冷却率を計算
rates = [get_cooling_rate(hden_val, z_val, T) for T in temps]

plt.figure()
plt.loglog(temps, np.abs(rates))
plt.xlabel("Temperature [K]")
plt.ylabel("Cooling Rate [abs(erg cm^3/s)]")
plt.title("Cooling Rate vs Temperature (hden=0, z=0, log-log)")
plt.grid(True, which="both", ls="--")
plt.tight_layout()
plt.show()

少し範囲は異なりますが、値はまぁ悪くないでしょう。

次にcooling timeを計算するコードに移ります。

cooling time (ノート)

概要

Cloudyで生成された冷却テーブル（HDF5ファイル）を用い、 水素密度・温度・金属量・電子存在率を自由に指定して

• 冷却率 \([erg~cm^{-3}~s^{-1}]\)
• 冷却タイム \([{\rm Gyr}]\)

を高速・大量データ向けにNumPyベクトル化で算出するPythonコードの実装例です。

前提

• Cloudy等で生成したHDF5形式の冷却率データを/mnt/data/CloudyData_UVB=HM2012.h5として用いる
• 対応する冷却テーブルの軸は
  • 水素数密度（log10[cm\(^{-3}\)]）
  • 赤方偏移（z）
  • 温度（[K]）
• 金属量（Z）はテーブル内で太陽値を基準とし、外部スケールとして扱う

必要なパッケージ

import numpy as np
import h5py
from scipy.interpolate import RegularGridInterpolator

定数定義

m_p = 1.6726219e-24     # プロトン質量 [g]
X_H = 0.76              # 水素質量比（標準的宇宙組成）
k_B = 1.380649e-16      # ボルツマン定数 [erg/K]
sec_in_Gyr = 3.1536e16  # 1 Gyr [秒]

冷却テーブルの読込と補間器の構築

file_path = '/mnt/data/CloudyData_UVB=HM2012.h5'

with h5py.File(file_path, 'r') as f:
    data = np.array(f['CoolingRates/Metals/Cooling'])
    attrs = f['CoolingRates/Metals/Cooling'].attrs
    hden_grid = attrs['Parameter1']
    z_grid = attrs['Parameter2']
    temp_grid = attrs['Temperature']

# グリッドの昇順化と重複除去（scipyの補間仕様に合わせる）
z_order = np.argsort(z_grid)
z_grid_sorted = z_grid[z_order]
data_sorted = data[:, z_order, :]
z_grid_unique, unique_indices = np.unique(z_grid_sorted, return_index=True)
data_unique = data_sorted[:, unique_indices, :]

cooling_interp = RegularGridInterpolator(
    (hden_grid, z_grid_unique, temp_grid),
    data_unique,
    bounds_error=False,
    fill_value=None
)

主要関数

1. 水素数密度への変換

def calc_hden(rho):
    """
    質量密度 [g/cm^3] → 水素数密度 log10 [cm^-3] へ変換
    rho: 質量密度 [g/cm^3]（NumPy配列/float対応）
    戻り値: log10(n_H [cm^-3])
    """
    n_H = (X_H * rho) / m_p
    return np.log10(n_H)

2. 冷却率の計算

def get_cooling_rate_general(rho, T, Z, x_e, redshift=0.0):
    """
    Parameters
    ----------
    rho : ndarray
        ガス質量密度 [g/cm^3]
    T : ndarray
        温度 [K]
    Z : ndarray
        金属量（太陽値=1.0で正規化）
    x_e : ndarray
        電子存在率 (n_e / n_H)
    redshift : float, optional
        赤方偏移（z、既定値は0）

    Returns
    -------
    cooling_rate : ndarray
        体積あたり冷却率 [erg/cm^3/s]
    """
    hden = calc_hden(rho)
    pts = np.stack([hden, np.full_like(hden, redshift), T], axis=-1)
    Lambda = cooling_interp(pts)
    Lambda_scaled = Lambda * Z  # 金属量スケール
    n_H = 10 ** hden
    n_e = n_H * x_e
    return n_H * n_e * Lambda_scaled

3. 冷却タイムの計算

def cooling_time(rho, T, Z, x_e, redshift=0.0):
    """
    Parameters
    ----------
    rho : ndarray
        ガス質量密度 [g/cm^3]
    T : ndarray
        温度 [K]
    Z : ndarray
        金属量（太陽値=1.0で正規化）
    x_e : ndarray
        電子存在率 (n_e / n_H)
    redshift : float, optional
        赤方偏移（z、既定値は0）

    Returns
    -------
    cooling_time : ndarray
        冷却タイム [Gyr]
    """
    hden = calc_hden(rho)
    n_H = 10 ** hden
    n_e = n_H * x_e
    n_tot = n_H + n_e  # He/金属の厳密計算が必要なら適宜調整
    L = get_cooling_rate_general(rho, T, Z, x_e, redshift)
    E = (3 / 2) * n_tot * k_B * T  # 熱エネルギー密度 [erg/cm^3]
    tcool_sec = E / np.abs(L)
    tcool_gyr = tcool_sec / sec_in_Gyr
    return tcool_gyr

使い方例

# サンプル：10,000個の状態について一括計算
N = 10000
rho = np.full(N, 1.67e-24)       # [g/cm^3]（n_H=1/cm^3相当）
T = np.logspace(2, 8, N)         # 100K～10^8K
Z = np.ones(N)                   # 太陽金属量
x_e = np.full(N, 1.0)            # 完全電離（参考値）
redshift = 0.0

cool_rate = get_cooling_rate_general(rho, T, Z, x_e, redshift)
tcool = cooling_time(rho, T, Z, x_e, redshift)

print("冷却率shape:", cool_rate.shape)
print("冷却タイムshape:", tcool.shape)
print(f"T=10^4Kでの冷却タイム[Gyr]: {tcool[np.abs(T-1e4).argmin()]}")

Grackleとの比較

• 黒線がGrackleです
• 赤線が補完ver

全然違います。この図は一応、ノートにプラスして加熱についても考慮しましたが、駄目でした…

ちゃんとGrackleで計算しないとまずいのですかね…

やっぱりGrackleを使おう・・・✅