日震学（星震学）データで遊んでみる

日震学とは

太陽内部を「波」を用いて調べる学問を「日震学（helioseismology）」といいます。地震波を使って地球内部の構造探査をするのと同様の手法で、まさに地震学（seismology）の太陽（helio-）版というわけです。

太陽は圧縮性の高いプラズマガスでできており、その内部では乱流熱対流によって励起された音波で満たされていると考えられます。このようにランダムに励起された音波が無数に重ね合わさった結果、おそらく全体としては定在波を形成していると思われます。これは音波の固有振動とも考えられ、その固有関数や固有振動数は太陽内部構造やパラメータによって決定されています。日震学では、逆に太陽表面の音波を観測し、その固有振動数を測ることで、（光学観測が不可能な）太陽内部情報を調べます。これは数学的には、逆問題を解くことに相当します。

日震学は既にこれまで数多くの成果を上げており（代表的なものは太陽内部の差動回転観測など）、その知見は近年星震学（asteroseismology）として一般の恒星に応用されており、星の音波モードを調べることで（従来の天文学的手法とは完全に独立に）質量・半径・年齢といった構成パラメータを推定することができる強力なツールとしても認識されています。

太陽観測データ

日震学の観測で主に用いられているのは、太陽表面（光球）の視線方向速度です。これは、光球で形成される吸収線のドップラー遷移を測定することで得られます。
下図（左）に示したのが、観測された太陽表面でのドップラー速度マップです。ここから、測定機器の統計誤差や自転・平均流成分を除去したデータが下図（右）になります。日震学では主にこのデータを解析していくこととなります。

実は、この速度マップには太陽の音波モードだけでなく超粒状斑といった対流モード成分がランダムに重なり合っています。しかしながら、これらは時間変化のスケールが大幅に異なる（音波は約５分周期 / 対流は数日周期）ため、時間方向にフーリエ変換し数分スケールの現象のみに注目することによってこれらを区別できます。

音波のパワースペクトル

下図に示したのは、ドップラー速度マップを時間・空間方向にそれぞれスペクトル変換したパワースペクトル図です。一般に波の解析によく用いられる手法で、こういったパワースペクトル図は「$ k-\omega $ ダイアグラム」と呼ばれます。ここでは、空間方向には球面調和関数変換（フーリエ変換の球面バージョン）をしていて、その際音波は等方的であると仮定して、位数 $m$ 成分については足し合わせて、次数（球面波数）$l$ のみの関数としています。時間方向には通常のフーリエ変換をしており、音波の典型的な周波数付近のみをプロットしています。

これをみると、たくさんの赤い筋状の線（これをパワーリッジと呼ぶ）が確認できます。これら一本一本が太陽内部の音波の異なるモードに対応していて、それぞれ異なる分散関係式に従って分布しています。これらのモードは、便宜上動径方向の節点の数 $n$ で分類されています。

f モード（表面重力波）

最も下に位置しているリッジは $n=0$ に対応します。実はこのモードだけ少し特殊で、太陽内部で反射されずに表面のみを伝わる波で、表面重力波（ $f$ モード）と呼ばれています。復元力は重力で、非圧縮性波動です。後述するように、太陽半径の測定にも使われています。

p モード（音波モード）

対して $n>0$ のモードは圧力傾度力を復元力とする圧縮性の音波モードで $p$ モードと呼ばれています。低波数成分（$l$ が小さい / 大スケール）に注目すると、3-4 mHz の振動数でパワーが全体的に大きくなっています。これは、周期に変換すると約５分に対応しており、「太陽の５分振動」と呼ばれています。歴史的には５分振動が先に発見されており長らく謎だったのが、後に太陽内部の音波の固有振動の集合であると理解されるようになったという経緯があります。 $p$ モードは水平波数と振動周波数に応じて太陽深部まで伝播するため、太陽内部の診断に利用されています。

公開されている太陽振動数データ

それでは実際に SDO/HMI によって測定された太陽の振動数データをダウンロードして、プロットしながら遊んでみたいと思います。データの多くは http://jsoc.stanford.edu にて公開されています。

wget http://jsoc.stanford.edu/SUM61/D664163469/S00000/m10qr.6328

上のコマンドでダウンロードしたファイルには、観測的に特定された 2000 個以上の太陽の固有モードの情報がまとめられています。これらの情報は、すべて逆問題を解く際の制約条件であるため、条件数が多ければ多いほど精度良く太陽の内部構造が推定できるようになります。ここからも日震学が如何に高精度か実感できるかと思います。
今回は、最初の３列のみ使用します。それぞれ（#1）球面水平波数 $l$（#2）動径節数 $n$（#3）振動周波数 $\nu$ に対応しています。さっそく以下の python スクリプトでプロットしてみます。

plot_freq.py

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes

#read data
data = np.loadtxt('m10qr.6328')
l = data[:,0]
n = data[:,1]
freq = data[:,2]*1.0e-3 #from micro Hz to mHz

#--- figure ---#
fig = plt.figure(figsize=(6,5))
plt.rcParams["font.size"] = 13

#plot mode frequencies
ax = fig.add_subplot(111)
ax.scatter(l,freq,s=4,c=n,cmap='jet_r',vmin=0,vmax=20)
ax.set_xlim(0,300)
ax.set_ylim(0.6,4.9)
ax.set_xlabel(r'Spherical degree $l$')
ax.set_ylabel(r'Frequency $\nu$ [mHz]')

#add a colorbar inside the plot
plt.rcParams["font.size"] = 12
axins = inset_axes(ax,width="3.5%",height="80%",loc='upper right'
                   ,bbox_to_anchor=(0.35,0.425,0.5,0.5),bbox_transform=ax.transAxes)
norm = mpl.colors.Normalize(vmin=0,vmax=20.0)
sm = plt.cm.ScalarMappable(cmap='jet_r',norm=norm)
sm.set_array([])
cb = plt.colorbar(sm,cax=axins,ticks=[0,5,10,15,20])
cb.set_label(r"$n$")

#output
plt.subplots_adjust(top=0.95,bottom=0.125,left=0.13,right=0.95)
plt.show()

各点が各固有モードに対応しています。異なる色はそれぞれ異なる $n$ のモードに対応させています。

f モード周波数を用いた太陽半径の推定

手始めに、測定された $f$ モード周波数を用いて太陽半径を計算してみます。次数 $l$ の時の伝播周波数を$\omega_{l}$ とすると、分散関係式は $\omega_{l}=\sqrt{gk_{h}}=(G M_{\odot}/R_{\odot}^{3})\sqrt{l(l+1)}$ で与えられます。重力定数 $G$ は定数なので、太陽質量 $M_{\odot}$ が既知とすると、測定値（$l$, $\omega_{l}$）から太陽半径 $R_{\odot}$ が求まります。
試しに、$l=100$ の値を使ってみると、$R_{\odot}=689$ Mm となります。国際天文学連合が定めている太陽の視半径 $R_{\odot}=696$ Mm よりも約１％小さくなってしまいましたね。では、$l=300$ としてみましょう。そうすると $R_{\odot}=693$ Mm となって標準値に近づきました。$l$ が大きいほど、半径の決定精度が上がるようです。なるほど、$l$ が十分大きければ、太陽の曲率が無視できるようになるので、表面の $f$ モードの伝播は上記の平面波の分散関係式で近似できるようになるのですね。
このように周波数解析で半径が推定できるというのは非常に面白いのですが、この手法は高い空間分解能（$l>100$）が要求されるため、他の恒星には応用できないのが残念です。

デュヴァル（ Duvall ）の法則

上図では、単純に 振動数 $\nu_{nl}$ を 次数 $l$ に対してプロットして、リッジが $n$ に応じて無数に分かれることを確認しました。それでは、（天下り的ですが）少し変数変換して、$(n+\alpha)/\nu_{nl}$ を $\nu_{nl}/\sqrt{l(l+1)}$に対してプロットするとどうなるのでしょうか。ただし $\alpha=1.5$ とします。さっそくやってみます。

plot_duvall.py

#newly-defined variables
x = freq/np.sqrt(l*(l+1.0))*1.0e3
y = (n+1.5)/freq*1.0e3

#plot
ax = fig.add_subplot(111)
ax.scatter(x,y,s=4,c=n,cmap='jet_r',vmin=0,vmax=20)
ax.set_xscale('log')
ax.set_xlabel(r'$\nu_{nl}/\sqrt{l(l+1)}$ [$\mu$Hz]')
ax.set_ylabel(r'$(n+\alpha)/\nu_{nl}$ [s]')
plt.show()

なんと！驚くべきことに、これまで $n$に応じて無数に何本にも分かれていた $p$ モードリッジが、全て $n$ に依らず１本の曲線の上に乗っています。この驚くべき観測事実は、第一発見者の名前をとって「デュヴァル（Duvall）の法則」と呼ばれています。

その後、すぐにこの法則に対する理論的裏付けがなされ、太陽内部の音波の伝播を記述する方程式に高周波数漸近近似を適用することで関係式が導出されました。簡単に物理的背景を説明すると、太陽内部は内側に行くこど高温なので音波は内部を伝播するほど強く屈折して表面に戻ってきます。実は、横軸にプロットした $\nu_{nl}/\sqrt{l(l+1)}$ は音波が屈折され反射された深さに関連した物理量になっています。一方、縦軸にプロットした $(n+\alpha)/\nu_{nl}$ はその反射地点から表面まで音波の位相を積分したような量になっており、 反射される深さが同じであれば、水平波数によらず動径積分した際の波の位相は等しくなる（節点が $n$ 個あれば位相は $n$ 倍になる）ということを示しています。

Duvall の法則は、太陽内部の音速分布を推定（これをインバージョンと呼ぶ）するのに使われる（最初に使われた / 現在ではより精度の高いインバージョン法がある）という点で非常で重要です。というのも、もし上図のプロットをフィッティングすることで曲線の関数が求まれば、残る未知数は太陽内部の音速分布だけになり、積分方程式を解くことでこれを求めることが原理的に可能になるのです。こうして日震学によって太陽内部の音速分布が測定され、内部構造（対流層の厚さなど）が観測的に決定されました。

星震学への応用（星としての太陽）

では続いて、与えられた太陽の固有振動数データを使って、もう少し踏み込んだ「星震学的な」解析をして遊んでみましょう。

低波数近似

そもそも星震学とは日震学の恒星バージョンです。しかし、残念ながら一般の星は太陽と違って細かく空間分解できません（基本的に光源 / １ピクセルからの情報）。もし様々な空間スケールの音波が励起されていたとしても、我々に届くのはそれらが全て足し合わさった情報であり、この過程で高波数成分の音波の情報はほとんど失われてしまうと考えられます。その結果、光源観測から我々が知ることができるのはせいぜい $l=0, \ 1, \ 2, \ 3$ のモードに限られます。
そこで、星震学への応用を見据えて、太陽の固有モードの低次数成分（low $l$ ）のみに注目します。理論的考察から、$n$ が比較的大きい（従って比較的高周波）領域では、近似的に $\nu_{nl}\propto (n+l/2)\Delta\nu$ が成り立つことが示されます。実際に先ほどダウンロードした HMI データで確認してみましょう。

確かにまあまあ成り立っていますね。次数 $l$ を固定したとき、振動数は $n$ に応じて $\Delta\nu$ の間隔で増えていっています。このような $\Delta\nu$ を星震学の用語では 大間隔（large separation） と呼びます。上のプロットより太陽の大間隔は約 $135.3 \ \mu$Hz と求まります。
また、$\nu_{n+1,l} \approx \nu_{n,l+2}$ であることも確認できます（厳密には成り立たない）。実際、$l=0$ と $l=2$ は１個違いの $n$ 同士で非常に近い固有振動数を持っています。これは、$l=1$ と $l=3$ のペアに対しても同様です。 星震学的観測において、スペクトル解析した際に非常に隣接した２つのピークが存在すれば、それは $l=0$ と $l=2$ のペアかもしくは $l=1$ と $l=3$ のペアのどちらかということになります（経験則として、$l=3$ のピークは $l=2$ に比べて非常に小さいので、多くの場合 $l=0$ と $l=2$ のペア）。

エシェル図（ Eshelle ダイアグラム）

さて、音波の低波数モードが大間隔おきに並んでいるのであれば、大間隔で割った余りの振動数（mod $\Delta\nu$）に対してパワースペクトルをプロットしたくなるのが、人の性というものでしょう。
では、実際にやってみましょう。以下の python スクリプトを用います。

plot_echelle.py

#measured large separation of the Sun
delnu=135.3 #[micro Hz]

#--- figure ---#
fig = plt.figure(figsize=(6,5))

plt.rcParams["font.size"] = 12
ax = fig.add_subplot(111)

#plot echelle diagram
ax.plot(np.mod(nu[np.where(l==0)],delnu),nu[np.where(l==0)],'ko')
ax.plot(np.mod(nu[np.where(l==1)],delnu),nu[np.where(l==1)],'ro')
ax.plot(np.mod(nu[np.where(l==2)],delnu),nu[np.where(l==2)],'bo')
ax.plot(np.mod(nu[np.where(l==3)],delnu),nu[np.where(l==3)],'go')
ax.set_xlim(0,delnu)
ax.set_ylim(1250,4250)
ax.set_xlabel(r'Reduced frequency mod($\Delta\nu$) [$\mu$Hz]')
ax.set_ylabel(r'Frequency $\nu$ [$\mu$Hz]')

#add text
ax.text(38,3800,r'$l=0$',color='k')
ax.text(100,4000,r'$l=1$',color='r')
ax.text(20,3800,r'$l=2$',color='b')
ax.text(80,4000,r'$l=3$',color='g')

plt.subplots_adjust(top=0.95,bottom=0.135,left=0.19,right=0.96)
plt.show()

横軸を大間隔に対する振動数剰余（Reduced frequency）、縦軸を振動数として $p$ モードをプロットしたものが下図（左）になります。それぞれの次数 $l$ に対して、さまざまな動径節点数 $n$ のモードが縦に一列に並んだような模様ができています（周波数が上がるにつれて $n$ が増加している）。また、$l=1$ と $l=3$ の２本と $l=0$ と $l=2$ の２本ラインがそれぞれ互いに隣接していることも確かめられます。
このような図を一般に星震学では、 エシェル図（Echelle diagram） と呼びます。参考のため、「星としての太陽」観測によって得られたエシェル図のパワースペクトルプロットを右図に示しました。

さて、先ほどはざっくり $\nu_{n,l} \approx \nu_{n-1,l+2}$ として説明しましたが、エシェル図を見ると $l=1$ モードと　$l=3$ モードの振動数剰余にはわずかに隔たりが存在していることに気が付きます。この差異 $\delta\nu_{n,l}=\nu_{n,l}-\nu_{n-1,l+2}$ を 小間隔（small separation） と呼びます。星震学では、エシェル図からこの小間隔を測定することが非常に重要な手続きとなっています。太陽では、$\delta\nu_{l=0}\approx 9 \ \mu$Hz、$\delta\nu_{l=1}\approx 15 \ \mu$Hz となることから、後述する $D_{0}$ パラメータは 1.5 $\mu$Hz と算出されます。

大間隔と小間隔の物理的意味

大間隔 $\Delta\nu$ は音波が星の直径を横断するのに要する時間の逆数と捉えることができます。いわば星の力学的タイムスケールを決めるパラメータで、その星の平均密度や質量と密接に関係しています。

一方、小間隔はもう少し複雑なのですが、結論から言うと星の年齢に関係した量と捉えられています。以下詳しい説明。まず、 $\delta\nu_{l} \approx (4l+6)D_{0}$ によって $l$ に依らない定数 $D_{0}$ を定義します。この定数はざっくり $D_{0}\propto -\int(dC_{s}/dr)(dr/r)$ となっていて、コア付近の音速勾配にマイナス比例した量となっています。星は誕生した時は中心部に行けば行くほど高温で至る所 $dC_{s}/dr<0$ となるので、$D_{0}>0$ となります。しかしながら、星が年齢を重ねていくにつれて中心コアでの核融合反応が進行し水素が重いヘリウムに置き換わっていくため中心の音速が遅くなり、結果的に中心部で $dC_{s}/dr>0$ となる領域が形成されます。これは $D_{0}$ を減少させる効果があるため、$D_{0}$ は星が年齢を重ねると共に小さくなります。

上図は、横軸大間隔 $\Delta\nu$、縦軸小間隔 $D_{0}$ として様々な質量の星の進化過程を示した通称 JCD 図です。今回求めた太陽パラメータ（$\Delta\nu=135 \ \mu$Hz、$D_{0}=1.5 \ \mu$Hz）は JCD 図上では赤点で示したようにちゃんと太陽（$M=M_{\odot}$）の進化トラックに乗っていますね。年齢は太陽寿命の約半分ほどでしょうか。
今回は太陽データを使いましたが、一般にこのような星震学的解析によって $\Delta\nu$ および $D_{0}$ が求まれば、JCD 図より星の質量や年齢などの基本パラメータが推定可能となるのです。星震学の威力が実感できますね。

まとめ

今回は、公開されている太陽の固有モード振動数データを使って、まず最初に日震学的観測事実（Duvall の法則等）を確かめた後、低波数成分のみを取り出すことで星震学的アプローチを試みてみました。
正直日震学は煩雑な方程式をひたすら計算するイメージが強いため、僕のように苦手意識を持っていた人も多いのではないでしょうか。授業や教科書で色々な図を見せられても、結局どのようなデータ処理をしてそうなったのかよく理解できないため、挫折するケースも多いと思います。そんな時は、実際に公開されている生のデータを使って自分で色々遊んでみることで理解が一段と深まると思います。皆さんも是非試してみてください。