土屋 拓也

スライドショーには JavaScript が必要です。

講演者：土屋拓也（早稲田大学）
題　目：Einstein方程式の数値計算法　—その手法と数値安定性について—
日　時：2015年3月30日(月) 18:00-
場　所：早稲田大学西早稲田キャンパス63号館数学応用数理会議室（102号室）

講演

Ａ：それでは数理人セミナーの第5回ということで、本日は早稲田大学の理工学術院数学科の助教の土屋さんにお願いしました。土屋さんは早稲田大学で米田先生の所で学位を取得されました。では、本日は「Einstein方程式の数値計算法」ということでご講演をよろしくお願いいたします。

土屋：今、紹介された土屋と申します。僕はこの早稲田の数学科にいらっしゃる米田先生のもとで学位を取りました。米田先生がやっている研究は、数学科なんですけど一般相対論で、特にその分野ではEinstein方程式というのが支配方程式になっていて、その方程式についての数値解析をやっています。

　今日の内容ですが、7項目ありますが1個1個は重たいわけではないので最初にサラッと概要を言います。
　最初はまず数値相対論というのがどういうものなのかについて話します。2つ目にEinstein方程式の時空分解というもので、これがたぶん一番見慣れない内容だと思うんですけど、その動機から話していきたいと思います。3つ目の話が、拘束伝播方程式という方程式の話になります。これが、今回の中心的な役割を果たしてくれます。4つ目の話というのが、数値計算するときの不安定性という、要するに数値エラーみたいなものが増大していくというのを、どうにかして抑えなければいけないという話についてです。それをどうやって修正していけばいいのかというのが5番目と6番目の内容で、それぞれ補正システム、 $C^2$ -adjusted Systemという名前になっています。最後に、これらの内容を全て考えたうえで、数値計算で実際にやってみるとどうなるのかということを確認するというかたちになっています。

　それではまず最初に数値相対論の概要から話していきたいと思います。これはどういうものなのかというと、一般相対論における方程式に対して、数値計算によって解析を行う分野です。一般相対論に出てくる方程式の中で、測地線方程式とよばれているものとEinstein方程式が主要な方程式になっています。測地線方程式というのが何なのかというと、これは運動方程式です。Newtonの運動方程式の時空の曲がり具合の効果が入ったものになっています。次にEinstein方程式というのが場の方程式になっていて、電磁気を知っている方であれば、測地線がLorentz力の入った質点の運動方程式になっていて、Einstein方程式がMaxwell方程式という関係になっています。つまり、質点の運動を記述してくれる方程式と、場自体がどうなっているかを教えてくれる方程式の2種類あるという構造になっています。これらの式の数値計算全部を含めて数値相対論と呼んでもいいと僕は考えています。ただ、今回の話に登場する数値相対論という言葉自体は、Einstein方程式の数値計算で、その中でも初期値構成と発展計算の2つに関するものを指すことにします。

　まず、Einstein方程式自体があまりなじみがないと思うので、この説明からしたいと思います。Einstein方程式はどういうものなのかというと、時空の曲がり具合と物質の配置のつり合いの式になっています。時空の曲がり具合を解析・記述するために何が必要なのかというと、Riemannテンソル呼ばれる微分幾何学で用いる量で表現できることが知られていて、これを使ってやります.これはどんなのかというと、全部成分が0というのが時空が曲がっていない状態というのを表しています。これは計量 $g_{\mu\nu}$ というものによって全部成分が書き表されるので、このRiemannテンソルというのを全部計量を用いて表してやると、こういうかたちで書くことができます.
　ここの中に出てくる ${}^{(4)}\Gamma^{\lambda}{}_{\mu\nu}$ はChristoffel記号とよばれます。この記号の積やその微分した量で定義される ${}^{(4)}R_{\mu\nu}$ がRiemannテンソルとよばれる量です。Riemannテンソルは計量の2階微分というものになっています。
　あと補足事項としては、よく一般相対論の中で出てくる話なんですが、上の数字と下の数字が一致した場合には和を取る規則があります。これは「Einsteinの規約」とかよばれているものです。だいたい一般相対論の本には必ず書いてある内容なので、気になった方は見てみてください。

　Einstein方程式はどうなっているのかというのと、こんな式になっています。
　まず、 $G_{\mu\nu}$ という記号は $R_{\mu\nu}$ というRicciテンソルとよばれているものと $R$ というスカラー曲率とよばれる量であらわされます。RicciテンソルはさっきのRiemannテンソルの添え字を上下1つずつの和を取って減らしている。要するに添え字が4から2に減っている状態になっています。スカラー曲率という記号は、Ricciテンソルの添え字の和をとったものです。
　 $\Lambda$ というのが宇宙定数と呼ばれるもので、宇宙定数と計量で表されている項は「宇宙項」と呼ばれます。この宇宙項が入った方程式が一般的なEinstein方程式です。
　右辺の $8\pi$ 。これは円周率の8倍です。それから右辺の $T_{\mu\nu}$ は物質の配置を表すエネルギー・運動量テンソルとよばれる量です。
　この方程式自体を解くことでさまざまな宇宙の現象を知ることができます。ただ、この式はさっき言ったとおり2階の偏微分のかたちになっています。さらにこれは非線形性がかなり強いです。なので一般解は得ることはたぶん無理です。それで数値計算によるダイナミカルな計算というのが一般的に行われる手法となっています。
　このEinstein方程式を解く動機として、重力波の数値計算というのがあります。荷電粒子が動くと電磁波が放出されるのと同様に、質量があるものが動けば重力波は出てくるのではと考えられています。Einsten方程式に平坦な解の線形摂動として解の形を仮定してやると, 確かに波動方程式が導出されます。昔から重力波の研究は行われているんですが、まだ直接観測はされていません。間接的な発見は、TaylorとHulseという2人が、PSR B1913+16というパルサーの軌道の観測結果に対して、重力波によるエネルギー減少を計算したところ、観測結果と計算結果が一致したということで、重力波の存在の間接的証明となっています。ただ、直接的に観測はまだできていないので、その観測をしようということでKAGRAという名称で大型低温重力波望遠鏡を今、日本の神岡鉱山の地下に建設を行っています。
Ａ：どこどこ？
土屋：岐阜県の神岡鉱山です。
Ａ：神岡鉱山？
土屋：あそこには小柴昌俊先生がNobel賞を受賞した研究のときに用いたカミオカンデというニュートリノ観測装置があるんですけど、地下に新たにまた穴を掘りました。Michelson-Morleyの実験というのがあって、それは光を何回も反射もさせて、その反射した距離を測ってズレが出るか出ないかというのでこの観測というのを行うんですけど、距離が長くなきゃいけないというのがあるので、そこで地下に3kmの穴を掘りました。感度を良くするためには震動が来ない所で、かつ熱が入ると熱による雑音が入るので、どうしても低温のところが欲しいということで地下でやろうということになったのだと思います。

　数値計算をしたいという話に戻るんですが、結構困った状況にあって、Einstein方程式は4次元共変といわれている形式になっています。4次元共変はどういうものなのかというと、時間と空間が同じものとして扱われています。数値計算するということは、時間軸が1個あって空間がそれ以外あって、そこを発展させていくというのが普通です。そのため、時間1次元と空間3次元に分かれてないといけないんですけど、もともとごちゃ混ぜになっているのをどうやって分ければいいのかというのが昔から問題になっていました。その問題をまず解決するために、時空分解ということをやらなきゃいけない。要するに、Einstein方程式そのままだと数値計算ができないという構造になっています。Einstein方程式は時間軸というのが固定されていない。もうちょっと言うと、時間軸だと思っていたものを取って計算していくと、いつのまにか空間軸に変わっちゃったりするということが起こり得ます。だから計算がそこで破綻してしまうという構造になっているので、何とかして時間軸を取り出すいうのを最初にやる必要があります。
　話が少し変わりますが、20世紀の前半ぐらいに量子力学が構築されました。その時に一般相対論だけは含まれていなかったので、一般相対論も量子化したかった。量子力学のほうは、時間と空間というのがきれいに分かれた状態で構築されていたので、合わせるためには一般相対論のほうも何とかして時間と空間を分けなきゃいけないというのがありました。つまり、もともとは量子重力理論というのをやるために研究が40年ぐらいされていて、1960年ぐらいになってようやくArnowitt、Deser、Misnerという3人によって、うまく時間軸を取る方法というのが考案されたというふうになっています。この3人の頭文字をとって「ADM formulation」とよんでいるのが、今回出てくるメインの式になっています。数値相対論で用いる時空分解の方法はその時の遺産を使っています。
　このADM formulationという方程式系は、正準形式になっているのでそれはそれでよい性質を持っていますが、ただ数値計算がうまくいかなかったときに、いったい何の物理量がどううまくいかないのかということが良く分かりません。感覚的に捉えやすい幾何な表現にできないかということで、その後10年ぐらいした時にYorkとかSmarrという人たちが、もともとの4次元のRiemann多様体を時間一定の3次元Riemann多様体という葉相構造と捉え直して、新しい方程式系を作りました。YorkとSmarrによって作られた形式は正準形式ではないですが, 方程式が3次元のテンソル構造が保たれている状態できれいになっています。今では、こっちのほうが普通ADM formulationとよばれています。

　ここから先は、実際に幾何的にどうやって4次元の多様体中のEinstein方程式を3次元に分解していくのかという話をします。
　まず時空の分解という話ですが、3次元の超曲面が重なり合っている状態というのが4次元多様体だという概念の下で、4次元のものを3次元に分解するのがこの時空分解という方法になっています。　まず、超曲面上任意の場所に法線が取れます。その法線に対してちょっと時間がたったところの面というのを考えたときに、その動きというのを追ってあげます。どうなっているのかというと、この場合法線の方向と時間軸が異なっているのが一般です。そのずれというのをシフトベクトルというもので表していて、要するにこのシフトベクトルが0だったら法線の方向と時間軸は一致してくれます。だからシフトベクトルが存在していると時間軸と法線の方向が一致しない。
　また、 $\alpha$ というのが時間軸の幅というのを表しラプス関数と呼ばれています。さっきでてきたシフトベクトル $\beta^i$ というのが、時間軸と法線方向のずれを表し, シフトベクトルと呼ばれています。さらにこの $\gamma_{ij}$ というのが3次元曲面上の計量となっています。その反変量 $\gamma^{ij}$ が $\gamma_{ij}$ の逆行列という関係になっています。こういうふうに4次元計量 $g_{\mu\nu}$ を分解するのが時空分解になります。
　さらに、3次元計量 $\gamma_{ij}$ を使ってあげて外部曲率というのを定義します。 $\gamma_{ij}$ が法線方向に沿ってどう動くかというのを表しているのが、この外部曲率と思ってください。ここに出てくる $\mathcal{L}$ というのがLie微分を表しています。Lie微分というのはどういうものなのかというと、ここに出てくる添字のベクトルの方向に沿って動くと変化しないという量になっています。力学的な見方をすれば、 $\gamma_{ij}$ というものに対する速度みたいなものが $K_{ij}$ だと思ってください。

　ここから先はEinstein方程式の時空分解をやっていきます。どうやって分解を行っていくのかということなんですが、さっき言った法線ベクトル $n_{\mu}$ というものと、それに対して直交する射影テンソルというのを使います。要するに法線の方向と法線じゃない方向にきれいにスプリットするというのが、今からやっていく話です。
　そうすると2階のテンソルは、時間・時間方向というものと、時間・空間方向と、空間・空間方向というように分解ができます。

　Einstein方程式というものは2階のテンソルなので、今言った3つの成分に分けることができます。最初にまず、Riemannテンソルの分解から説明します。4次元Riemannテンソルというのを3次元部分に分解します。Riemannテンソルをそれぞれ時間・時間方向と空間・時間方向へ2つの成分に分けると、これは微分幾何の中で出てくるものなんですが、Gauss-Codazzi方程式という式とCodazzi-Mainardi方程式というのがでてきます。さらにRiemannテンソルは空間・空間方向にも分けることができます。結果、Riemannテンソルを3次元の量で表してあげると、この3種類が出てくるというかたちになっています。

　今度は今やったRiemannテンソルの分解式からを、Ricciテンソルの分解式を作ります。　Riemannテンソルから作った3種の式から添字を2つ全部消していってあげると3つの式が出てきます。これらの式が、後々Einstein方程式の時空分解で出てきます。　
　あとは、残っているスカラー曲率も3次元の量で書き表します。そうすると4次元のスカラー曲率が3次元スカラー曲率や外部曲率などで表せます。これらを使ってやってEinstein方程式をスプリットします。

　まず、Einstein方程式に対して時間・時間方向でスプリットする。時間・時間成分だけ取ってきてあげるとHamilton拘束方程式とよばれているものがでてきます。これは、エネルギーみたいなものを表していると思ってください。時間に対応したものがエネルギー、空間に対応したものが運動量という構造になっているので、これは時間成分を取り出しているので、エネルギーみたいなものとして扱われています。
　今度は空間部分が入ったところを取ってきてあげる。空間・時間という部分になっているんですが、これも恒等的に0になる量として知られていて、運動量拘束方程式とよばれている。
　最後に、空間・空間方向へ分解をやってあげると、計量の時間発展の式になっていて、結局、2種類の拘束方程式と1種類の時間発展方程式というふうに分けることができました。

　まとめると、Hamilton拘束方程式1本と、運動量拘束方程式の添え字が1から3まであるので、3本の式があります。空間・空間から得られる発展方程式は時間2階になっているので、時間1階の連立の形に直してあげます。それには、 $\gamma_{ij}$ の発展方程式というのが、外部曲率というところの定義式があったんですが、それを使います。この4種のセットで、拘束方程式が1本と3本と、 $\gamma_{ij}$ の発展方程式が6本。 $K_{ij}$ の発展が同様に6本あり計12本なので、4本の拘束条件と12本の発展方程式というかたちになっています。
　実際に数値計算する時は、当然、拘束方程式を解くわけじゃなくて、この $\gamma_{ij}$ の発展方程式と $K_{ij}$ の発展方程式を解くことによって次の $\gamma_{ij}$ 、次の $K_{ij}$ を求めて、その2つの量からその後の時間発展した後の時空構造がどうなっているのかというのを見ます。拘束方程式はその時に計算が満たされているかどうかをチェックするための量として一般的に使われています。

　Einstein方程式の時空分解なんですが、さっき説明した時空分解の方法で発展方程式系が一応得られるんですが、一意的に決まらない。実際にこれは数学的に言うと、この方程式系を同値変形しても新しい方程式系ができてしまう。例えば発展方程式にHamilton拘束条件をくっつけるとか運動量拘束方程式もくっつけるとかいうのをやって、新しく式を作っていっても別に全然かまわないので、そういう意味で同値変形することが可能です。だから数学的に全く変わらないものをいくらでも作れるんですが、ただ、数値計算させるときに発展方程式のかたちがどうなっているかというのが結局、数値安定性に絡んできちゃうので、どういうふうな発展方程式を作っていくと数値安定なものができるかというのが問題としてあります。　こんなことを言いだすのは、さっき作ったADM formulationという系自体が、実は数値的に全然安定じゃなくてうまくいかないというのがよく知られているからです。計算させていくと理想はエラーが単調増加しないものです。なのに、数値計算をやってあげるとどうにもこうにも時間とともに発散しちゃって、全然欲しいところの時間まで計算が続かないというのが昔からよく知られている話で、これは「formulation問題」とよばれています。
　数値スキームをRunge-Kuttaなり何なりというのをいろいろ変えても、発散してしまいどうも数値スキームと発展方程式の相性が悪い。数値スキームをどうやって開発するかというのは結構ヘビーな話なので、それをやるのも一つなんですけど、まず方程式のかたちも悪いだろうということで、それをどうにかしましょうというので1990年代ぐらいにその研究がずいぶんされてきました。今言ったADM formulationというもの以外の、数学的には同値変形なものをいろいろ作って、何か新しいものを作らなきゃいけないというふうになりました。それをうまく作ったのが、京都大学の中村先生と柴田先生という2人で、95年ぐらいに作りました。それをBaumgarteとShapiroという二人がすこし改良し、今では「BSSN formulation」とよばれているものが広く使われています。

　ただ、個人的には数値安定的な方程式というのをどうやって構築すべきかという指針をつけるべきだと考えています。それを見るために、拘束伝播方程式というのがキーになってきます。この拘束伝播方程式というのが何なのかということですけど、拘束方程式はHamilton拘束方程式と運動量拘束方程式の2種類がありました。この拘束方程式の時間発展を求めてやると、その後どういうふうにして拘束条件の破れが発展していくのかというのを見ることができます。これを得るには、さっき言ったHamilton拘束方程式と運動量拘束方程式というのを普通に時間微分して計算してあげればよくて、2種類の発展方程式が出てきます。これら式でまず分かることというのは、時間微分したときにHamilton拘束条件と運動量拘束条件の線形和の形で書き表すことができます。なので、初期にHamilton拘束方程式と運動量拘束方程式というのが0というふうになっているんだったら、発展させたもの自体が0になっているので、その後もずっと0であるということが、解析的に保証はされています。ただし、解析的に保証されているからといっても数値的には必ずエラーが出るわけで、そのエラーというのがどうやって蓄積されていくのかというのをこれらの式を解析することで知ることができます。

　ここでまず最初にする話というのが全然数値相対論に限らない話で、ただの常微分的な話になります。概念から話すと、こういう式で $a$ は定数だと思ってください。こういう式があったときに、解としては定数 $C_0$ として $C_0e^{at}$ となります。この解の挙動というのは、 $a$ が正になっていれば当然発散、 $a$ のところが負になっていれば $t$ が無限大になってくれれば0になっていくので、 $C_0$ が何であろうともuは0になっていくという構造になっています。
　今度、これを行列のパターンにした場合はどうなるかというと、 $A$ というものが係数行列で与えられたとして、その行列の固有値を $\lambda_i$ だとしたときに、一応こういうふうに式を書くことができます。これの解を出してあげると、こんなかたちになります。解に固有値の部分が肩に乗ってきてというふうになって、同様にこれも $t \rightarrow \infty$ にしてあげると、 $\lambda_i$ のところが正だったら発散、 $\lambda_i$ が負だったら解が0に収束していくというかたちになっています。

　この話を、さっき言った拘束伝播方程式のときの解析に使えないかというふうに米田先生らが思ったようです。　拘束方程式の時間発展の係数行列がもしその固有値がすべて負になっていてくれるんだったら、時間が進んだとしても0に収束するということが期待できます。要するに、何か少し発散していく様子が見られたとしても、そのうちまた0に収束していってくれるということで、漸近的に安定していってくれるんじゃないかということになっています。たぶんこれは数値相対論に限らず、他の所もよくやられている話だとは思うので、それを数値相対論でやってみたという話になっています。
　この話は結局拘束方程式の係数行列の固有値の符号と大きさを見ているので、Constraint Amplification Factors（CAF）という名前で呼んでいます。このCAF解析というのを米田先生と共同研究の真貝先生という人が提案しました。要するに、ADM formulationというのは悪いので、何とかして同値変形して方程式が作り変える。そうすると、さっきの拘束伝播方程式というのも影響を受けるので、その影響を受けた結果が漸近的に0にいくような構造をもつように作ってあげれば、うまくいっている式になるんじゃないのかというふうな指針です。

　やり方としては、拘束伝播方程式というのの中には空間微分が入っているのでこのままだと固有値が求められない。代数的にやるためにFourier変換してあげて、微分を代数的な値に変更してあげることで解析ができます。
　見方としては、固有値の実部が負であれば安定、正だったら不安定という構造になっているというふうに見ます。要するに、固有値の部分が正であるということはさっきのところで言った固有値が正になっている状態なので、それがそのまま時間発展していくとどんどん時間とともに増えていくのでダメ。負になっていれば、時間発展させていってもだんだんそれが抑えられていくので安定するという感じです。では虚数部分はどうなっているかということなんですけど、虚数部分自体はあれは結局回転を表しているだけなので、だからずれるというだけであって、安定性に直接的な寄与はしないというふうになっています。でも一応、無いよりはあったほうがよいというふうに考えています。虚部というのは、あれば誤差を少し散らしてくれるんじゃないかというぐらいなイメージになっています。
　係数行列には $\alpha$ とか $\gamma_{ij}$ とかさっきやった変数が含まれているので、普通は背景時空というのを固定する必要があります。背景時空というのは、数値計算する時に最初に近い時空を仮定してやる。例えばブラックホールの時空を計算したいとなったら、ブラックホールの時空を最初に解として与えてあげて、その解がその後で時間発展のどこまで行くのかを計算します。これはどういう時空に対する解析をやっているのかというのをまず考えた上で、背景時空を設定する必要があります。これは一般性がなくなるので本当はやりたくないんですけど、やらないと解析ができないので。
　実際に行列のかたちでさっき言った拘束伝播方程式というのを作ってあげると、こんなかたちで書くことができます。これに対してFourier変換をやってあげます。一番簡単な場合でフラットな時空を考えます。フラットな時空というのは要するに何もない状態というふうになっているんですけど、その場合だと、 $\alpha$ が1で、 $\beta^i, K_{ij}$ が0、 $\gamma_{ij}$ が $\delta_{ij}$ となっています。要するに、4次元計量 $g_{\mu\nu}$ が対角行列でその成分が $-1,1,1,1$ になっているという状態になっています。これを入れてあげると、ここの部分の係数行列はこんな感じになっています。この固有値を出してあげると $(0,0,\pm{\it i}|k|)$ と虚部のみが出てきている。虚部しかないという状態なので、「不安定ではない」ということは分かるんですけど、「安定でもあるとはいえない」ということになっています。

　今度は、よく数値計算でやられているペナルティ法という方法があるんですけど、それをやったときにどうなるかというのを考えてあげます。要は発展方程式に拘束条件をを付けてあげて、それでブレーキをかけてあげるというイメージになっています。さっきのADMという式をいろいろといじるときに拘束値をくっつけてあげて、式を変えてあげるという感じになっています。問題なのは、どうやってその項を付ければいいのか、付けた後、どうやってそれがいいか悪いかを判断すればいいのかということです。これがここからの話になります。方法としては、発展方程式に拘束値の微分が0階や1階、2階で作った項をくっつけてあげるとこのようになって、拘束伝播方程式が影響を受けて、元のものに対してここの付けた部分からそれの影響値というのが出てきます。さっき登場したこの部分のCAFというのを求めてあげたときに、この部分からの影響が解を安定にする方向にさえ働いてくれれば、その系はうまくいっていると考えられます。

　ADM formualtionに対して、例えばこんなふうにしていじくってあげます。拘束値自体は別に変えるものじゃないので変えなくて、発展方程式のところにこんな拘束値から作った項を付けてあげます。

　そうすると、拘束伝播方程式をFourier変換して背景にフラットなものを入れてあげると、こういうふうに変わってきてくれます。これの固有値を計算するとどうなっているかというとこんな具合になります。ここで $\kappa_D$ と書いたのが、パラメータだと思ってください。これを全部0にしてあげれば元のADM formualtionでの場合に戻って、 $\kappa_D$ に何らかの値が入っていると、それとは違った状態になっていると思ってください。それをやってあげると、このかたちを見る限り、 $\kappa_D$ が正であれば固有値の実部が負ということになるので、ADM formulationよりも安定であるというのが解析的にいえます。こうやってやれば実際にできるんですけど、どうやってさっき言った付加項を作ればいいのかというのがよく分かりません。どうやって作ったのかは論文にも書いてないです。こうやってやるとうまくいくしか書いてないので。
Ａ：「うまくいく」というのは、固有値が全部0にならないという意味での「いっている」という意味？
土屋：固有値の実部が負になっているという意味でうまくいっている。実部が正のものよりは0のもののほうがいい。実部だけだと負のものがいい。虚数部分は無いよりあったほうがいいという感じなので。ADM formulationの場合は固有値に0があって、残りの部分に純虚数が入ってきているだけでした。今の場合だと実部に負の部分の成分が追加されている。だからこっちのほうがいいという感じの式になっています。

　さっきも言いましたが、付加項の作り方がよくわからないというのがまず問題にあって、さらに背景時空というのが要するに元の初期値みたいなものを与えているんですけど、それは実際に数値計算によって時間発展していくと必ずその背景時空からずれてきます。実際の数値計算ではずれていくので、その後のずれた場合でもこの解析の結果は大丈夫なのかという問題があります。
　なので、背景時空とかに依存しないように付加項を設定する必要というのがある。要するに、この解析自体が背景を固定しているというのもあるので、背景を固定した上でこの話はできているので、実際に背景を固定しないでうまくいく方法とかはないかというので考案されたのが、 $C^2$ -adjusted Systemとよばれているものです。

　どういうふうにやっているのかというのをちょっと説明していくと、さっき言ったとおり、拘束条件付きの発展方程式はこういうふうに書けます。発展方程式があって、拘束条件の式があって、この拘束条件の式を使ってやって発展方程式をこういうふうに補正してあげます。発展方程式に対して何か拘束条件の2乗したものを発展変数で変分しているという構造になっています。こういうふうにしてやってやると何がいいのかというと、2乗したものの時間微分というのを考えてあげるときに、元のものに対してこの部分の影響というのは拘束伝播方程式にこんなかたちで出てきます。なので、この部分 $\kappa^{ij}$ がもし正定値になっているのであれば、付加の寄与は全体が負。要するに、元のもの2乗の値に対して必ず負のものというのが出てくるので、それが負の方向に抑えられていくという構造になっている。あとはこの $\kappa^{ij}$ のところを正定値でうまく取ってあげさえすれば、この部分がうまく負の方向に走っていってくれる。時間発展したら乱れたものを戻してくれるというのができているという構造になっています。
Ａ：これは目的ですか？　これは新しく項を変えてますよね？　 $C^2$ -adjustedで。これは目的なんですか？　それともその結果なんですか？「以下のようになる」と書いてあるけど、そうしたかったということなんじゃないの？　そこがよく分からない。
土屋：たぶん最初に考えたのはおそらく拘束伝播方程式の形で、どう発展方程式に項を付けるのかが問題だったと思います。今、山中さんがおっしゃられたように、おそらくこれが最初にあったような気はします。いきなりこれをやったら偶然うまくいったとはちょっと思えないので。もしかしたら似たような式を見た時に「こうやってやればうまくいくのかな」と思って、こっちを付けたのかなという気はします。

　今言った話というのを、さっき言ったADM formulationに適用してあげるというのが僕がやったやつで、実際にこのCの2乗というのはこんなかたちで表されています。拘束値の2乗の和みたいなもので作られています。これらに対して発展変数というのでそれぞれ変分してあげたものをくっつけてあげて、方程式を変えてあげます。そうすると、当然さっきの話からすると付加の部分がうまく負の方向に行ってくれるはずだというふうになっています。

　拘束伝播方程式の変化を見てみます。背景をフラットなものにして、このときに $\kappa_\gamma{}_{ijmn}$ と $\kappa_K{}_{ijmn}$ を正定値にすると、拘束伝播方程式は、元の式に対してこんな量がくっついていてくれます。この式を見てあげると、上式右辺第2項目と下式右辺第2項目にdamping term、要するに自分自身の所のマイナスみたいなものが出てきてくれているので、自分自身の破れを抑えるという構造になっています。この辺が安定化のキーになっていて全体的に拘束値の破れが抑えられていくのかなと思います。

　これをCAF解析という係数行列の固有値の話に持っていってあげます。このときFourier変換した式としてはこんなかたちで出てきます。これの部分の固有値を求めてあげると、こんなかたちで出てきます。この部分はどうなっているのかと見てあげると、 $\kappa$ というものが正に与えられているんだったら、固有値として出てくるものの実部は負になっているという状況になっています。なので、安定だということがいえます。

　実際に今言ったものを実際に数値計算するとどうなるかというのをこれから先ちょっと見ていきます。数値計算する例としてよく使われているGowdy waveという式があるんですけど、これはEinstein方程式の厳密解の一つになっています。どういうものなのかというのをちょっと説明すると、これ自体は初期に特異点があって、そこから時間とともに宇宙が膨張していくという様子を表している例になっています。要するに、今の宇宙論というのは初期に何かあって、その後インフレーションが起きて膨張していって、その膨張していったものが現在の宇宙になっているというふうになっているので、それを一応模したいということでたぶんGowdyという人が作ったんであろうと思うんですけど、宇宙の膨張を表してくれている式になっています。実際にこれを数値計算するときには、膨張していく方向に見ているんではなく、時間反転させて初期の特異点の方向に向かうように計算します。これは、数値相対論でやりたいものはブラックホールの計算であって、ブラックホールの計算というのは必ずsingularityを含むということで、そのsingularityが発生した時にうまく計算がいってほしいということの確認がしたいためです。逆方向に時間発展していくとsingularityに近づいていくので、だんだんと数値エラーというのが積算しやすくなっていてエラーが出やすいという構造になっています。
　あと数値計算するために決めなきゃいけないのは座標条件とよばれているもので、これはゲージ条件なんですけど、よく取られるのは調和条件といわれているものなんですが、こうやって取ると時間発展がうまくいくというのが提案されていて、これを使っています。この説明をし始めると結構いろいろとあるので、これはこういうふうに $\alpha$ というのを追っているというふうに思ってください。
　境界条件もいろいろと変わると大変なので、周期境界条件で今回はやっています。離散化の手法としてはCrank-Nicolsonを使ってやっています。

　数値計算の結果を見てあげると、こんな図が出来上がります。この横軸というのは時間発展の方向。ただし、さっき言ったとおり時間を逆方向にしているので、時間としてはマイナスの時間だと思ってください。縦のほうがCの2乗のlog10になっているんですけど、要するにこれはエラーです。エラーがどうなっているのかというのが縦軸で、横軸が時間の方向になっています。
　この式を見てあげると、Standard ADMと書いてあるのが最初に出したADM formulationで、この赤いライン。Detweilerと書いてあるのが、補正の例というところで紹介したもので、この緑色のラインで、こんなふうな感じになっています。最後にやった $C^2$ -adjustedというやつがこの青いラインになっています。
　この図を見てあげる限りでは、赤も緑も時間発展とともにどんどんどんどん伸びていって、これはlogでとっているのでエラーが10の何乗みたいな感じでどんどんどんどん飛んでいくので、もう計算としてはおかしな状態になっている。この $C^2$ -adjusted ADM formulationでやったものだけがエラーが下がっている状態になっているので、これからするとまだこっちのほうがいいという状態になっています。つまり、この $C^2$ -adjusted ADM formulationのほうが拘束値の破れが減少したというふうになっています。
　一応こういうふうな解析をすることによって、どうやってADM formulation、つまり、元の方程式を「良いもの」に変えていけるかというのをお話ししました。一応、スライドとしてはこれで終わりです。

質疑

Ａ：なんで途中がこうなるの？　時間がどっちに見たらいいのか分からないけど、左から右に見るか、右から左に見るのか分からない。
土屋：左から右です。
Ａ：なんでこう……
土屋：こういうふうになっている理由ですか？
Ａ：なんでこうなるの？
土屋：それは、 $\alpha$ って何だったかというのを思い出してみると、 $\alpha$ というのは最初の面と次の面の時間の幅を表しています。実はこの調和条件というゲージ条件で $\alpha$ の時間発展をやってあげると、ここに図は無いですが、 $\alpha$ の時間発展の様子を表示をしてあげると、ある時刻に増えて、次の時刻にあまり進まなくて、その次に増えてというふうになっているので、α自体の進み方が一様じゃないんです。その一様じゃないものを使って時間を作っているので、こんな感じで増えたり減ったりみたいな感じになると思います。
Ａ：提案された $C^2$ というのが分かんないけど……
土屋：青いやつ。
Ａ：青いやつがまた誤差が途中から小さくなるというのは、安定性だけが理由なの？　これを見ると、「安定性、良くなってるじゃん」というふうに見ていいんですか？
土屋：見ていいんです。というより、ここで出している量自体は数値エラーなんですけど、エラーの全てを表しているわけじゃないと思ってます。それは、この縦軸は拘束条件を満たしているかどうかというものだけ見ています。拘束条件が0からどのぐらい離れているかというのだけ見てる。なので、拘束条件が破れてたけど低くなっているというのは、上がるものよりはいいという。数値計算のエラーがどこかで減っていくというのはかなり不自然なので、おそらくですけど、エラーが減っているというのはたぶん他の所にエラーを散らしているだけで、どこかの飛び出たエラーのせいで計算が止まるので、それをうまく散らすことによって何とか計算を生き延びさせようとやっているというのが、この手法かなという気がします。だから実際のエラー積算自体は減っているはずはなく、必ず増えているはずです。要するに数値相対論ではどういう感じになっているかというと、拘束条件は最低限満たすようにする。それを計測する以外に他にエラーを計るものが無いので。そこを追った上で最終的に数値計算のエラーって全部どうなっているのかというのを見るのはできると思うんですけど、そこまで研究はまだ行っていない。そういう分野ができるところではまだないので。そこまでできれば素晴らしいんですけど、なかなかそれができない。
Ａ：確かに誤差は減るか、誤差みたいなものは持っていて突然また落ち始めたら、確かに悪くなるよりかはいいですよね。
土屋：要するに、ここを超えるということは、拘束条件の式が満たされないということを表しています。これが満たされないということは、この4種の式が全て満たされた状態がEinstein方程式が満たされている状態なので、満たされていないということはEinstein方程式じゃないところの計算をやっちゃっている。そうするとまずそれがまずい。その部分ができた上で、もうちょっと詳しい話ができるけど、拘束条件が満たされていないところでやるものじゃないという感じですね。
Ａ：少なくとも青いやつは満たされているっていう？
土屋：そう。まだ満たされている。まだいいというような。本当はもうちょっと解析的な話をきちんとやるべきではあるんですけど、いかんせん方程式があまりにも複雑過ぎるので大変ですね。
Ａ：先ほどの図を作るのに、数値計算結果を作るのに、どのぐらいの実行時間がかかっているんですか？
土屋：赤い線が10分ぐらい、緑が15分ぐらい、青が半日ぐらいです。
Ａ：この後どうなっていくかというのは、だいたいご存じなんですか？
土屋：まず、緑はちょっと忘れましたけど、赤はずっと行きます。そのままズーッと上まで行って、どこまで行ったかな。表記している所は $t=1,600$ じゃないですか。4万ぐらいにまでは行きましたね、とりあえずエラーが増大はしつつそのままの状態でずっと計算が続きます。
　青はたぶんどこかで破綻していたような気がします。たぶんエラーを散らしている状態の方法が悪いか何かで正常じゃない可能性があると思っています。だから、これもやっぱり結局、完全な方法じゃないんです。もっとうまい方法があるとは思うんですけど。Crank-Nicolsonなんか使わないでRunge-Kuttaを使ったりすると、たぶんもうちょっと良くなったりする気もするんですが、ただ、たぶん構造的にCrank-NicolsonとかRunge-Kuttaでやっても、ベストになっているかちょっと分からない。たぶん何かしらの意味でエラーは蓄積されます。
Ｂ：たぶん数値計算をしようとしたら、ADM formulationとかの固有値が全部0になっちゃうというのは、それは解けないという印象が僕はあるんですけど、そこでペナルティをこう入れるのはすごく一般的というかデータがあるので、それでいいスキームができたんだというのはすごい面白いなと思いました。ここでやっぱり元の方程式系が崩れていると見てもいいかどうか……。
土屋：そこなんですよね。僕もレフリーにその辺を文句つけられました。他の補正の仕方がたぶんあるとは思うんですけど、僕がやった補正だと拘束条件の2乗というのがあるわけで、このHamiltonian拘束条件に $\gamma_{ij}$ が2階微分が入っています。それの2乗を変分していることになるので、補正項の部分に $\gamma_{ij}$ の4階微分が出てきます。運動量拘束条件には $K_{ij}$ の1階微分が入っています。それを2乗しているので、変分してやると $K_{ij}$ の2階微分が出てきます。もともと $\gamma_{ij}$ の発展方程式には $\alpha$ の項に微分なんか入ってなかったにもかかわらず、補正項に $\gamma_{ij}$ の高階が出てきちゃったりとか、 $K_{ij}$ の発展方程式に $K_{ij}$ の2階が入ります。
　Einstein方程式は双曲型で、そのため重力波というものが出てくる。それからするとこんなものをくっつけてあげるとむしろ放物っぽい形になっちゃってるので、「大丈夫なの？」というのは言われたんですけど。そうですね、その辺が問題といえば問題ともいえます。
　ただ、これは厳密に双曲型じゃないので、双曲と放物の混合型になっているので、要するにきれいな波動にはなっていない。そういう意味で発展の部分を波動型に変えてあげるというのをやっているグループがいろいろと別の形式というのを提案しています。ただ、そうすると今度は双曲型になるよういろいろと発展方程式の構造をいじくっています。それも言い換えてあげれば、波動の部分だけ取り出してきて他の部分を全部捨てるようなものだと思っています。重力波ということで波動が伝播するように作った方が良いというけど、実はそれが本当に必要なところだけ取り出してきて要らないところを切っているのか、要るところまで切ってしまっていないのかがちょっと分からない。あまりにもいろいろなやり方があり過ぎてよく分からない。そこで、CAFのような解析でやっていくんであれば、発展方程式の形と拘束伝播の形からどのようなものであるのがいいのかが統一的にわかると考えています。
Ｂ：ペナルティ項の選び方はこちらの自由ですよね？
土屋：自由です。ADM formulationと同値であればいいので、何でもありです。あと問題になってくるのは、この辺のパラメータをどういうふうなサイズにすべきかとか、どうとるべきかというのが残ってはいますが、その辺も経験則でしかわからないので、ちょっとやり様が無いぐらいな感じにはなっています。
Ａ：ありがとうございました。他に質問とかありますか。特には無いですか。分かりましたか。それではどうもありがとうございました。（拍手）
土屋：ありがとうございました。

プロフィール

2012年3月早稲田大学大学院基幹理工学研究科数学応用数理専攻中途退学。博士（理学）。早稲田大学理工学術院助手、早稲田大学理工学術院助教を経て、現在千葉工業大学工学部未来ロボティクス学科、東京女子医科大学医学部、日本女子大学理学部数物科学科、早稲田大学基幹理工学部非常勤講師。Einstein方程式に関する研究を主要研究分野として、特に初期値問題を数値的手法や解析的手法を用いて行っている。日本数学会、日本応用数理学会会員。