1. 要旨
中性子星は超新星爆発時の非対称性によって大きな固有速度を獲得すると考えられており、その速度分布は爆発機構を探る手がかりとなる。本記事では、ATNF Pulsar Catalogue を psrqpy で取得し、固有運動と距離から孤立パルサーの横方向速度 $v_\perp$ を計算した。
若い視差距離パルサーに限ると、得られた速度分布は Hobbs et al. (2005) の標準的な Rayleigh 分布よりも低速側に寄って見える。log-normal 分布との比較に加え、カット条件のロバスト性と、線形ビン・対数ビンでピーク位置の解釈が変わる Jacobian の効果を確認する。結論を急がず、公開カタログからパルサーキック速度分布の議論をどこまで再現できるかを検証してみる。
2. なぜパルサー速度分布を見るのか
中性子星は、重力崩壊型超新星爆発の後に残される高密度天体である。爆発が完全に球対称であれば、残された中性子星に大きな反跳速度は生じにくい。しかし実際には、多くのパルサーが数百 km/s 程度の固有速度を持つ。この速度は、中性子星が誕生時に非対称な力を受けた痕跡と考えられ、一般に natal kick と呼ばれる。
パルサーの速度分布は、この natal kick を統計的に調べるための基本的な観測量である。個々の天体の速度だけでなく、多数のパルサーを集めたときにどのような速度スケールを持つのか、低速成分や高速成分がどの程度存在するのかを見ることで、中性子星形成過程に対する制約を与えることができる。
Figure1 : Hobbs et al. (2005) によるパルサー速度分布の解析。1次元・2次元速度の観測分布から、3次元速度分布が Maxwellian でよく記述されることを示している。
この分野で特によく参照されるのが Hobbs et al. (2005) の結果である。同研究では、パルサーの 3 次元速度分布が 1 成分の マクスウェル分布で近似できるとされ、1 次元速度分散として $\sigma = 265 \ \mathrm{km \ s^{-1}}$ が与えられた。この値はその後、中性子星キック速度の標準的なスケールとして多くの研究で用いられてきた。
一方で、観測から直接得られるのは、通常は視線方向を含む 3 次元速度ではなく、固有運動と距離から計算される横方向速度 $v_\perp$ である。また、距離推定の方法やサンプル選択によって、得られる速度分布は変化しうる。ただ、近年の視差距離を用いた解析では、若い孤立パルサーの速度分布が Hobbs et al. (2005) の分布よりも低速側に寄る可能性も議論されている（Disberg & Mandel 2025）。本記事はこの問題意識を出発点として、公開カタログから同様の傾向が再現されるかを確認する。
この記事ではATNF Pulsar Catalogue から固有運動と距離を取得し、孤立パルサーの横方向速度 $v_\perp$ を計算する。その上で、Hobbs et al. (2005) に基づく Rayleigh 分布と log-normal 分布を比較し、得られる速度分布がサンプル選択や表示方法にどの程度依存するのかを確認する。
3. 横方向速度とRayleigh分布
この記事で扱う速度は、中性子星の真の3次元空間速度ではなく天球面上の運動から得られる横方向速度である。パルサーの固有運動を $\mu$、距離を $d$ とすると、横方向速度 $v_\perp$ は
$$
v_\perp =
4.74047
\left(\frac{\mu}{\mathrm{mas \ yr^{-1}}}\right)
\left(\frac{d}{\mathrm{kpc}}\right)
\ \mathrm{km \ s^{-1}}
$$
で与えられる。ここで係数 $4.74047$ は、1 AU/yr を km/s に換算した値である。固有運動の大きさは、赤経・赤緯方向の固有運動成分から
$$
\mu =
\sqrt{
\mu_\alpha^2 + \mu_\delta^2
}
$$
として計算した。したがって、本解析で得られる $v_\perp$ は、視線速度成分を含む真の速度
$$
v =
\sqrt{
v_\perp^2 + v_r^2
}
$$
そのものではない。多くのパルサーでは視線速度 $v_r$ が測定されていないため、本記事では観測的に直接構成しやすい $v_\perp$ の分布を比較対象とする。
Hobbs et al. (2005) では、パルサーの3次元速度分布が1成分の マクスウェル分布で近似できるとされ、1次元速度分散として
$$
\sigma = 265 \ \mathrm{km \ s^{-1}}
$$
が与えられた。3次元速度の大きさ $v$ に対する マクスウェル分布は
$$
f_{\mathrm{Maxwell}}(v)=
\sqrt{\frac{2}{\pi}}
\frac{v^2}{\sigma^3}
\exp\left(
-\frac{v^2}{2\sigma^2}
\right)
$$
である。一方、本記事で直接比較するのは横方向速度 $v_\perp$ である。速度の各成分が同じ分散 $\sigma^2$ を持つ等方的な Gaussian に従うなら、2次元射影速度 $v_\perp$ は Rayleigh 分布
$$
f_{\mathrm{Rayleigh}}(v_\perp)=
\frac{v_\perp}{\sigma^2}
\exp\left(
-\frac{v_\perp^2}{2\sigma^2}
\right)
$$
に従う。したがって本記事では、Hobbs et al. (2005) の速度分布を横方向速度に射影したものとして、$\sigma=265 \ \mathrm{km \ s^{-1}}$ の Rayleigh 分布を基準モデルに用いる。
比較対象として、log-normal 分布も考える。これは、速度そのものではなく速度の対数が正規分布に従うと仮定するモデルであり、
$$
\ln v_\perp \sim \mathcal{N}(m \ s^2)
$$
と書ける。このとき確率密度関数は
$$
f_{\mathrm{LN}}(v_\perp)=
\frac{1}{v_\perp s\sqrt{2\pi}}
\exp\left[
-\frac{(\ln v_\perp - m)^2}{2s^2}
\right]
$$
である。log-normal 分布の最頻値は
$$
v_{\mathrm{mode}}=
\exp(m-s^2)
$$
で与えられる。本記事では、若い視差距離サンプルに対して $(m,s)$ を最尤推定し、Hobbs 型 Rayleigh 分布との比較を行う。
分布の一致度は、主に経験的累積分布関数とモデル累積分布関数の差によって評価する。観測値から得られる経験的 CDF を $F_{\mathrm{emp}}(v)$、モデル CDF を $F_{\mathrm{model}}(v)$ とすると、KS 統計量は
$$
D=\sup_v
\left|
F_{\mathrm{emp}}(v)-
F_{\mathrm{model}}(v)
\right|
$$
で定義される。さらに、最尤値に基づく
$$
\mathrm{AIC}=
2k - 2\ln \mathcal{L}_{\max}
$$
も補助的な指標として用いる。ここで $k$ はモデル中でデータから推定したパラメータ数である。この記事では、Hobbs 型 Rayleigh 分布については $\sigma=265 \ \mathrm{km \ s^{-1}}$ を固定し、本解析のデータにはフィットしないため $k=0$ とする。一方、log-normal 分布では $(m,s)$ を最尤推定するため $k=2$ とする。
この扱いは、Hobbs 型 Rayleigh 分布に対してパラメータ罰則を課さないという意味で、むしろ Rayleigh モデルに有利な保守的比較である。したがって、AIC において log-normal 分布が同等または優位になる場合、それは2つの自由パラメータの罰則を受けてもなお、観測分布をよく記述していることを意味する。ただし、サンプル数は限られているため、KS 検定や AIC のみから速度分布の真の形を断定することは避ける。
なお、速度分布のピークを議論する際には、表示方法にも注意が必要である。線形速度 $v$ に対する密度は $f(v)$ であるが、対数速度 $\ln v$ のビンで分布を見る場合、対応する密度は
$$
\frac{dP}{d\ln v}=
v f(v)
$$
となる。したがって同じ物理モデルであっても、線形ビンで見るか対数ビンで見るかによって、見かけのピーク位置は変化する。この Jacobian の効果については、後の節で改めて可視化する。
4. ATNF + psrqpy とサンプル選別
本解析では、パルサーの固有運動・距離・年齢などの情報を ATNF Pulsar Catalogue から取得した。データ取得には Python パッケージ psrqpy を用いた。今回使用したキャッシュは 2026-06-11 14:14 に取得したもので、カタログバージョンは v2.8.0 である。
ATNF カタログは、電波パルサーを中心に多様な中性子星種族を含む非常に有用なデータベースである。一方で、本記事の目的は「若い孤立パルサーの natal kick に由来する速度分布」を調べることであり、カタログ全体をそのまま一つの母集団として扱うことは適切ではない。特に、連星パルサーやミリ秒パルサーは、通常の若い孤立パルサーとは進化史が異なる。
Figure2 : Manchester et al. (2005) Figure 4。ATNF Pulsar Catalogue に含まれるパルサーの周期分布。通常の孤立パルサー、連星パルサー、AXP などが異なる周期領域に分布することが分かる。
この図が示すように、パルサーとして同じカタログに登録されていてもその集団は一様ではない。ミリ秒パルサーは長い進化の中でリサイクルを受けた天体であり、連星パルサーは伴星との相互作用や連星進化の影響を受けている。また、AXP や XINS などは通常の電波パルサーとは観測的・物理的性質が異なる。したがって、若い中性子星の誕生時キックを調べるためには、これらを同じ速度分布の母集団に含めないようにサンプルを選別する必要がある。
そこで本記事では、まず固有運動の赤経・赤緯成分がともに測定されている天体を抽出し、その後、連星、ミリ秒パルサー、球状星団内のパルサー、AXP・XINS などを除外した。さらに、横方向速度の計算に必要な距離情報を持つ天体のみを残し、固有運動がゼロとして登録されている placeholder 的な天体も除外した。
距離の扱いについても注意が必要である。本記事では、全サンプルの概観やカット条件の確認には、ATNF カタログの DIST に基づく横方向速度を用いた。ATNF の標準パラメータでは、これは VTRANS に対応する量である。ここで DIST は単純な DM 距離ではなく、DIST_A などの独立距離推定、有意な視差がある場合の $1/\pi$、距離上下限を用いた推定などを優先して設定される derived quantity である。
一方、本記事の主解析である young PX sample では、視差 PX を持つ天体に限定し、距離を単純に
$$
d = \frac{1}{\pi}
$$
として計算した。ここで $\pi$ は mas 単位の視差であり、$d$ は kpc 単位である。この距離は最も単純な点推定であり、非対称な距離事後分布や Lutz-Kelker 型のバイアス補正を取り込んだものではない。したがって、視差サンプルで見える低速側への寄りは、物理的な速度分布だけでなく、視差距離推定に伴うバイアスの影響を受けている可能性がある。この点は考察で改めて述べる。
このようにして、最終的に 116 天体からなる孤立パルサーサンプルを構成した。以後、このサンプル全体に加えて、特性年齢 $\tau < 10 \ \mathrm{Myr}$ の若いサンプル、およびその中で視差距離を持つサンプルを主に用いて、横方向速度分布を比較する。
5. P–Ṗ図によるサンプルの確認
速度分布の比較に進む前に、本解析で選んだ天体がパルサー集団の中でどのような位置を占めているかを確認する。ここでは、ATNF カタログ中のパルサーを背景に置き、本解析で用いたサンプルをP–Ṗ図上に重ねて表示した。
特性年齢は、磁気双極子放射による単純なスピンダウンを仮定すると
$$
\tau =
\frac{P}{2\dot{P}}
$$
で与えられる。なお、観測される Ṗ には、パルサーの横方向運動によって生じる見かけの増加、いわゆる Shklovskii 効果が含まれうる（Shklovskii 1970）。この効果は特性年齢の評価に影響する可能性があるため、本解析では補正量を別途見積もった。ただし、本サンプルでは影響は多くの天体で小さく、基準的な young / old 分類を大きく変えるものではなかったため、主解析では観測値に基づく特性年齢を用いた。本記事では、$\tau < 10 \ \mathrm{Myr}$ を young sample、さらに $\tau < 100 \ \mathrm{kyr}$ を very young sample として区別した。なお、very young sample は young sample の部分集合である。
Figure 3: P–Ṗ図。灰色の点は ATNF カタログ中のパルサー、色付きの点は本解析で用いた孤立パルサーサンプルを表す。淡青は old sample（$\tau > 10 \ \mathrm{Myr}$, $N=36$）、青は young sample のうち $\tau < 10 \ \mathrm{Myr}$ かつ $\tau \geq 100 \ \mathrm{kyr}$ の天体（$N=69$）、黒は very young sample（$\tau < 100 \ \mathrm{kyr}$, $N=11$）を示す。破線は等特性年齢線である。
図から分かるように、本解析で用いた young sample は、主に通常パルサーが分布する領域に位置しており、ミリ秒パルサーが集中する左下の領域とは明確に分かれている。これは、前節で行った MSP 除外や連星除外によって、リサイクルを受けた古いパルサーではなく、若い孤立パルサーを主に抽出できていることを示している。
また、very young sample は P–Ṗ図上で高いṖを持つ領域に集中しており、等特性年齢線の若い側に分布している。これは、特性年齢による分類が図上でも自然な分離として現れていることを確認するものである。
したがって、本解析で構成したサンプルは、単に固有運動と距離が測定されている天体を機械的に集めたものではなく、P–Ṗ図上でも若い孤立パルサーとして解釈しやすい集団になっている。次節では、このサンプルに対して横方向速度分布を比較する。
6. Hobbs Rayleigh と log-normal MLE
前節で構成したサンプルに対して、横方向速度 $v_\perp$ の分布を調べる。ここでは、まず全サンプルと若いパルサーの速度分布を確認し、その後、視差距離を持つ若いサンプルに対して Hobbs et al. (2005) に基づく Rayleigh 分布と log-normal 分布を比較する。
本記事では、特性年齢
$$
\tau < 10 \ \mathrm{Myr}
$$
を満たす天体を young sample と呼ぶ。また、P–Ṗ図上の表示では $\tau < 100 \ \mathrm{kyr}$ の天体を very young として分けているが、速度分布の比較ではこれらも young sample に含める。最終的に、全サンプルは $N=116$、young sample は $N=80$、そのうち視差距離を持つ young PX sample は $N=45$ である。
Figure 4: 孤立パルサーの横方向速度分布
(a) 全サンプル $N=116$ と young sample（$\tau < 10 \ \mathrm{Myr}$, $N=80$）の比較。淡青が全サンプル、青が young sample、赤破線が Hobbs et al. (2005) の $\sigma=265 \ \mathrm{km \ s^{-1}}$ Rayleigh 分布を表す。
(b) 視差距離を持つ young sample（young PX sample, $N=45$）に対するモデル比較。灰色のヒストグラムが観測された $v_\perp$ 分布、赤破線が Hobbs et al. (2005) Rayleigh 分布、青実線が young PX sample に対して最尤推定した log-normal 分布を表す。
全サンプルおよび young sample の横方向速度は、数百 $\mathrm{km \ s^{-1}}$ 程度に広く分布している。一方で、視差距離を持つ young sample に限ると、分布の中心は Hobbs et al. (2005) の $\sigma=265 \ \mathrm{km \ s^{-1}}$ Rayleigh 分布よりも低速側に位置しているように見える。これは、距離推定の不確かさが比較的小さい視差サンプルでは、従来よく用いられてきた速度スケールよりも低い速度成分が強調される可能性を示唆する。
この印象を定量化するため、young PX sample に対して log-normal 分布を最尤推定した。Hobbs 型 Rayleigh 分布では $\sigma=265 \ \mathrm{km \ s^{-1}}$ を固定し、本解析のデータにはフィットしない。一方、log-normal 分布では $(m,s)$ をデータから推定する。したがって、AIC の計算では Rayleigh 分布に対して $k=0$、log-normal 分布に対して $k=2$ とした。
経験的累積分布関数を用いた比較を以下に示す。
Figure 5: young PX sample の経験的 CDF と、Hobbs et al. (2005) Rayleigh 分布および log-normal MLE の比較。
young PX sample の経験的 CDF は、Hobbs 型 Rayleigh 分布よりも低速側で早く立ち上がる。すなわち、同じ累積割合に到達する速度が Hobbs 型分布より小さい。log-normal MLE はこの立ち上がりを比較的よく追っており、特に分布の中心付近では経験的 CDF とのずれが小さい。
今回のモデル比較を表にまとめる。
| モデル | KS 統計量 $D$ | KS $p$ 値 | AIC |
|---|---|---|---|
| Hobbs+05 Rayleigh（$\sigma=265 \ \mathrm{km \ s^{-1}}$ 固定） | 0.178 | 0.103 | 597.3 |
| log-normal MLE | 0.124 | 0.460 | 596.3 |
KS 統計量は log-normal MLE の方が小さく、経験的 CDF との最大ずれは Hobbs 型 Rayleigh 分布より小さい。また、AIC でも log-normal MLE の方がわずかに小さい。ただし、その差は大きくない。さらに、KS 検定において Hobbs 型 Rayleigh 分布も棄却されていないため、この結果だけから速度分布の真の形を断定することはできない。
重要なのは、young PX sample では Hobbs et al. (2005) の標準的な速度スケールよりも低速側に分布が寄って見える、という点である。
7. カット条件ロバスト性
ここまでの解析では、基準サンプルとして「$\tau < 10 \ \mathrm{Myr}$、$P_0 > 30 \ \mathrm{ms}$、視差 S/N $>3$、視差距離のみ」を満たす young PX sample を用いた。しかし、速度分布の比較はサンプル選択に依存しうる。特に、若いパルサーを定義する特性年齢の閾値、MSP を除外する周期閾値、視差距離の信頼度、DM 距離を含めるかどうかは、得られる速度分布に直接影響する。
そこで、基準条件から一つずつカット条件を変更し、サンプル数 $N$、横方向速度の中央値、log-normal MLE の最頻値、KS 検定の $p$ 値、AIC の差を比較した。ここでは $\Delta \mathrm{AIC} = \mathrm{AIC}{\mathrm{Hobbs}} - \mathrm{AIC}{\mathrm{lognormal}}$ と定義する。と定義する。したがって、$\Delta \mathrm{AIC} > 0$ であれば log-normal モデルの方が AIC は小さく、観測分布をよりよく記述していることを意味する。
Table1 : サンプル選択条件を変えたときの速度分布比較。Baseline は $\tau < 10 \ \mathrm{Myr}$、$P_0 > 30 \ \mathrm{ms}$、視差 S/N $>3$、視差距離のみを用いた young PX sample である。Mode は log-normal MLE の最頻値を表す。
表から分かるように、log-normal MLE の最頻値は多くの条件で $100 \ \mathrm{km \ s^{-1}}$ 前後に現れる。基準条件では mode は $124 \ \mathrm{km \ s^{-1}}$ であり、年齢閾値を $\tau < 3 \ \mathrm{Myr}$ や $\tau < 40 \ \mathrm{Myr}$ に変えても、それぞれ $96 \ \mathrm{km \ s^{-1}}$、$120 \ \mathrm{km \ s^{-1}}$ と大きくは変化しない。
MSP 除外閾値を $P_0 > 10 \ \mathrm{ms}$ や $P_0 > 20 \ \mathrm{ms}$ に変更しても結果は基準条件と同一であり、この条件は今回の young PX sample にはほとんど影響しない。また、視差 S/N 閾値を $2$ または $5$ に変えても、中央値や mode は大きく変化しない。
一方でDM 距離を含めると $\Delta \mathrm{AIC}$ は大きく正になり、log-normal 側がより強く好まれる。ただし、DM 距離は電子密度モデルに依存するため、距離推定の扱いが速度分布の形に大きく効くことを示すものとして読むべきである。
8. Jacobianとピーク位置の解釈
ここまでの議論では、速度分布のピークが Hobbs et al. (2005) の速度スケールより低速側に見えることを確認してきた。しかし、速度分布のピーク位置を議論する際には、もう一つ注意すべき点がある。それは、同じ分布であっても、線形速度 $v$ に対する密度として見るか、対数速度 $\log v$ に対する密度として見るかによって、ピーク位置が変わるという点である。
線形速度 $v$ に対する確率密度を $f(v)$ とする。このとき、ある範囲 $[v, v+dv]$ に天体が入る確率は
$$
dP = f(v) \ dv
$$
で与えられる。一方で、速度を対数軸で見る場合、変数は $v$ ではなく $x=\log_{10}v$ である。このとき、確率密度は $f(v)$ そのものではなく、変数変換に伴う Jacobian を含んだ
$$
p(x)=
\frac{dP}{dx}=
f(v)\frac{dv}{dx}
$$
で与えられる。$x=\log_{10}v$ なので、
$$
\frac{dv}{dx}=
\ln(10) \ v
$$
であり、したがって
$$
p(\log_{10}v)=
\ln(10) \ v f(v)
$$
となる。定数 $\ln(10)$ は規格化にだけ効くため、ピーク位置を考える上では本質的には $v f(v)$ が重要である。
これは、ヒストグラムを対数ビンで描くときにも同じである。線形ビンのヒストグラムでは、縦軸はおおまかに単位速度あたりの確率密度 $f(v)$ を表す。一方、対数ビンのヒストグラムでは、縦軸は1 dex あたりの確率密度 $p(\log_{10}v)$ を表す。したがって、対数ビンのヒストグラムに線形密度 $f(v)$ をそのまま重ねると、比較している量が一致しない。
具体的に Hobbs 型 Rayleigh 分布を考える。線形速度に対する Rayleigh 分布は
$$
f(v)=
\frac{v}{\sigma^2}
\exp\left(
-\frac{v^2}{2\sigma^2}
\right)
$$
であり、この $f(v)$ のピークは
$$
v = \sigma
$$
にある。したがって、Hobbs et al. (2005) の $\sigma=265 \ \mathrm{km \ s^{-1}}$ を用いると、線形密度としてのピークは $265 \ \mathrm{km \ s^{-1}}$ になる。
しかし、対数速度に対する密度では $v f(v)$ を見ることになる。この場合、ピークは
$$
v = \sqrt{2}\sigma
$$
に移動する。したがって同じ Hobbs 型 Rayleigh 分布であっても、対数ビンで見たピークは
$$
\sqrt{2}\times 265
\simeq
375 \ \mathrm{km \ s^{-1}}
$$
となる。
log-normal 分布についても同様である。線形速度に対する密度 $f(v)$ の最頻値は $\exp(m-s^2)$ であるが、対数速度に対する密度では $v f(v)$ を見るため、ピークは $\exp(m)$ に移動する。つまり、log-normal 分布の「ピーク速度」は、線形密度を見ているのか、対数密度を見ているのかによって異なる。
この違いは単なる作図上の細部ではない。速度分布のピークを比較するとき、Hobbs 分布のピークはどこかlog-normal 分布はどの速度に集中しているのか、という議論そのものに影響する。そこで次の図では、同じ young PX sample と同じ2つのモデルを、線形ビンと対数ビンの両方で表示し、Jacobian によってピーク位置がどのように変わるかを確認する。
Figure 6: Jacobian によるピーク位置の違い。
(a) 線形速度 $v_\perp$ に対する密度 $f(v)$。灰色のヒストグラムは young PX sample、赤破線は Hobbs et al. (2005) Rayleigh 分布、青実線は log-normal MLE を示す。線形密度ではピーク位置はそれぞれ $265 \ \mathrm{km \ s^{-1}}$ と $124 \ \mathrm{km \ s^{-1}}$ である。
(b) 同じ分布を対数ビンで表示したもの。対数速度に対する密度では $dP/d\log v \propto v f(v)$ となるため、ピーク位置は線形表示から移動し、Hobbs Rayleigh 分布では $375 \ \mathrm{km \ s^{-1}}$、log-normal MLE では $226 \ \mathrm{km \ s^{-1}}$ となる。
この図から分かるように、同じモデルであっても線形ビンで見るか対数ビンで見るかによってピーク位置が変化する。特に Hobbs 型 Rayleigh 分布では、線形密度 $f(v)$ のピークは $265 \ \mathrm{km \ s^{-1}}$ であるが、対数速度に対する密度では $375 \ \mathrm{km \ s^{-1}}$ に移動する。
一方、log-normal MLE でもピークは $124 \ \mathrm{km \ s^{-1}}$ から $226 \ \mathrm{km \ s^{-1}}$ に移動する。ただし、どちらの表示でも log-normal MLE のピークは Hobbs 型 Rayleigh 分布より低速側にある。したがって、本解析で見えている低速側ピークは、単にビンの取り方だけで生じたものではない。
重要なのは、速度分布のピークを議論するときには、何に対する密度を見ているのかを明示する必要があるという点である。線形密度 $f(v)$ と対数密度 $dP/d\log v$ を混同すると、同じ分布に対して異なるピーク速度を比較してしまう。
9. 考察・まとめ
本記事では、ATNF Pulsar Catalogue v2.8.0 を psrqpy で取得し、固有運動と距離から孤立パルサーの横方向速度 $v_\perp$ を計算した。その上で、若い視差距離パルサーを中心に、Hobbs et al. (2005) の $\sigma=265 \ \mathrm{km \ s^{-1}}$ Rayleigh 分布と log-normal MLE を比較した。
基準サンプルである young PX sample では、経験的 CDF が Hobbs 型 Rayleigh 分布よりも低速側で早く立ち上がり、log-normal MLE の方が分布の中心付近をややよく再現した。また、カット条件を変えた場合でも、log-normal MLE のピークは多くの場合 $100 \ \mathrm{km \ s^{-1}}$ 前後に現れた。したがって、若い孤立パルサーの横方向速度分布が、Hobbs et al. (2005) の標準的な速度スケールより低速側に寄って見えるという傾向は、基準条件だけに依存したものではなさそうである。
ただし、この結果は「パルサーキック速度分布は log-normal である」と断定するものではない。young PX sample は $N=45$ と小さく、KS 検定では Hobbs 型 Rayleigh 分布も棄却されない。また、AIC の差も基準条件では大きくない。さらに、本解析で扱ったのは真の3次元速度ではなく横方向速度 $v_\perp$ であり、視線速度成分は含まれていない。加えて、視差サンプルそのものが距離推定に関する選択効果を持つ点にも注意が必要である。視差測定には Lutz--Kelker 型のバイアスがあり、単純な視差距離は距離を過小評価し、その結果として $v_\perp = 4.74\mu d$ も過小評価する方向に働きうる（Lutz & Kelker 1973; Verbiest et al. 2010, 2012）。したがって、本記事で見えた低速側への寄りの一部は、物理的なキック速度分布だけでなく、視差サンプル固有のバイアスによって説明される可能性がある。
本記事で特に重要だった点は、速度分布の比較が、サンプル選択と表示方法の両方に依存することである。距離サンプルに DM 距離を含めるかどうかで AIC 差は大きく変化し、また線形密度 $f(v)$ と対数密度 $dP/d\log v$ ではピーク位置の解釈も変わる。したがって、速度分布の「ピーク」や「モデルの良さ」を議論する際には、どのサンプルに対して、どの変数の密度を比較しているのかを明示する必要がある。
以上より、本解析は Hobbs et al. (2005) の速度分布を否定するものではなく、公開カタログから再現可能な範囲で、若い孤立パルサーの横方向速度分布が低速側に重心を持って見えることを確認した、という位置づけである。公開データを用いた単純な再解析であっても、サンプル選別、統計モデル、表示方法を丁寧に扱うことで、論文の主張を自分の手で検証するところまで進めることができた。
10. 参考文献
-
Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. (2005).
A statistical study of 233 pulsar proper motions.
Monthly Notices of the Royal Astronomical Society, 360, 974--992.
arXiv:astro-ph/0504584, doi:10.1111/j.1365-2966.2005.09087.x -
Disberg, P., & Mandel, I. (2025).
The Kick Velocity Distribution of Isolated Neutron Stars.
arXiv preprint.
arXiv:2505.22102 -
Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. (2005).
The Australia Telescope National Facility Pulsar Catalogue.
The Astronomical Journal, 129, 1993--2006.
arXiv:astro-ph/0412641, doi:10.1086/428488 -
ATNF Pulsar Catalogue.
Australia Telescope National Facility, CSIRO.
本記事ではpsrqpyのキャッシュにより、2026-06-11 14:14 に取得された ATNF Pulsar Catalogue v2.8.0 を用いた。
ATNF Pulsar Catalogue -
Pitkin, M. (2018).
psrqpy: a python interface for querying the ATNF pulsar catalogue.
The Journal of Open Source Software, 3(22), 538.
arXiv:1806.07809, doi:10.21105/joss.00538 -
psrqpydocumentation.
The psrqpy package: A Python tool for interacting with the ATNF pulsar catalogue.
https://psrqpy.readthedocs.io/ -
Cordes, J. M., & Lazio, T. J. W. (2002).
NE2001. I. A New Model for the Galactic Distribution of Free Electrons and its Fluctuations.
arXiv preprint.
arXiv:astro-ph/0207156 -
Cordes, J. M., & Lazio, T. J. W. (2003).
NE2001. II. Using Radio Propagation Data to Construct a Model for the Galactic Distribution of Free Electrons.
arXiv preprint.
arXiv:astro-ph/0301598 -
Yao, J. M., Manchester, R. N., & Wang, N. (2017).
A New Electron-density Model for Estimation of Pulsar and FRB Distances.
The Astrophysical Journal, 835, 29.
arXiv:1610.09448, doi:10.3847/1538-4357/835/1/29 -
Shklovskii, I. S. (1970).
Possible Causes of the Secular Increase in Pulsar Periods.
Soviet Astronomy, 13, 562.
ADS -
Lutz, T. E., & Kelker, D. H. (1973).
On the Use of Trigonometric Parallaxes for the Calibration of Luminosity Systems: Theory.
Publications of the Astronomical Society of the Pacific, 85, 573.
doi:10.1086/129506 -
Verbiest, J. P. W., Lorimer, D. R., & McLaughlin, M. A. (2010).
Lutz--Kelker bias in pulsar parallax measurements.
Monthly Notices of the Royal Astronomical Society, 405, 564--572.
arXiv:1002.1213, doi:10.1111/j.1365-2966.2010.16488.x -
Verbiest, J. P. W., Weisberg, J. M., Chael, A. A., Lee, K. J., & Lorimer, D. R. (2012).
On Pulsar Distance Measurements and Their Uncertainties.
The Astrophysical Journal, 755, 39.
arXiv:1206.0428, doi:10.1088/0004-637X/755/1/39 -
Akaike, H. (1974).
A new look at the statistical model identification.
IEEE Transactions on Automatic Control, 19, 716--723.
doi:10.1109/TAC.1974.1100705
11. Append
本文の結論は変わらないが、注意が必要な点を以下にまとめる。
A.1 log-normal 側の KS p 値はそのままでは正しくない（Lilliefors 問題）
本文の KS 検定には非対称性がある。Hobbs 型 Rayleigh 分布は $\sigma$ を外部の文献値に固定しており、データから何も推定していないため、標準的な KS 検定の $p$ 値がそのまま使える。一方、log-normal 分布は $(m, s)$ を検定対象と同じデータから最尤推定している。この場合、KS 統計量 $D$ の帰無分布は標準のものより小さい側に偏るため、標準の $p$ 値（本文の 0.460）はモデルの適合を過大評価する。正規分布に対するこの問題は Lilliefors 検定として知られる。
A.2 測定誤差を伝播させていない
本解析の $v_\perp$ は、固有運動と距離の点推定値から計算したものであり、測定誤差を伝播させていない。視差 S/N $>3$ というカットは、視差の相対誤差が最大 30% 程度のサンプルを許容しており、これは距離・速度の相対誤差としてもほぼそのまま効く。
これには2つの含意がある。第一に、観測されたヒストグラムは真の速度分布に測定誤差が畳み込まれたものであり、真の分布より広がって見える。第二に、正の量に乗法的な誤差が乗ると分布は右に歪むため、log-normal 的な形が誤差構造によって部分的に「作られて」いる可能性も排除できない。厳密には、誤差を含めた階層ベイズモデルで真の速度分布を逆畳み込みすべきであり、これは本記事の範囲を超える。本記事の比較は、あくまで「観測値の経験分布」に対する記述的なものである。
A.3 視差距離のバイアスと選択効果
視差サンプルには、本文の結果と同じ向きに働きうるバイアスが少なくとも2つある。
第一に Lutz–Kelker バイアス（Lutz & Kelker 1973; Verbiest et al. 2010, 2012）である。視差の測定誤差と天体の空間分布の効果により、観測された視差は系統的に過大になりやすく、$d = 1/\pi$ による距離は過小評価、したがって $v_\perp = 4.74,\mu d$ も過小評価されうる。つまり「低速側に寄って見える」という本文の傾向の一部は、このバイアスで説明される余地がある。S/N $>3$ のカットでもこの効果は完全には消えない。
第二に選択効果である。視差が測定できるのは基本的に近傍（高々数 kpc）の明るいパルサーに限られ、young PX sample は近傍天体に偏ったサブサンプルである。さらに、高速度の天体は誕生後すみやかに銀河面から離れるため、面内のサーベイで検出されにくくなる方向の選択も原理的にはありうる。これらが速度分布に対して中立である保証はない。視差サンプルの速度分布を真のキック分布として解釈するには、これらの補正を陽に扱った解析（例えば Disberg & Mandel 2025）との比較が必要である。
A.4 特性年齢の限界と Shklovskii 効果
サンプル分類に用いた特性年齢 $\tau = P/2\dot{P}$ は、制動指数 $n=3$（磁気双極子放射）、初期周期 $P_0 \ll P$、磁場一定という仮定に依存する推定量であり、超新星残骸との対応がつく天体では真の年齢と数倍ずれる例が知られている。本記事の $\tau < 10,\mathrm{Myr}$ という閾値は、厳密な年齢カットではなく「リサイクルを受けていない通常パルサー」を選ぶための操作的な基準として読むべきである。閾値を変えても結果が大きく変わらないことは第7節で確認した。
また、観測される $\dot{P}$ には固有運動に由来する見かけの寄与（Shklovskii 効果; Shklovskii 1970）
$$
\frac{\dot{P}_{\mathrm{Shk}}}{P} = \frac{\mu^2 d}{c}
$$
が含まれる。本サンプルは固有運動が測定された（しばしば大きい）天体に偏っているため原理的には注意が必要だが、若い通常パルサーでは固有スピンダウンの $\dot{P}$ が十分大きく、この補正が特性年齢の分類を変えることは実質的にない。この効果が深刻になるのは主にミリ秒パルサーであり、それらは本サンプルから除外済みである。
A.5 モデル比較の解釈について
第一に、$N=45$ と小サンプルなので、有限サンプル補正をした AICc $= \mathrm{AIC} + 2k(k+1)/(N-k-1)$ を使う方が丁寧である。ただし log-normal（$k=2$）への補正は $+0.29$ 程度で、$\Delta\mathrm{AIC} \approx 1$ という本文の結論は変わらない。そもそも経験則として $\Delta\mathrm{AIC} < 2$ は実質的な差なしとされる範囲であり、本文の表は「log-normal が勝つ」ではなく「2つの自由パラメータを使っても Hobbs 型と同等以上」程度に読むのが正しい。
第二に、本記事は1成分の Rayleigh と log-normal だけを比較したが、文献では低速・高速の2成分 Maxwellian でキック分布を記述する研究系列がある。本記事で見えた「低速側への寄り」は、log-normal でなくとも、低速成分を持つ2成分モデルで説明できる可能性がある。モデル選択の候補を2つに絞ったのは簡単のためであり、分布の真の形を制約する意図ではない。
A.6 コード
"""
Pulsar kick velocity analysis — standalone reproduction script
==============================================================
Fetches the ATNF Pulsar Catalogue via psrqpy, computes transverse kick
velocities for young isolated pulsars, and compares the Rayleigh distribution
of Hobbs et al. (2005) with a log-normal MLE fit.
Outputs
-------
fig_histogram.png velocity histogram + model PDFs
fig_cdf.png empirical CDF vs model CDFs with KS statistics
stdout cut-flow table, MLE parameters, KS/AIC summary
Tested with
-----------
psrqpy 1.3.2 / numpy 2.4.4 / scipy 1.17.1 / matplotlib 3.10.9 / Python 3.14.5
Note: the first run downloads ~5–10 MB from the ATNF server and caches it
locally (psrqpy default cache location). Results depend on catalogue version;
numbers in the article are based on v2.8.0.
References
----------
Hobbs, G. et al. 2005, MNRAS, 360, 974
Disberg, P. & Mandel, I. 2025, ApJL (arXiv:2505.22102)
Manchester, R. N. et al. 2005, AJ, 129, 1993
"""
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import psrqpy
# ── Constants ─────────────────────────────────────────────────────────────
C_VT = 4.74047 # 1 AU/yr → km/s
YR_S = 365.25 * 86400
SIGMA_HOBBS = 265.0 # Hobbs+05 1-D velocity dispersion [km/s]
AGE_THRESH_YR = 1e7 # characteristic age cut: 10 Myr
P0_MIN_S = 0.030 # exclude recycled pulsars below this period [s]
PX_SNR_MIN = 3.0 # parallax S/N threshold (ω / σ_ω)
N_BOOT = 10_000 # parametric bootstrap iterations
SEED = 42 # random seed (fixed for reproducibility)
C_DATA = "#2C3E50"
C_HOBBS = "#C0392B"
C_LN = "#1F618D"
# ── 1. Fetch catalogue ────────────────────────────────────────────────────
print("Fetching ATNF Pulsar Catalogue...")
q = psrqpy.QueryATNF(
params=["NAME", "PMRA", "PMDEC", "PX", "P0", "P1",
"AGE", "TYPE", "BINARY", "ASSOC"],
cache=True, # PX_ERR is added automatically by psrqpy
)
df = q.pandas.copy()
print(f" Catalogue version : {q.get_version}")
print(f" Total pulsars : {len(df)}\n")
# ── 2. Sample selection ───────────────────────────────────────────────────
n_all = len(df)
# proper-motion measurements in both components; remove (0, 0) placeholders
df = df[df["PMRA"].notna() & df["PMDEC"].notna()]
df = df[~((df["PMRA"] == 0.0) & (df["PMDEC"] == 0.0))]
n_pm = len(df)
df = df[df["BINARY"].isna()] # isolated only
n_binary = len(df)
df = df[df["P0"].isna() | (df["P0"] >= P0_MIN_S)] # exclude MSPs
n_msp = len(df)
df = df[~(df["ASSOC"].fillna("") == "GC")] # no globular clusters
n_gc = len(df)
df = df[~df["TYPE"].fillna("").isin(["AXP", "XINS"])] # no AXP / XINS
n_type = len(df)
# parallax sub-sample (PX_ERR auto-added by psrqpy when PX is queried)
px_snr = df["PX"] / df["PX_ERR"]
has_px = (df["PX"] > 0) & df["PX_ERR"].notna() & (px_snr > PX_SNR_MIN)
df_px = df[has_px].copy()
n_px = len(df_px)
# characteristic age τ_c < 10 Myr; fall back to P0/(2P1) if AGE is missing
mask_calc = df_px["AGE"].isna() & df_px["P1"].notna() & (df_px["P1"] > 0)
df_px.loc[mask_calc, "AGE"] = (
df_px.loc[mask_calc, "P0"] / (2.0 * df_px.loc[mask_calc, "P1"]) / YR_S
)
df_young = df_px[df_px["AGE"].notna() & (df_px["AGE"] < AGE_THRESH_YR)].copy()
n_young = len(df_young)
print("=== Cut-flow ===")
print(f" All pulsars : {n_all:5d}")
print(f" Proper motion (both components) : {n_pm:5d}")
print(f" Isolated (no BINARY flag) : {n_binary:5d}")
print(f" Non-recycled (P0 >= 30 ms) : {n_msp:5d}")
print(f" No globular cluster (ASSOC != GC) : {n_gc:5d}")
print(f" No AXP / XINS : {n_type:5d}")
print(f" Parallax S/N > 3 : {n_px:5d}")
print(f" tau_c < 10 Myr [young PX sample] : {n_young:5d}\n")
# ── 3. Transverse velocity ────────────────────────────────────────────────
# d [kpc] = 1 / parallax [mas]; v_perp = C_VT * mu [mas/yr] * d [kpc]
df_young["DIST_PX"] = 1.0 / df_young["PX"]
df_young["MU_TOTAL"] = np.sqrt(df_young["PMRA"]**2 + df_young["PMDEC"]**2)
df_young["VT"] = C_VT * df_young["MU_TOTAL"] * df_young["DIST_PX"]
vt = df_young["VT"].dropna().values
N = len(vt)
print("=== Velocity statistics (young PX sample) ===")
print(f" N = {N}")
print(f" Median v_perp = {np.median(vt):.0f} km/s")
print(f" Mean v_perp = {np.mean(vt):.0f} km/s")
print(f" Min / Max = {vt.min():.0f} / {vt.max():.0f} km/s\n")
# ── 4. Statistical analysis ───────────────────────────────────────────────
# Log-normal MLE: lognorm.fit(data, floc=0) returns (s, 0, scale)
# mu_ln = log(scale), sigma_ln = s
# mode (peak of linear-bin PDF) = scale * exp(-s^2)
s_ln, _, scale_ln = stats.lognorm.fit(vt, floc=0)
mu_ln = np.log(scale_ln)
mode_ln = scale_ln * np.exp(-s_ln**2)
print("=== Log-normal MLE ===")
print(f" mu_ln = {mu_ln:.3f}")
print(f" sigma_ln = {s_ln:.3f}")
print(f" mode = {mode_ln:.0f} km/s")
print(f" median = {scale_ln:.0f} km/s\n")
# KS tests (one-sample, two-sided)
ks_ray = stats.kstest(vt, lambda x: stats.rayleigh.cdf(x, scale=SIGMA_HOBBS))
ks_ln = stats.kstest(vt, lambda x: stats.lognorm.cdf(x, s_ln, loc=0, scale=scale_ln))
D_obs = ks_ln.statistic
# Parametric bootstrap p-value for the log-normal KS test.
# The standard KS p-value is conservative when parameters are estimated from
# the same data. We resample N_BOOT synthetic datasets from the fitted
# log-normal, refit each, compute D*, and report P(D* >= D_obs) as p_boot.
print(f"Running parametric bootstrap (N_boot={N_BOOT}, seed={SEED})...",
end="", flush=True)
rng = np.random.default_rng(SEED)
D_boot = np.empty(N_BOOT)
for i in range(N_BOOT):
x_b = rng.lognormal(mean=mu_ln, sigma=s_ln, size=N)
s_b, _, sc_b = stats.lognorm.fit(x_b, floc=0)
D_boot[i] = stats.kstest(
x_b,
lambda x, s=s_b, sc=sc_b: stats.lognorm.cdf(x, s, loc=0, scale=sc),
).statistic
p_boot = float(np.mean(D_boot >= D_obs))
print(" done\n")
# AIC: Rayleigh uses fixed sigma (k=0); log-normal fits mu and sigma (k=2)
ll_ray = float(np.sum(stats.rayleigh.logpdf(vt, scale=SIGMA_HOBBS)))
ll_ln = float(np.sum(stats.lognorm.logpdf(vt, s_ln, loc=0, scale=scale_ln)))
aic_ray = 0 * 2 - 2 * ll_ray
aic_ln = 2 * 2 - 2 * ll_ln
daic = aic_ray - aic_ln
print("=" * 63)
print(f"{'Model':<26s} {'D':>6s} {'KS p':>7s} {'AIC':>8s} {'k':>2s}")
print("-" * 63)
print(f"{'Rayleigh sigma=265 km/s':<26s} "
f"{ks_ray.statistic:>6.3f} {ks_ray.pvalue:>7.3f} {aic_ray:>8.1f} 0")
print(f"{'Log-normal MLE':<26s} "
f"{ks_ln.statistic:>6.3f} {ks_ln.pvalue:>7.3f} {aic_ln:>8.1f} 2")
print(f"{' (bootstrap p)':<26s} {'':>6s} {p_boot:>7.3f}")
print("-" * 63)
print(f"Delta AIC = AIC_Rayleigh - AIC_lognormal = {daic:+.1f}")
print("=" * 63)
print()
# ── 5. Figures ────────────────────────────────────────────────────────────
plt.rcParams.update({
"font.family" : "serif",
"font.size" : 11,
"axes.labelsize" : 12,
"xtick.labelsize" : 10,
"ytick.labelsize" : 10,
"legend.fontsize" : 9.5,
"legend.framealpha" : 0.92,
"legend.edgecolor" : "#AAAAAA",
"figure.dpi" : 150,
"savefig.dpi" : 300,
"savefig.bbox" : "tight",
"axes.spines.top" : False,
"axes.spines.right" : False,
"axes.grid" : False,
"lines.linewidth" : 1.8,
"xtick.direction" : "in",
"ytick.direction" : "in",
"xtick.minor.visible": True,
"ytick.minor.visible": True,
})
v_grid = np.linspace(0.5, 1200, 3000)
pdf_ray = stats.rayleigh.pdf(v_grid, scale=SIGMA_HOBBS)
pdf_ln = stats.lognorm.pdf(v_grid, s_ln, loc=0, scale=scale_ln)
# Histogram + model PDFs
fig, ax = plt.subplots(figsize=(6, 4.5))
bins = np.linspace(0, 1000, 20)
ax.hist(vt, bins=bins, density=True,
histtype="stepfilled", color=C_DATA, alpha=0.28, linewidth=0)
ax.hist(vt, bins=bins, density=True,
histtype="step", color=C_DATA, linewidth=1.0,
label=rf"Young PX ($N={N}$, $\tau_c < 10\,\mathrm{{Myr}}$)")
ax.plot(v_grid, pdf_ray, color=C_HOBBS, lw=1.8, ls="--",
label=rf"Hobbs+05 Rayleigh ($\sigma=265\,\mathrm{{km\,s^{{-1}}}}$)")
ax.plot(v_grid, pdf_ln, color=C_LN, lw=1.8,
label=rf"Log-normal MLE (mode $={mode_ln:.0f}\,\mathrm{{km\,s^{{-1}}}}$)")
ax.set_xlabel(r"$v_\perp\;(\mathrm{km\,s^{-1}})$")
ax.set_ylabel(r"Probability density $(\mathrm{km^{-1}\,s})$")
ax.set_xlim(0, 1000)
ax.set_ylim(bottom=0)
ax.legend(loc="upper right", handlelength=1.8)
ax.xaxis.set_minor_locator(ticker.MultipleLocator(50))
ax.tick_params(which="both", direction="in")
fig.savefig("fig_histogram.png")
plt.close()
print("Saved: fig_histogram.png")
# CDF comparison
fig, ax = plt.subplots(figsize=(6, 4.8))
vt_sorted = np.sort(vt)
ecdf_y = np.arange(1, N + 1) / N
ax.step(vt_sorted, ecdf_y, where="post",
color=C_DATA, lw=1.8, label=rf"Young PX ($N={N}$)")
ax.plot(v_grid, stats.rayleigh.cdf(v_grid, scale=SIGMA_HOBBS),
color=C_HOBBS, lw=1.7, ls="--", label=r"Hobbs+05 Rayleigh")
ax.plot(v_grid, stats.lognorm.cdf(v_grid, s_ln, loc=0, scale=scale_ln),
color=C_LN, lw=1.7, ls="-.", label=r"Log-normal MLE")
ax.text(0.03, 0.97,
rf"Rayleigh: $D={ks_ray.statistic:.3f}$, $p={ks_ray.pvalue:.2f}$",
transform=ax.transAxes, ha="left", va="top",
fontsize=8.5, color=C_HOBBS)
ax.text(0.03, 0.88,
rf"Log-norm: $D={ks_ln.statistic:.3f}$, $p={ks_ln.pvalue:.2f}$"
rf" ($p_{{\mathrm{{boot}}}}={p_boot:.2f}$)",
transform=ax.transAxes, ha="left", va="top",
fontsize=8.5, color=C_LN)
ax.axhline(0.5, color="#AAAAAA", lw=0.8, ls=":")
ax.set_xlabel(r"$v_\perp\;(\mathrm{km\,s^{-1}})$")
ax.set_ylabel("Cumulative fraction")
ax.set_xlim(0, 1050)
ax.set_ylim(0, 1.02)
ax.legend(loc="lower right", handlelength=2.0,
handletextpad=0.5, labelspacing=0.35)
ax.xaxis.set_minor_locator(ticker.MultipleLocator(100))
ax.tick_params(which="both", direction="in")
fig.savefig("fig_cdf.png")
plt.close()
print("Saved: fig_cdf.png")