原始惑星系円盤の密度進化を解く（拡散方程式の無次元化）

これは，**[湧源クラブ Advent Calendar](https://qiita.com/advent-calendar/2019/yugen) 20日目** の記事です！

**この記事を読むとできるようになること**
`一見難しそうな拡散方程式を数値計算できる`

学部の授業でやった計算を紹介します．
出てくる式は，修論の研究でもめっちゃ使っています．

## 原始惑星系円盤(Protoplanetary disk)とは
![fig1-5-1078x516.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/403437/86a10f8c-0bd5-aaa1-fe76-8f586571d494.jpeg)
credit: NASA/JPL-Caltech 

できたてほやほやの恒星の周りにできる，ガスと塵からなる円盤のことです．
質量のほとんど（99%くらい）は水素とヘリウムが主体のガスであると考えられています．
原始惑星系円盤の中で，塵やガスが集まって惑星が誕生します．
現在の太陽系には原始惑星系円盤は残っていませんが，
他の恒星の周りに存在するものは電波望遠鏡で観測されています．
今回は，この原始惑星系円盤のガスの密度進化を解いてみましょう．

## 拡散方程式
$r$方向（円盤の中心から外側に向かう向き）のガス面密度の拡散方程式は次のようになります．

```math:
\frac{\partial\Sigma_{\rm g}}{\partial t}-\frac{1}{r}\frac{\partial}{\partial r}\left[3r^{1/2}\frac{\partial}{\partial r}\left(\Sigma_{\rm g}\nu r^{1/2}\right)\right]=0
```

$\Sigma_{\rm g}$はガスの面密度と呼ばれる量で，普通の質量密度を$\rho_{\rm g}$として次で定義されます（次元は${\rm ML^{-2}}$）．

```math:
\Sigma_{\rm g} = \int_{-\infty}^{\infty} \rho_{\rm g} dz
```

$\nu$は粘性係数です（次元は${\rm L^2T^{-1}}$）．
拡散方程式の導出は記事の後ろに書いてあるかもしれないです．

## 無次元化
では，この式を計算していきたいのですが，明らかに複雑でめんどくさそうですよね．

数値計算を行うときによく用いられる手法として，**無次元化**というものがあります．
様々な物理量（今回の場合，$r, \Sigma_{\rm g}, \nu, t$）をそれぞれ適当な値で規格化することで，式をスッキリさせて計算しやすくします．

円盤の典型的な半径，円盤の典型的な寿命，中心（$r=0$）での面密度をそれぞれ
$r_1, t_0, \Sigma_0$とします．
そして，次のように無次元の数$x, \sigma, \tau$を用いてそれぞれの物理量を変数変換します．

```math:
r \equiv x^2 r_1\\
\Sigma \equiv \sigma \Sigma_0\\
t \equiv \tau t_0
```
$\nu$についてはとりあえず置いておきましょう．
これらを元の方程式に代入すると，

```math:
\frac{\partial \sigma}{\partial \tau} - \frac{3t_0}{4r_1^2} \frac{1}{x^3} \frac{\partial^2}{\partial x^2}(\sigma \nu x) = 0.
```

ここで，$\nu \propto r$ と仮定し，
$\nu = 4r_1^2/(3t_0) x^2$ とおきます．

すると，次のようになります．

```math:
\frac{\partial \sigma}{\partial \tau} - \frac{1}{x^3} \frac{\partial^2}{\partial x^2}(\sigma x^3) = 0
```

$y = \sigma x^3$とおくと，

```math:
\frac{\partial y}{\partial \tau} - \frac{\partial^2 y}{\partial x^2} = 0
```

となって，よく見る拡散方程式の形になります．
これを計算して，最後に変数変換したものを元に戻せばOKです．

## 計算結果
Fortranで書いた計算コードを記します（学部時代に書いたのそのまま）．

<details><summary>サンプルコード</summary><div>

```Fortran:sample.f90
program main
  implicit none

  ! r方向の刻み
  INTEGER, parameter :: jmax = 200
  ! 時間ステップ数
  INTEGER, parameter :: nmax = 10000000
  ! 出力時間ステップ
  INTEGER, parameter :: nint = 10000

  REAL, dimension(0:jmax) :: u
  REAL, dimension(1:jmax-1) :: unew
  REAL delx, delt, r, x

  INTEGER j,n
  INTEGER i

  CHARACTER(len = 20) fname

  i = 0

  delx = 1.0 / real(jmax)
  delt = (delx**2) /3
  ! クーラン数
  r = delt / delx*2

  !初期条件
  u(0) = 0.0
  DO j = 1, jmax-1
     x = real(j)*delx
     u(j) = x*exp(-(x/0.1)**2)
  END DO
  u(jmax) = 0.0

  DO n = 0, nmax
     IF (mod(n, nint) == 0) THEN
        write(fname,'(i5.5,a)') , i, '.dat'
        open(10, file = fname)
        DO j = 1, jmax
           x = real(j)*delx
           ! 元々の変数に変換して出力
           write(10,'(2f12.5)') x**2, u(j)/(x**3)
        END DO
        CLOSE(10)
        i = i +1
     END IF

     ! 拡散方程式の差分方程式
     DO j = 1, jmax-1
        unew(j) = u(j+1)*r + u(j)*(1.0 - r*2.0) + u(j-1)*r
     END DO

     ! 値の更新
     DO j = 1, jmax-1
        u(j) = unew(j)
     END DO
  END DO

END program main
```
</div></details>

計算すると，こんな感じで面密度が拡散していきます
（横軸縦軸ともに，適当な値で規格化しています）．
![anime.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/403437/c9f49cf8-dad7-e676-b207-c9c5e1c0fc5b.gif)

（昔書いたコードを使ったので，動画はgnuplotで作りました）

## 自己相似解
さて，上の動画をキャプチャするとこんな感じになります．
色が濃くなるにつれて時間が進んでいくことを表しています．
![sample.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/403437/96d3e239-61e8-7f6a-e935-af478ab4a7f4.png)

よく見ると，各時刻の面密度のグラフは
他の時刻の面密度のグラフを伸び縮みさせた形になっています．
このような解を**自己相似解**といいます．

## 付録：拡散方程式の導出
力尽きたので，読者への宿題とする．

## 付録：拡散方程式の解析解
ここまで頑張って書いておいてなんですが，
実はこの方程式，特別な場合には解析解が存在します．

```math:
\nu = \nu_1(r/r_1)^{\gamma}\qquad （\nu_1 は r=r_1 での \nu）
```

と仮定すると，

```math:
\Sigma_{\rm g} = \frac{C}{3\pi \nu_1 (r/r_1)^{\gamma}} T^{-\frac{5/2-\gamma}{2-\gamma}} \exp\left[ -\frac{(r/r_1)^{2-\gamma}}{T} \right].
```
ここで，$C$は積分定数，

```math:
T = \frac{t}{t_{\rm dif}} +1,\\
t_{\rm dif} = \frac{1}{3(2-\gamma)^2}\frac{r_1^2}{\nu_1}
```

半径Rより内側の円盤の質量は，

```math:
M_{\rm d}(R, t) = \int_0^R 2\pi r \Sigma_{\rm g} {\rm d}r
```

と書けます．円盤全体の質量は，$R\rightarrow \infty$として求まり，
$\Sigma_{\rm g}$の式を使って頑張って計算すると，
$C$が求まります．

```math:
\Sigma_{\rm g} = \frac{M_{\rm d}(0)}{2\pi r_1^2 (r/r_1)^{\gamma}} T^{-\frac{5/2-\gamma}{2-\gamma}} \exp\left[ -\frac{(r/r_1)^{2-\gamma}}{T} \right]
```

## 参考文献
- [Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385](http://doi.org/10.1086/305277)