ms.dvi

絶対温度が負となる点渦系に関する力学的考察

∼専用計算機を用いたダイレクトシミュレーション結果

および解析的結果∼

靜岡大学 八柳祐一

神奈川大学 羽鳥尹承

研究協力者 京都大学

(名誉教授)

 際本泰士

2010

年

11

月

21

日

1

はじめに

温度をパラメタとして，点渦系のダイナミクスを理解するための研究である。 本原稿は，研究室の学生向けの review にも活用しようと思い，これだけである程 度完結した形にしたかったため，本来は「参考文献」として参照を示すに留めるべ きと思われる所も，冗長ではあるが，内容を再掲載した部分がある。また，背景な どについても，比較的丁寧に説明したつもりである。よって，他の文献と一部内容 が重なる部分もあるが，ご容赦願いたい。 負温度点渦系の研究は，Onsager [1] に始まると思って良い。それ以降，様々な理 論的，数値的な研究がされてきたが，特に数値的な結果において，点渦系の数値計 算はそれなりに計算パワーを必要とすること，および負温度性を議論するためには 系が境界に閉じ込められていることが必須であるが，境界がある系の扱いは数値的 にも煩雑であること等が理由となり，研究を始めた当時調べた範囲で，大規模なシ ミュレーション結果は，周期的境界条件下の矩形領域にとどまっていた。また，そも そも我々は，円筒容器に閉じ込められた純電子プラズマの解析に対応した 2 次元点 渦シミュレーションを行いたかったため，境界条件は円形にする必要があった。以 上の理由から，我々は，円筒境界に閉じ込められた点渦系のシミュレーション研究 を始めることを決意した。 本稿では，まず非中性純電子プラズマを用いた渦実験の成果を紹介する。次に，こ の実験から派生して行われた点渦系シミュレーションについて，なぜ，2 次元点渦 系では絶対温度が負となり得るのか，という点を中心に紹介する。そして，我々が 行ってきた数値的，解析的結果の紹介を行う。

図 1: 実験装置の図。写真中央にある円筒形の容器が超電導コイルで囲まれた真空 容器で，この中に電子を閉じ込める。左側の四角い箱の中には光学系が設置されて おり，電子の像を箱左端の CCD カメラへ導く役割を担う。このようにして超電導マ グネットから CCD カメラを遠ざけないと，強磁場の影響を受け正しい撮像が行え ない。

2

本研究の背景

本研究のスターティングポイントは，京都大学の際本研究室で行われていた非中 性プラズマを用いた渦実験である [2–4]。図 1 に実験装置の写真を示す。普通，プラ ズマ状態といえば，元々電気的に中性だった原子／分子が電離してイオンと電子に 分離した状態なので，巨視的に見れば電気的には中性となる。非中性プラズマとは， この電気的な中性条件が破れた状態のプラズマを指し，電気的に負の粒子である電 子のみから構成されたプラズマ (純電子プラズマ) などが該当する。 この実験装置では，純電子プラズマを円筒容器の軸方向に超電導コイルで印加し た最大 2T の磁場により径方向に，両端の電極に印加した負電位により軸方向に，そ れぞれ閉じ込める (図 2)。すると，導体の円筒容器と電子の間に径方向内向きの電 場が発生する。この電場と軸方向磁場により，電子はE × B 方向に回転運動をする (理由は後述)。一方，電子は磁場に平行方向には自由に運動することができるので， 両端の負電位の間を高速に往復運動する。すなわち，電子は全体として軸方向に高 速な往復運動を繰り返しながら，軸の周りにゆっくり回転していく。 2T における典型的パラメタを表 1 に示す。ここで注目すべき点は，磁場に垂直面 内での運動の時間スケールが 100ms のオーダなのに対して，磁場方向の往復運動の 周期が 3ms と約 2 桁短いことである。これは，すなわち，磁場に垂直な方向のどこ

B

E

Rotation

Trap and hold

図 2: 電子が閉じ込められる様子。電子は，円筒容器の軸方向に超電導コイルで印 加した最大 2T の磁場により径方向に，両端の電極に印加した負電位により軸方向 に，それぞれ閉じ込められ，導体でできた円筒境界との間に内向き電場を発生させ る。電子は，軸方向の磁場と径方向の電場に垂直なE × B 方向に回転運動をする。 表 1: 2T の磁場をかけた場合の代表的パラメタ 磁場に垂直面内での運動の時間スケール ∼100 ms その時の電子密度 ∼1013/m3 閉じ込め領域の平均的軸方向長さ 360mm 最大長 600mm 磁場方向の往復運動の周期 ∼3ms その時の電子温度 ∼0.3 eV 電子フィラメントの 1/e 半径 ∼0.15 m その時の電子数 ∼ 109

で電子分布を測っても，金太郎飴状の分布をしていることが期待できることを表す。 お互いに斥力が働く電子がなぜ閉じ込められるのだろうか？ 電子系の正準角運 動量Pθは，次の式で与えられる。 Pθ = N  i pθi = N  i mer 2 iθ˙i+ eB₂0ri2 (1) 電子質量はme，電荷はe，軸方向の一様磁場の強さは B0，i 番目の電子の正準角運 動量はpθi，軌道半径はri，角速度は ˙θiである。右辺第 1 項と第 2 項の比をとると， mer2iθ˙i eB0r2_i 2 = 2 mev⊥ eB0 ri = 2 rL ri (2) ここで，riθ˙i = v⊥とした。ラーマー半径rLは電子分布のスケールriに比べて十分 小さい。すなわち，(1) 式の第 1 項は第 2 項に比べて十分小さく無視できるので，正 準角運動量は近似的に Pθ ≈ N  i e B0 2 r 2 i = const. (3) となる。この式が意味するところは，1 個の電子が飛び出ようとすると，必ず他の 電子が中心方向に寄らなくてはいけないということである。この拘束条件は，電荷 が単一符号の系だから達成され，異符号の電荷が混在する系では，正負ペアでの脱 走ならば，トータルの正準角運動量が不変に保たれてしまう。つまり，正準角運動 量を使っての閉じ込めは不可能である。 先ほど述べたように，電子分布は金太郎飴状の分布になっているはずなので，磁 場に垂直な面で電子分布のどこを切っても，同じ模様が観測出来るはずである。す なわち，電子分布については磁場に垂直な断面内の運動を追跡すれば，ほぼ全体像 を推定できる。2 次元断面内での電子の運動方程式は次のように導かれる。軸方向 の単位ベクトルを ˆz とする。軸方向磁場は以下のとおりである。 B = B0zˆ (4) 電子の運動方程式は， medv dt =−e(E + v × B) (5) である。運動方程式をサイクロトロン運動の周期で平均化 (時間部分＝ 0) すると， v = E × B |B|2 (6) が導かれる。これは，電子が電場E と磁場 B に垂直な方向に運動することを表す。 これが“E × B ドリフト ”と呼ばれる運動である。(4) 式と静電ポテンシャル E = −∇φ (7) を使って，(6) 式を書き直すと， v = 1 B0z × ∇φˆ (8)

(a) mode 2 (b) mode 3 (c) mode 4 図 3: (a) モード 2 (b) モード 3 (c) モード 4 の Kelvin-Helmholtz 不安定性の例 が得られる。これの回転微分をとると渦度が得られる。 ωzz = ∇ × v =ˆ zˆ B0∇ 2_{φ =} en 0B0zˆ (9) ここでポアソン方程式の関係を使った。(9) 式は，渦度が電子の数密度n に比例する ことを表す。また，速度の発散微分から， ∇ · v = ∇ ·  ˆ z × ∇φ B0  = 0 (10) より，磁場に垂直な面内の電子運動は非圧縮となる。すると，電子の数密度に関す る連続の式 ∂n ∂t +v · ∇n = 0 (11) の両辺にe/(0B0) をかけた式は ∂ωz ∂t +v · ∇ωz = 0 (12) となり，非粘性非圧縮の 2 次元オイラー方程式 (に対応した渦度方程式) に書き換え られることがわかる1。すなわち，磁場に垂直な断面内の電子の運動は，2 次元非圧 縮非粘性流体に対するオイラー方程式で記述可能な事がわかり，磁場と電場で流れ 場を高度に制御した流体実験が，純電子プラズマにより可能となることがわかる。 際本研究室の実験装置で観測された Kelvin-Helmholts(Diocotron) 不安定性により 電子分布が塊に分解される様子を図 3 に示す。この図は，円筒容器の軸上に設置さ れた CCD カメラで撮影した電子分布を表す。明るい部分が電子であり，濃淡が電子 密度を表す。モード 2，3，4 の渦塊が準平衡分布として生成され，その後，熱平衡 分布である右端の分布へ遷移する。なお，本実験の観測は破壊計測なので，たとえ ば 6 枚のスナップショットを撮影するためには，7 回実験を行わなくてはいけないこ とを申し添えておく。 1_{学生諸君，以上の式をすべて導出してみよう!}

3

負温度点渦系

3.1

点渦系

対象とする系は，N/2 個の正の点渦と，N/2 個の負の点渦が，半径 R の円形境界 内に閉じ込められている 2 次元点渦系である。i 番目の点渦 (i = 1, 2, · · · N) の位置 ベクトルをri = (xi，yi)，循環を Ωiで表す。渦度は Dirac のデルタ関数の和で表現 される。 ωz(r, t) =  i Ω_{iδ(r − ri(t))} (13) この系の流れ関数，及び速度場は，G(r) をポアソン方程式に対する Green 関数と して， ψ(r, t) =  i Ω_{iG(r − ri(t))} (14) u(r, t) = −ˆz × ∇ψ(r, t) (15) と与えられ，相互に， ωz(r, t)ˆz = −∇2ψ(r, t)ˆz = ∇ × u(r, t) (16) u(r, t)) = ∇ × (ψ(r, t)) (17) という関係で結ばれている2。循環として取り得る値は Ω₀および−Ω₀(Ω₀は正定数) の 2 種類に限定される。系の保存量は，エネルギーH と慣性モーメント I である。 H = − 1 4π N  i N  j=i Ω_iΩ_jln|r_i− r_j| + 1 4π N  i N  j Ω_iΩ_jln|r_i− ¯r_j| − 1 4π N  i N  j Ω_iΩ_jln R |rj|, (18) I = N  i Ω_i|r_i|2 (19) (18) 式の右辺第 3 項は，境界での流れ関数の値をゼロにするために導入した。また， 半径R の円筒境界の効果を取り入れるため，右辺第 2 項に鏡像渦 ¯rjを導入した。 ¯ rj = R 2 |rj|2 (20) 循環の符号が正負逆の鏡像渦が境界外にあることにより，円筒境界に沿った流れ関 数の値を一定にすることができる。すなわち，円筒境界に対して垂直方向の速度成 分は自動的にゼロとなる。 各点渦の運動方程式は，(18) 式のH を用いて， Ω_idxi dt = ∂H ∂yi, (21) Ω_idyi dt = − ∂H ∂xi (22) 2_{学生諸君，検算をしてみよう!}

図 4: 円形境界に閉じこめられた点渦系の時間発展シミュレーションを元にしたベ ンチマーク結果。個数には鏡像渦の個数は含まれていないので，実際には 104個な らば 2× 104個の点渦の積分を行っている。実機スペックは，表 2 参照。 または，具体的に上記の式を書き下し dri dt =− 1 2π N  j=i Ω_j(ri − rj)× ˆz |ri− rj|2 + 1 2π N  j Ω_j(ri− ¯rj)× ˆz |ri− ¯rj|2 (23) と表される。(23) 式右辺は，離散点渦系に対する 2 次元ビオサバール積分である。 境界からの寄与は，(23) 式右辺第 2 項に含まれる鏡像渦のビオサバール積分により 表現される。 この積分は系の粒子数の 2 乗に計算時間が比例するため，数値計算上は「重たい」計 算に分類される。ここをいかに速く計算するかが，数値計算屋としての腕の見せどこ ろな訳だが，我々はこの計算を速くするために，分子動力学専用計算機 MDGRAPE-2，MDGRAPE-3 を使用した。分子動力学専用計算機は，クーロン力やファンデル ワールス力などの 2 体力を高速に計算する仕組みをもっており，ポアソン方程式の グリーン関数にあたる部分をユーザが自由に変更出来るという特徴をもっている。 この特徴により，クーロン力もファンデルワールス力も統一的に扱うことができる。 我々は，この部分をビオサバール積分用に変更して，シミュレーションに使用した。 少し古いデータになるが，専用計算機と Core 2 Quad の速度比較を行った結果を図 4 に示す。最近では，計算機の CPU もかなり速くなり，GPU(Graphic Processing Unit) も目を見張るほどの計算速度となるが，いまだに MDGRAPE-3 の計算速度に 一日の長があると言える。

表 2: ベンチマークに使用した機材のスペック一覧

名称 CPU Memory FSB Clock

Dual Core Core 2 Duo E6750 (2.66GHz) DDR2 800 4GB 1333MHz Quad Core Core 2 Quad Q6600 (2.4GHz) DDR2 800 4GB 1066MHz MDGRAPE-2 Pentium4 2.4GHz DDR 266 512MB 533MHz MDGRAPE-3 Pentium4 660 (3.6GHz) DDR2 667 2GB 800MHz

W(E)

E

図 5: 状態密度とエネルギーの関係

3.2

統計力学的に定義された絶対温度が負となる状態

系の逆温度β は，統計力学的に β = dS dE = d log W (E) dE (24) と定義される。ここで，W (E) は状態密度，E はエネルギーである。通常，状態密度 はエネルギーの上昇とともに単調増加するためβ が負となることはない (図 5)。こ こで，全相空間体積が有限という条件を与える。すると，E → ∞ の極限で W → 0 とならなければならない。すなわち，状態密度は，あるE0でピークとなり，E > E0 でd ln W/dE < 0 となる (図 6)。ただし，状態密度が最大となるエネルギーは一つ と仮定した。このような系では，E > E0において，系の温度は負となる。エネル ギーが上がるほど状態数が少なくなるような系，例えば，エネルギーが最大となる 唯一無二の状態が存在するような系には，負温度状態が現れる。 上記で説明した負温度状態が有限領域に閉じ込められた点渦系に現れうることを 初めて指摘したのは Onsager である [1, 5]。点渦系の運動方程式を，点渦系のハミル トニアンを用いて表すと，(21)，(22) 式のような形式で表現できる。Onsager は，い わゆるハミルトン方程式との類推から，点渦系の相空間は配位空間と一致すると考 えた。ここで，点渦系が閉じ込められている領域が面積A の領域に限定されたとす ると，相空間の微小体積要素dΓ が dΓ = dx1dy1· · · dxNdyN (25)

W(E)

E

E0

β㸺

β㸼

図 6: 負温度状態がある場合の状態密度とエネルギーの関係 なので，相空間体積は， _ dΓ =  dxdy _N = AN (26) となり，確かに有限となる。よって，ある有限領域に閉じ込められた点渦系は，相空 間体積が有限，つまり全状態数は有限となり，負温度状態が現れうることがわかる。 Onsager が提唱した負温度点渦系については，理論的にも，数値的にも，精力的 に調べられてきた。矩形の平衡分布に関する議論は，Joyce らの結果が有名である。 矩形領域について数値的に平衡分布を求め，同符号点渦が凝集する様子を示してい る。また，理論的な考察から，平衡分布は，sinh-Poisson 方程式 ∇2_{ψ = λ sinh(βψ) (27)} に従うことを予言した [6, 7]。ここで，λ は，最も確かな平衡分布を求める際のラグ ランジュの未定乗数から決まる値，ψ は流れ関数である。ほかにも，円形境界系の軸 対称平衡解 [8]，非軸対称平衡解 [9]，点渦系の状態密度の直接サンプリング [10, 11]， 乱流的な視点での統計性質 [12, 13] などの結果がある。 しかし，数値的な結果を見てみると，計算パワーの限界から，あまり大規模な計 算は行われてこなかったことが読み取れる。そこで，我々は，分子動力学専用計算 機を用いて，点渦系の大規模シミュレーションを始めるに至った。

4

シミュレーション結果

我々は，温度によって点渦系がどのように特徴付けられるかに焦点をあてて数多 くのシミュレーションを行ってきた。ここでは，代表的な結果について紹介をする。

4.1

状態密度の直接数値計算

小正準集団に従ったモンテカルロシミュレーションを行った。次に述べる時間発 展シミュレーションと同じ条件の 6724 個の点渦 (正負がそれぞれ 3362 個ずつ) をラ ンダムに配置し，エネルギーE，慣性モーメント I を計算し，座標 (E, I) に点を打

図 7: 直接数値計算で求めた状態密度。左側の軸が慣性モーメント，右側の軸がエネ ルギーを表す。ピークとなるところを境に温度の符号が入れ替わる。 㻝㻜 㻞㻜 㻟㻜 㻠㻜 㻡㻜 㻢㻜 㻣㻜 㻤㻜 㻥㻜 㻜 㻝 㻞 㻟 㻺㼡㼙㼎㼑㼞㻌㼛㼒㻌㻿㼠㼍㼠㼑㼟㻌㻔㼤㻌㻝㻜 㻡 㻕 㻱㼚㼑㼞㼓㼥㻌E 図 8: 直接数値計算で求めた状態密度。I = 0 での断面図。確かにエネルギーが上が ると状態密度が下がる状態が見て取れる。

β<0

β>0

E

図 9: エネルギー (温度) をパラメタとした平衡分布の相図。青い点は循環が負の点 渦，赤い点は正の点渦を表す。 つ，という作業を 108回繰り返し作成した度数分布図が，図 7 である。 正負同数の ため，I = 0 について対称な形となる。一方，E についても，ある E = E0でピー クとなる分布になっており，少なくとも今回のサンプリングによって，状態密度が ピークとなるE0が少なくとも一つあり，E > E0 では絶対温度が負となることが確 認できた。I = 0 の面で切り取った図から，状態密度がピークとなるエネルギーの 値E0は 29.1 と読み取れる。

4.2

温度

(

エネルギー

)

をパラメタにした平衡分布の分類

今回対象としている系はエネルギー保存系のため，点渦の初期分布でエネルギー をコントロールすることにより，系の温度もコントロールすることができる。図 9 に示すシミュレーションでは，人為的に選んだ様々なエネルギーの初期分布からス タートして，(23) 式を用いた長時間に渡る時間発展シミュレーションを行い，時間 漸近的に各エネルギーにおける平衡分布を得た。温度が正の場合，正の点渦と負の 点渦はそれぞれ円筒境界内に一様に分布しているのに対して，系のエネルギーを上 げ温度が負になる領域に入ると，同符号の点渦同士が凝集し，小さな渦塊を作り始 める。更にエネルギーを上げていくと，正負でそれぞれ 1 個の大きな渦塊ができ，エ ネルギーの上昇とともに凝集の度合いが上がることが確認できる。シミュレーショ ンで追跡することは不可能だが，この系の唯一無二のエネルギー最大となる状態 (回 転方向の自由度は除く) は全ての正の点渦と全ての負の点渦がそれぞれ 1 点に凝集 した状態となり，これ以上エネルギーが高い状態は存在しない。実は，このような 状態が存在することも，負温度状態があり得ることの証拠の一つとなる。

図 10: ランダムな分布の渦度が初期条件として与えられる。左側は流れ関数，右側 は渦度をプロット下ものである。(Physica D, 51 (1991) 531-538 より引用) 図 11: 図 10 からスタートした時間発展シミュレーションで得られた平衡解。この平 衡分布は負温度点渦系の平衡解を記述する sinh-Poisson 方程式に合致する。(Physica D, 51 (1991) 531-538 より引用)

5

解析的結果

5.1

背景

従来，点渦系のシミュレーションでは，Leonard が点渦系のシミュレーションに 関する review において指摘したとおり [14]3，粒子数を上げると結果が思わしくな くなることが経験的に知られていた。 一方，Physica D, 51 (1991) 531-538 に Matthaeus らが投稿した論文では，減衰 性 2 次元 Navier-Stokes 系の平衡解と 2 次元点渦系の平衡解が類似している旨，指摘 があった。 図 10 に示すようなランダムな渦度分布から減衰性 2 次元 Navier-Stokes 方程式 ∂ωz ∂t +u · ∇ωz = ν∇ 2_ω z (28) R = ν−1 = 14286 (29) Δt = 1 2048 (30)

3_{Leonard の指摘は以下の通り：”It now appears that using an increased number of point vortices}

of decreased strength will not yield a converged solution. ... Ironically, best results with the point vortex method often are achieved by using only a few vortices with a diffusive time integration scheme.”

図 12: 渦度が時間発展していく様子の鳥瞰図。正負二つに分離した渦度が見て取れ る。(Physica D, 51 (1991) 531-538 より引用) で時間発展を追跡すると，図 11 や図 12 に見られるように，正の渦と負の渦が 1 個 ずつできあがる様子が見て取れる。この 2 極構造の分布は，先ほど触れた負温度点 渦系の平衡解を与える sinh-Poisson 方程式 [7] によりうまく表現できることが報告 された。 以上，Leonard の指摘と Matthaeus らの報告から，我々は，点渦系というのは非 粘性 2 次元 Euler 方程式の解になっているが，実はこの粒子性に起因するような拡 散項が点渦系には実質的に導入されているのではないか，という予想を行った。こ の予想に基づく解析を行った結果を次に示す。

5.2

ミクロ系とマクロ系

2 次元オイラー方程式 ∂ωz(r, t) ∂t +u · ∇ωz(r, t) = 0 (31) は，ミクロな粒子解である点渦解 (13) を持つ。 ∂ ∂tωz(r, t) = ∂ ∂t   i Ω_{iδ(r − ri(t))}  = − i Ω_i  ∂ ∂tri(t)  · ∇δ(r − ri(t)) = −u(r, t) · ∇ω_z(r, t) (32)

我々は，点渦解をみたすオイラー方程式は，実はミクロな粒子方程式であり，それ がたまたまマクロな流体方程式と同じ形をしているだけなのではないかと考えた。 この点渦解を有するミクロな粒子方程式に対し，何らかの平均操作を施すことによ り得られるマクロな流体方程式には，実は，ミクロな点渦という粒子同士の衝突に 類する効果に由来する粘性が現れるだろうと予想した。 ∂ωz

∂t +u · ∇ωz = viscous term (33)

ここで，· は平均操作を表す演算子である。 これ以後，ミクロな物理量とマクロな物理量を区別するため，ミクロな物理量に は，ˆ· を付すことにする。すなわち，ミクロな渦度は，ˆωz(r, t) と表される。マクロ な渦度は，ミクロな渦度の平均 ωz(r, t) ≡ ˆωz(r, t)SE (34) と定義される。ここで，平均操作を表す演算子の添え字 S，E は，それぞれ空間平 均，アンサンブル平均を表し，空間平均は，微小面積要素 Λ に対して， ˆωz(r, t)S = _|Λ|1  Λ(r)dr _ωˆ z(r, t) (35) と定義される。なお，以降，基本的にrは空間積分の積分変数として用いることに する。一方，ミクロな渦度は，マクロな渦度＋揺らぎ ˆ ωz(r, t) =  i Ω_{iδ(r − ri(t))} = ωz(r, t) + δωz(r, t) (36) と定義される。ミクロな渦度に併せて，ミクロな流れ関数，速度場も次のように表 され，それぞれ，マクロな平均部分と揺らぎから成る。 ˆ ψ(r, t) =  i Ω_{iG(r − ri(t))} (37) ˆ u(r, t) = −ˆz × ∇ ˆψ(r, t) (38)

5.3

拡散係数の導出

出発点は，点渦解を厳密解として持つミクロな渦度方程式である。 ∂ ˆωz ∂t + ˆu · ∇ˆωz = 0 (39) この式に，(36) 式で定義される揺らぎを含んだミクロな渦度を代入する。すると， ∂ ∂tωz(r, t) + ∇ · [u(r, t)ωz(r, t)] = −∇ · δu(r, t)δωz(r, t)SE (40) を得る。左辺は全てマクロな量である。一方，右辺には揺らぎの積の平均が残ると 考える。この項が，粒子性 (ミクロ性) に由来する拡散項であると考える。

δωz(r, t) の具体的表式を得るため，(39) 式に (36) 式で定義される揺らぎを含んだ ミクロな渦度を代入し，揺らぎに関して 1 次の項を集めた線形化方程式を用いる。 ∂ ∂tδωz(r, t) + u(r, t) · ∇δωz(r, t) = −δu(r, t) · ∇ωz(r, t) (41) この式の左辺第 2 項に現れるu(r, t)，および右辺に現れる ω_z(r, t) は，マクロな量 であることに注意する。マクロな量は，ミクロなスケールでは変化しない定数と見 なすことができると考えると，上式は積分可能で， δωz(r, t) =  _t −∞dτ δu (r − (t − τ )u, τ ) · ∇ωz(r, t) (42) を得る。ここで，δωz(r, t = −∞) = 0 を仮定した。この結果を (40) 式右辺に代入す ると， −∇ · δu(r, t)δωz(r, t)SE =−∇ · (↔η · ∇ωz) ↔ η = t

−∞dτ δu(r, t)δu(r − (t − τ )u, τ )SE (43)

となり，拡散係数に対応するテンソルの具体的表式が得られる。これが粒子性に由 来する拡散項の表式である。よく知られた久保公式には時間相関しか含まれていな いが，我々の結果には流れに伴う位置のずれを考慮した相関が含まれていることが 特徴となっている。

6

まとめ

2 次元点渦系について，温度によってその動力学的特性を特徴付けることを目標 とした研究である。シミュレーションでは，状態密度の直接数値計算から，実際に 負温度状態があり得ることを示した。また，負温度において同符号の点渦同士が凝 集したエネルギーの高い状態が平衡状態となることを示した。解析的結果として， 粒子性に起因する粘性が表れることをプラズマの世界でよく知られる Klimontovich formula に従い導出した。今回の我々の結果は，Leonard が述べた経験則が実は正し く，粒子数に依存した粘性が点渦シミュレーションには実効的に導入されているこ とを表している。また，連続系を点渦のような離散的モデルで離散化する場合には， 注意が必要であるということを示唆しているとも言えるだろう。 なお，本結果については，シミュレーションにより数値的に検討することも可能 であり，現在，課題として取り組んでいる最中である。

参考文献

[1] L. Onsager: Nuovo Cimento Suppl. 6 (1949) 279. [2] 際本泰士: 日本物理学会誌 56 (2001) 253.

[3] Y. Kiwamoto, K. Ito, A. Sanpei and A. Mohri: Phys. Rev. Lett. 85 (2000) 3173.

[4] Y. Kiwamoto, N. Hashizume, Y. Soga, J. Aoki and Y. Kawai: Phys. Rev. Lett.

99 (2007) 115002.

[5] G. L. Eyink and K. R. Sreenivasan: Rev. Mod. Phys.78 (2006) 87. [6] G. Joyce and D. Montgomery: J. Plasma Phys. 10 (1973) 107. [7] D. Montgomery and G. Joyce: Phys. Fluids 17 (1974) 1139. [8] S. Kida: J. Phys. Soc. Jpn. 39 (1975) 1395.

[9] R. A. Smith and T. M. O’Neil: Phys. Fluids B 2 (1990) 2961. [10] D. J. Johnson: Phys. Fluids 31 (1988) 1856.

[11] O. B¨uhler: Phys. Fluids 14 (2002) 2139.

[12] T. S. Lundgren and Y. B. Pointin: J. Stat. Phys. 17 (1977) 323. [13] Y. B. Pointin and T. S. Lundgren: Phys. Fluids 19 (1976) 1459. [14] A. Leonard: J. Comput. Phys. 37 (1980) 289.