円柱状自己重力流体の平衡形状の導出

基礎方程式

　重力と圧力勾配力が釣り合っている静水圧平衡状態を考えます。 静水圧状態の式は、 \begin{align} \nabla\Phi = -\frac{1}{\rho}\nabla p \label{hseq} \end{align} と書けます。 ここで \(\Phi\) は重力ポテンシャル、\(\rho\) は流体の密度、\(p\) は流体のガス圧です。 ここでは円筒座標 \(\left(r,\,\phi,\,z\right)\) を考えますが、今回は一様で軸対称の構造を仮定します。 つまり、円筒の長さ方向 (\(z\) 方向) に一様で無限に長く、軸の周りの回転 (\(\phi\) 方向) に対して対称であるとします。 そのため、密度や圧力などは、円柱の半径方向 (\(r\) 方向) のみの関数になります。

　また、重力ポテンシャルと流体の密度分布は、以下のポアソン方程式で表されます。 \begin{align} \nabla^{2}\Phi = 4\pi G\rho \label{poisson} \end{align} \(G\) は重力定数 (万有引力定数) です。

　さらに、流体が等温の理想気体であると仮定します。 この場合、状態方程式は次のように書けます。 \begin{align} p = \frac{k_{\rm B}T}{\mu m}\rho \equiv K\rho \label{eos} \end{align} ここで、\(k_{\rm B}\) はボルツマン定数、\(T\) は温度、\(\mu\) は平均分子量、\(m\) は陽子の質量です。 なお、ここで等温音速を \(\cs\) を置くと \begin{align} \cs^{2} = \frac{k_{\rm B}T}{\mu m} \label{kcs} \end{align} であるため、式 \eqref{eos} は \begin{align*} p=\cs^{2}\rho \end{align*} と書けます。

　ここで、式 \eqref{hseq} と \eqref{eos} から \(p\) を消去すると、 \begin{align*} \nabla\Phi= -\frac{1}{\rho}\nabla K\rho = -\nabla K\ln\rho \end{align*} が得られます。 この式の両辺の発散をとる (両辺に \(\nabla\) をかける) と、 \begin{align*} \nabla^{2}\Phi = -\nabla^{2}K\ln\rho \end{align*} となります。 これに式 \eqref{poisson} を用いて \(\Phi\) を消去すると、 \begin{align*} -\nabla^{2}K\ln\rho = 4\pi G\rho \end{align*} となります。

　まずは以下のような変数変換を行います。

　これらの関係式を用いて、式 \eqref{eq} を \(\xi\) と \(\psi\) の式に直していきます。 式 \eqref{psi} を用いると、式 \eqref{eq} の左辺は \begin{align*} K\nabla^{2}\ln\rho &= K\nabla^{2}\ln\left(\rho_{0}e^{-\psi}\right) \\ &= K\nabla^{2}\left(\ln\rho_{0}+\ln e^{-\psi}\right)\\ &=K\nabla^{2}\left(\ln\rho_{0}-\psi\right)\\ &=-K\nabla^{2}\psi \end{align*} となります。 \(\rho_{0}\) は定数なので、ラプラシアン \(\left(\Delta=\nabla^{2}\right)\) を作用させると消えてしまいます。 よって、式 \eqref{eq} は \begin{align*} K\nabla^{2}\psi = 4\pi G\rho_{0}e^{-\psi} \end{align*} となります。 ここで式 \eqref{alpha} を用いることで、この式は以下のように書き直すことが出来ます。 \begin{align} \alpha^{2}\nabla^{2}\psi = e^{-\psi} \label{eq2} \end{align}

　ここで、軸対称の円筒座標系であることを考えると、微分演算子は \begin{align*} \nabla^{2} = \frac{1}{r}\frac{d}{dr}\left(r\frac{d}{dr}\right) \end{align*} となります。 また、式 \eqref{xi} から、 \begin{align*} dr=\alpha d\xi \end{align*} であることが分かるため、 \begin{align*} \frac{d\xi}{dr}=\frac{1}{\alpha} \end{align*} となることが分かります。 そのため、微分演算子には \begin{align*} \frac{d}{dr} = \frac{d}{d\xi}\frac{d\xi}{dr} = \frac{1}{\alpha}\frac{d}{d\xi} \end{align*} という関係があることも分かります。

　これらの関係を用いると、 \begin{align*} \nabla^{2}\psi & = \frac{1}{r}\frac{d}{dr}\left(r\frac{d}{dr}\psi\right)\\ & = \frac{1}{\alpha\xi}\frac{1}{\alpha}\frac{d}{d\xi}\left(\alpha\xi\frac{1}{\alpha}\frac{d}{d\xi}\psi\right)\\ & = \frac{1}{\alpha^{2}\xi}\frac{d}{d\xi}\left(\xi\frac{d\psi}{d\xi}\right)\\ & = \frac{1}{\alpha^{2}\xi}\left(\frac{d\psi}{d\xi}+\xi\frac{d^{2}\psi}{d\xi^{2}}\right)\\ & = \frac{1}{\alpha^{2}}\left(\frac{1}{\xi}\psi'+\psi''\right) \end{align*} と変形できます。 なお、\(\xi\) の 1 階微分を \(^\prime\)、2 階微分を\(^{\prime\prime}\) を付けて表しています。 これを用いると、式 \eqref{eq2} は以下のようになります。 \begin{align} \psi'' + \frac{1}{\xi}\psi' = e^{-\psi} \label{psieq} \end{align}

　式 \eqref{xi} より、\(r=0\) のときは \(\xi=0\) となります。 また式 \eqref{psi} で \(r=0\) の時に \(\rho=\rho_{0}\) であるためには、\(r=0\) すなわち \(\xi=0\) のときに \(\psi=0\) である必要があります。 さらに、式 \eqref{psieq} は \begin{align*} \xi\psi'' + \psi' = \xi e^{-\psi} \end{align*} と書け、\(\xi=0\) のときは \begin{align*} 0 + \psi' = 0 \end{align*} となるため \(\psi'=0\) を満たす必要があります。

　式 \eqref{psieq} を解くため、ここからさらに以下のような変数変換を行います。

　式 \eqref{t} より、\(dt = \sqrt{2}\displaystyle{\frac{d\xi}{\xi}}\) となるため、 \(\displaystyle{\frac{dt}{d\xi} = \frac{\sqrt{2}}{\xi}}\) となることが分かります。 また、\(\ln\xi=\displaystyle{\frac{t}{\sqrt{2}}}\) であるため、 \(\displaystyle{\xi=e^{t/\sqrt{2}}}\) となります。 これらの関係式から、 \begin{align*} \frac{dt}{d\xi} = \sqrt{2}\expt \end{align*} であることが分かります。

　従って、微分演算子の変換は以下のようになります。 \begin{align*} \frac{d}{d\xi}=\frac{d}{dt}\frac{dt}{d\xi}=\sqrt{2}\expt\frac{d}{dt} \end{align*} これらを用いて、式 \eqref{psieq} を変形します。

　次に、式 \eqref{zeq} を積分します。 \(e^{z}\) の項を右辺に移した \begin{align*} 2\frac{d^{2}z}{dt^{2}} = -e^{z} \end{align*} を考えます。 両辺に \(\displaystyle{\frac{dz}{dt}}\) をかけると \begin{align*} 2\frac{dz}{dt}\frac{d^{2}z}{dt^{2}} = -\frac{dz}{dt}e^{z} \end{align*} となります。 この式の両辺を積分すると、 \begin{align*} \left(\frac{dz}{dt}\right)^{2}=-e^{z}+C \end{align*} となります。 ここで \(C\) は積分定数です。 そのため、\(\displaystyle{\frac{dz}{dt}}\) についての以下の式が得られます。 \begin{align} \frac{dz}{dt} = \pm\left(C-e^{z}\right)^{1/2} \label{zinteg} \end{align}

　ここで、境界条件を用いて積分定数を決定します。 \(\displaystyle{\frac{dz}{dt}}\) を \(\xi\) および \(\psi\) で表すと、以下のようになります。 \begin{align*} \frac{dz}{dt}=\frac{dz}{d\xi}\frac{d\xi}{dt}&=\frac{1}{\sqrt{2}}e^{t/\sqrt{2}}\frac{d}{d\xi}\left(-\psi+2\ln\xi\right)\\ &= \frac{1}{\sqrt{2}}\xi\left(-\psi'+\frac{2}{\xi}\right)\\ &= -\frac{1}{\sqrt{2}}\xi\psi'+\sqrt{2} \end{align*} 従って、\(r=0\)、つまり \(\xi=0\) のときは、 \begin{align*} \left.\frac{dz}{dt}\right|_{\xi=0}=\left(-\frac{1}{\sqrt{2}}\xi\psi'+\sqrt{2}\right)_{\xi=0}=\sqrt{2} \end{align*} となります。

　式 \eqref{dzdt} を解くため、さらに変数変換を行います。

　式 \eqref{yinteg0} の両辺を積分すると、 \begin{align} \int\frac{dy}{y\left(2-y\right)^{1/2}}=\pm\int dt \label{yinteg} \end{align} となります。

　最後に、変数を元に戻します。 式 \eqref{psi} より、 \begin{align*} \psi=-\ln\frac{\rho}{\rho_{0}} \end{align*} であることが分かります。 これを式 \eqref{psisol} に代入することで、 \begin{align*} \ln\frac{\rho}{\rho_{0}} = \ln\left(1+\frac{1}{8}\xi^{2}\right)^{-2} \end{align*} が得られるため、 \begin{align} \rho=\rho_{0}\frac{1}{\left(1+\frac{1}{8}\xi^{2}\right)^{2}} \label{rhoxi} \end{align} となります。

　また、式 \eqref{xi} と \eqref{alpha} より、 \begin{align*} \xi^{2}=\frac{1}{\alpha^{2}}r^{2}=\frac{4\pi G\rho_{0}}{K}r^{2} \end{align*} となります。 ここで式 \eqref{kcs} より \(K\) を \(\cs^{2}\) に置き換え、さらに \begin{align} H_{0}^{2}\equiv\frac{2\cs^{2}}{\pi G\rho_{0}} \label{h0} \end{align} と置くと、 \begin{align*} \frac{1}{8}\xi^{2}=\frac{\pi G\rho_{0}}{2\cs^{2}}r^{2}=\frac{r^{2}}{H_{0}^{2}} \end{align*} となります。

　これを式 \eqref{rhoxi} に代入して整理すると、以下の式が得られます。 \begin{align} \rho\!\left(r\right)=\rho_{0}\left[1+\left(\frac{r}{H_{0}}\right)^{2}\right]^{-2} \label{cylinder} \end{align}

線密度の計算

　式 \eqref{cylinder} で計算した密度分布を元に、線密度 (line density) を計算します。 これは、線や円柱などの細長い物体の、単位長さあたりの質量を意味します。 密度にもいろいろな種類がありますが、一般に密度と言った場合は体積密度 (volume density) を意味することが多く、単位は SI であれば \({\rm kg\,m^{-3}}\) です。 これは単位体積あたりの質量を意味します。 その他には、単位面積あたりの質量である表面密度 (surface density) があり、こちらの単位は \({\rm kg\,m^{-2}}\) になります。 表面密度は面密度と言われることもあります。 今回求める線密度は、単位は \({\rm kg\,m^{-1}}\) になります。

　線密度 \(\ML\) は、以下のように計算します。 \begin{align*} \ML=\int_{0}^{\infty}2\pi r\rho\!\left(r\right)dr \end{align*} これに式 \eqref{cylinder} を代入すると、 \begin{align} \ML=2\pi\rho_{0}\int_{0}^{\infty}\frac{r\,dr}{\left[1+\left(\frac{r}{H_{0}}\right)^{2}\right]^{2}} \label{mlinteg} \end{align} となります。

積分計算その 1

　途中で、 \begin{align*} \int\frac{dx}{x\left(2-x\right)^{1/2}}=-\sqrt{2}\tanh^{-1}\left(\sqrt{1-\frac{x}{2}}\right)+{\rm const.} \end{align*} という積分計算が出てきました。 これは、以下のように変数変換を行うことで計算することが出来ます。

　まず、以下のように変形します。 \begin{align} \int\frac{dx}{x\sqrt{2-x}}=\int\frac{dx}{\sqrt{2}x\sqrt{1-\frac{x}{2}}} \label{A1} \tag{A1} \end{align} ここで、 \begin{align} t=\sqrt{1-\frac{x}{2}} \label{A2} \tag{A2} \end{align} と置きます。 この変換より、 \begin{align*} dt=\frac{1}{2}\left(1-\frac{x}{2}\right)^{-1/2}\left(-\frac{1}{2}\right)dx=-\frac{1}{4}\frac{dx}{t} \end{align*} となるため、 \begin{align*} dx=-4t\,dt \end{align*} が得られます。 これらを用いると、式 \eqref{A1} は以下のように変形できます。 \begin{align} \int\frac{dx}{\sqrt{2}x\sqrt{1-\frac{x}{2}}}=-\sqrt{2}\int\frac{dt}{1-t^{2}} \label{A3} \tag{A3} \end{align}

　この積分は、 \begin{align*} \int\frac{dt}{1-t^{2}} &= \int\frac{dt}{\left(1+t\right)\left(1-t\right)}\\ &= \frac{1}{2}\int\left(\frac{1}{1+t}+\frac{1}{1-t}\right)dt\\ &= \frac{1}{2}\left(\ln\abs{1+t}-\ln\abs{1-t}\right) + \rm{const.}\\ &= \frac{1}{2}\ln\frac{1+t}{1-t} + \rm{const.}\\ &= \tanh^{-1}t + \rm{const.} \end{align*} と計算できます。

　もう一つの積分計算、 \begin{align*} \int\frac{x\,dx}{\left(1+ax^{2}\right)^{2}}=-\frac{1}{2a^{2}x^{2}+2a}+{\rm const.} \label{B1} \tag{B1} \end{align*} についても紹介します。

　ここでは \begin{align*} t=1+ax^{2} \end{align*} という変換を行います。 この式より、 \begin{align*} dx=\frac{1}{2ax}dt \end{align*} となることが分かります。 これを用いると、積分は \begin{align*} \int\frac{x\,dx}{\left(1+ax^{2}\right)^{2}} &= \frac{1}{2a}\int\frac{dt}{t^{2}}\\ &= -\frac{1}{2a}\frac{1}{t}+\rm{const.}\\ &= -\frac{1}{2a}\frac{1}{1+ax^{2}}+\rm{const.} \end{align*} と計算できます。 従って、式 \eqref{B1} が得られます。

参考文献