天体の構造を記述する方程式

連続の式と静水圧平衡の式

　まずは天体内部の構造と質量の関係について立式します。 球対称な構造をしているガス球を考え、動径方向の座標を \(r\) とします。 ガス密度は半径の関数とし、\(\rho\!\left(r\right)\) と表します。 このとき、ある半径 \(r\) から内側にいる流体の全質量を \(\Mr\) と置くと、これは中心から \(r\) までの積分を使って \begin{align*} \Mr=\int_{0}^{r}4\pi r^{\prime2}\rho\!\left(r\right)dr^{\prime} \end{align*} と書くことができます。 この式は、半径が \(r^{\prime}\) で厚みが \(dr^{\prime}\) の薄い球殻の質量を 0 から \(r\) まで積分していることに対応しています。 つまり、中心から半径 \(r\) より内側に含んでいる質量を意味していて、enclosed mass と呼ばれることがあります。 日本語だと「内包質量」という表記がされる場合があります。 この式の両辺を \(r\) で微分して微分形で書き表すと、 \begin{align*} \frac{d\Mr}{dr}=4\pi r^{2}\rho\!\left(r\right) \end{align*} となります。 この式は連続の式 (continuity equation) の一種です。

　ここで、ポリトロープ関係と呼ばれる関係を導入します。 これは、密度と圧力の間に \begin{align*} p=K\rho^{1+\frac{1}{N}} \end{align*} という関係があるとするものです。 ここで \(K\) と \(N\) は定数のパラメータであり、\(N\) のことをポリトロープ指数 (polytropic index) と呼びます。 ポリトロープ指数は天体の構造を決める重要なパラメータとなります。 この関係式で表される流体の球のことをポリトロープ球と呼びます。

レーン=エムデン方程式の導出

　まずは式 \eqref{eq1} と式 \eqref{eq2} から、\(\Mr\) を消去します。 式 \eqref{eq2} を式 \eqref{eq1} に代入すると、 \begin{align} \frac{d}{dr}\left(\frac{r^{2}}{\rho}\frac{dp}{dr}\right)=-4\pi Gr^{2}\rho \label{eq4} \end{align} となります。

　一方、\eqref{eq3} の両辺を \(r\) で微分すると \begin{align*} \frac{dp}{dr}=K\left(1+\frac{1}{N}\right)\rho^{\frac{1}{N}}\frac{d\rho}{dr} \end{align*} が得られます。 これを式 \eqref{eq4} に代入すると以下のようになります。 \begin{align} K\left(1+\frac{1}{N}\right)\frac{d}{dr}\left(r^{2}\rho^{\frac{1}{N}-1}\frac{d\rho}{dr}\right)=-4\pi Gr^{2}\rho \tag{4'} \label{eq4p} \end{align}

　ここで、無次元化した密度を導入します。 中心部、つまり \(r=0\) での密度を \(\rhoc\) と置き、無次元化した密度 \(\theta\) をポリトロープ指数を使って \begin{align} \rho=\rhoc\theta^{N} \label{eq5} \end{align} と置きます。 これを式 \eqref{eq3} に代入すると \begin{align*} p=K\rhoc^{1+\frac{1}{N}}\theta^{N\left(1+\frac{1}{N}\right)}=K\rhoc^{1+\frac{1}{N}}\theta^{N+1} \end{align*} となります。 一方、同じく式 \eqref{eq3} において、中心での密度 \(\rhoc\) と中心での圧力 \(\pc\) を用いると \begin{align} \pc=K\rhoc^{1+\frac{1}{N}} \tag{3'} \label{eq3p} \end{align} であることが分かるため、圧力は無次元化した密度とポリトロープ指数を使って \begin{align} p=\pc\theta^{N+1} \label{eq6} \end{align} と書くことができます。

　無次元化した密度と圧力を用いて、式 \eqref{eq4p} を \(\theta\) の式に置き換えます。 式 \eqref{eq5} の両辺を \(r\) で微分すると、 \begin{align*} \frac{d\rho}{dr}=\frac{d}{dr}\left(\rhoc\theta^{N}\right)=\rhoc N\theta^{N-1}\frac{d\theta}{dr} \end{align*} となります。 また、 \begin{align*} \rho^{\frac{1}{N}-1}=\rhoc^{\frac{1}{N}-1}\theta^{1-N} \end{align*} です。 これらを用いると式 \eqref{eq4p} は \begin{align*} K\left(1+\frac{1}{N}\right)\frac{d}{dr}\left(r^{2}\rhoc^{\frac{1}{N}-1}\theta^{1-N}\rhoc N\theta^{N-1}\frac{d\theta}{dr}\right)=-4\pi Gr^{2}\rhoc\theta^{N} \tag{4''} \label{eq4pp} \end{align*} となり、整理すると \begin{align*} \frac{K\left(N+1\right)\rhoc^{\frac{1}{N}-1}}{4\pi G}\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{d\theta}{dr}\right)=-\theta^{N} \end{align*} となります。 あるいは、式 \eqref{eq3p} を用いて \(K\) を消去すると \begin{align} \frac{\left(N+1\right)\pc}{4\pi G\rhoc^{2}}\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{d\theta}{dr}\right)=-\theta^{N} \tag{4'''} \label{eq4ppp} \end{align} と書くこともできます。

　ここで、左辺の係数部分を \begin{align*} \alpha^{2}=\frac{\left(N+1\right)\pc}{4\pi G\rhoc^{2}} \end{align*} と置きます。 これは長さの 2 乗の次元を持つ (つまり \(\alpha\) は長さの次元を持つ) ことがすぐに確かめられます。 これを用いると \begin{align} \frac{\alpha^{2}}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{d\theta}{dr}\right)=-\theta^{N} \tag{4''''} \label{eq4pppp} \end{align} となります。

　ここまでで密度と圧力は無次元化してありましたが、ここで \(\alpha\) を用いて \(r\) も無次元化します。 無次元化した半径として、\(\xi\) を以下のように定義します。 \begin{align*} \xi=\frac{r}{\alpha} \end{align*} 無次元化した密度 \(\theta\) の \(r\) 微分は、次のように書き換えられます。 \begin{align*} \frac{d\theta}{dr}=\frac{d\theta}{d\xi}\frac{d\xi}{dr}=\frac{1}{\alpha}\frac{d\theta}{d\xi} \end{align*} またこの式から、微分演算子は \begin{align*} \frac{d}{dr}=\frac{1}{\alpha}\frac{d}{d\xi} \end{align*} となることも分かります。 これらを用いると、式 \eqref{eq4pppp} は以下のようにまとめることができます。 \begin{align*} \frac{1}{\xi^{2}}\frac{d}{d\xi}\left(\xi^{2}\frac{d\theta}{d\xi}\right)=-\theta^{N} \end{align*} この式は、無次元化した密度 \(\theta\)、無次元化した半径 \(\xi\) を結びつける式になっています。 従って、この式を解くことができれば密度と半径の関係が分かり、天体の構造を決定出来ることになります。

　この式はレーン=エムデン方程式 (Lane-Emden equation) として知られています。 また、この式の解はしばしばエムデン解と呼ばれます。

境界条件

　レーン=エムデン方程式は \(\xi\) の 2 階微分方程式になっているため、解を決定する為には境界条件が 2 つ必要です。 境界条件としては原点での値を与えます。

　まずは、式 \eqref{eq5} にあるように密度は中心で \(\rhoc\) となります。 つまり、中心では \(\theta=1\) となっている必要があります。 \(r=0\) は \(\xi=0\) に対応しているため、1 つ目の境界条件は 「\(\xi=0\) のとき \(\theta=1\)」 となることが分かります。

　また、中心では重力は 0 となり、圧力勾配力も 0 となります。 圧力勾配力が 0 であるということは、圧力勾配が 0 であることを意味しています。 密度と圧力は式 \eqref{eq3} のポリトロープ関係で結びつけられているので、圧力勾配が 0 であれば密度勾配も 0 になります。 当然、無次元化した密度の勾配も 0 となります。 そのため、2 つ目の境界条件は 「\(\xi=0\) のとき \(\displaystyle{\frac{d\theta}{d\xi}=0}\)」 となります。

方程式の解

　レーン=エムデン方程式については、\(N=0,\,1,\,5\) の時は解析解が存在することが知られています。 それぞれ、

(i) \(N=0\) のとき \begin{align*} \theta\left(\xi\right)=1-\frac{1}{6}\xi^{2} \end{align*}

(ii) \(N=1\) のとき \begin{align*} \theta\left(\xi\right)=\frac{\sin\xi}{\xi} \end{align*}

(iii) \(N=5\) のとき \begin{align*} \theta\left(\xi\right)=\frac{1}{\sqrt{1+\frac{1}{3}\xi^{2}}} \end{align*}

となります。

　これら以外の一般の \(N\) に対しては、数値計算をして構造を決定する必要があります。

　また、\(\xi\) は半径、\(\theta\) は密度に対応しているものであるため、それぞれ 0 以上の値を持ちます。 \(\theta\) は中心では 1 であり、そこから外側へ向かって減少して行く形になっていることが分かります。 \(\theta\) が 0 になる場所は密度と圧力が 0 になる距離であり、天体の "表面" に対応しています。 そこより先も関数は値を持ちますが、物理的に意味があるのは最初に \(\theta\) が 0 になる距離までです。