目次

1. 導入: 地球流体における線形波動論の基礎

2. Wigner変換と擬微分作用素

3. Wigner変換を用いた流体方程式の理論解析

4. シミュレーションデータ分析への利用

1. 導入: 地球流体における線形波動論の基礎

• Fourier変換に基づく分散関係, 偏光関係の解析

• 空間に依存した系におけるWKB近似

はじめに:

• 回転成層流体である, 大気や海洋には多種多様な波動が存在する.

  • 例: 音波, 慣性重力波, Rossby波, Kelvin波.

• 大気/海洋に見られる流体現象は波動力学の観点で説明されることが多い.

  • 代表的なものとしては, 波動-平均流相互作用による惑星規模循環の形成.

  • 波動間相互作用(海洋における内部重力波, 大気Rossby波)も重要なテーマ.

• 相互作用を考える際には, 非線形性が本質的に重要.

  • 大抵の場合は摂動法に基づいて議論される.

  • 線形波動論がすべての基礎になる.

GFDにおける線形波動論の基礎

回転系の浅水方程式:

流速を$(u, v)$, 水面変位を$\eta$, 重力加速度を$g$, コリオリパラメータを$f$, 静止状態での水深を$H$として,

• 境界条件: 無限遠で, $(u, v, \eta) \to 0$

この方程式にはどのような性質の波動が含まれるか?

初等的な考察: $(f, H)$ が定数の場合

指数関数解: $(u, v, \eta) = (\tilde{u}, \tilde{v}, \tilde{\eta}) e^{{\rm i} (kx + ly - \omega t)}$ を方程式に代入, 微分作用素を波数と周波数に置き換える:

• 固有値問題を解いて周波数が定まる: $\omega = 0, \pm \sqrt{f^2 + gH(k^2 + l^2)}$ ... 分散関係

• それぞれの分散関係に対応して, 従属変数の比である$(\tilde{u}, \tilde{v}, \tilde{\eta})$が定まる ... 偏光関係

独立した運動形態である, 慣性重力波モードと定常渦モードが得られる.

なぜ指数関数解を仮定して良い?

根底にあるのはFourier変換:

• それぞれの波数成分は独立に振る舞う.

• あらゆる状態が波数成分の重ね合わせで表現できる.

• 与えられた$(u, v, \eta)$に対し, Fourier変換 $\to$ 波数空間で成分分解 $\to$ Fourier逆変換 という手順を踏むことで, 独立したモード(慣性重力波と定常渦)を分離することができる.

$(f, H)$ が定数でない($x, y$の関数となる)場合

• 波長に比べて $f, H$の変化の空間スケールが十分大きいと仮定.

• 長い空間スケールを導入し, $\mu(x, y) \to (x, y)$, $\mu \ll 1$と置き換える.

  • 空間微分が$(\partial_x, \partial_y) \to \mu (\partial_x, \partial_y)$, と変更される.

• 解の形式を $(u, v, \eta) = (\tilde{u}(x, y), \tilde{v}(x, y), \tilde{\eta}(x, y)) e^{{\rm i}(\theta(x, y) / \mu - \sigma t)}$ と仮定して方程式に代入 (WKB近似), $O(\mu)$の項まで考慮して, $(\tilde{u}, \tilde{v}, \tilde{\eta}, \nabla\theta, \sigma)$の関係式を得る.

• 局所的な波数を$\boldsymbol{k} \equiv \nabla \theta$と定義すると, 分散関係式と偏光関係式が波数と空間位置に依存した形式で, $\sigma = \omega(\boldsymbol{x}, \boldsymbol{k})$, $(\tilde{u}, \tilde{v}, \tilde{\eta}) \equiv \tilde{\bm{u}} = \bm{u}(\boldsymbol{x}, \boldsymbol{k})$と得られる.

指数関数解とWKB近似解の比較

(WKB近似における波数は$\boldsymbol{k} = \nabla \theta$ として定義)

Q1. WKB近似解に対応する波動分離手法はないか？

空間に依存した系に特有の性質

浅水方程式系において, 係数の空間変化$(f(y), H(x,y))$を考慮した場合:
惑星/地形性ベータ効果が生じて渦モードが有限の周波数を得る. (Rossby 波)

例: $df / dy = \beta, H = const.$とした時のRossby波:

浅水方程式で単に指数関数解を代入しただけでは導けない.

WKB近似を使って摂動項をゴリゴリ計算すれば導ける.
背景流の影響等を考慮するとかなり複雑な計算になる.

Q2. 分散関係式の摂動項を見通し良く求める方法はないか？

波動平均流相互作用研究の観点から

エネルギー/擬運動量フラックスの解析. (群速度則: $\bm{F} = A \bm{C}_g$ を満たすように定式化)

• e.g. Takaya & Nakamura (2001), Kinoshita et al. (2010), Kinoshita & Sato (2013)

基本場に対する擾乱について線形化された方程式の理論解析が主な内容.

• WKB近似を前提とするが, 摂動項を無視することで実質的には指数関数解を仮定.

• Kinosita & Sato (2013) は, プリミティブ方程式を$v$のみの式にしてから指数関数解を代入することで, Rossby波の分散関係式を導いている ... ややAdhocな印象.

問題設定をより複雑にして, 背景流の勾配と$f$の変化のスケールが同程度になるような場合(近慣性波やRossby波の分散関係が背景流シアを含む)には, 既存の方法論では太刀打ちできない.

本研究: Wigner変換を用いて変数係数偏微分方程式の理論解析の見通しを良くする

複数種の波動を含む連立方程式を, 各波動成分が従う単独方程式に分解する.

Fourier変換に基づく波動分離では代数行列の固有値分解であったのに対し, Wigner変換に基づく波動分離では表象行列の対角化を行う.

ポイント: 物理量データを固有波動成分へ射影する作用素を導く.

2. Wigner変換と擬微分作用素

• 擬微分作用素と表象

• スター積による計算

• 表象行列に関するアルゴリズム

通常のフーリエ変換では, 微分作用素を波数に置き換える.

Wigner変換は, 一般の線形作用素を関数に変換する.

代表的な例: 量子力学のHamiltonianから古典力学のHamiltonianへの対応,

上記の例では, 一座標$x$を残したまま, 微分作用素$- {\rm i} \hbar \partial_x$を運動量$p$に置き換えている. Wigner変換によって, 変数の自由度は2倍になる: $x \to (x, p)$

変換前の作用素と変換後の関数の関係は, 擬微分作用素とその表象の関係に対応.

擬微分作用素(pseudo-diffential operator)は微分作用素を一般化した概念であり, フーリエ変換と逆変換を利用して, 積分変換の形式で

といった形で定義される. (上式で$\mu$は適当な定数であり, 量子論では$\mu = \hbar$とする)

• 被積分関数に含まれる$f(x, p)$は, 擬微分作用素$\hat{F}$の表象(symbol)と呼ばれる.

• 微分作用素も擬微分作用素の一種: $f = ({\rm i} p / \mu)^n$ $\leftrightarrow$ $\hat{F} = \partial_x^n$.

• 表象を$f = ({\rm i} p / \mu)^\alpha, \alpha \in \mathbb{R}$という形に設定することで、逆ラプラシアンや分数階微分なども擬微分作用素として定義できる.

• Wigner変換とは, 擬微分作用素$\hat{F}$をその表象$f(x, p)$へ対応させる変換である.

Wigner変換の具体的手順:

• 対象となる作用素を積分変換の形式で表現する (Schwartzの核定理に基づく):

• 核関数を相対位置座標についてFourier変換したものが作用素の表象になる:

例: 移流作用素 $u(x) \partial_x$ を考えると, 核関数は $F(x, x') = u(x) \delta'(x - x')$であり, これをフーリエ変換して$f(x,p) = - {\rm i} (p/\mu) u + \partial_x u / 2$を得る. 単に微分作用素を波数へ置き換えた場合に比べて追加の項が生じる.

• 作用素の線形演算(足し引きやスカラー倍)は, 表象の演算に自然に対応する.
  作用素の積は単純には対応付けられない.

積演算の定義が特殊なのは, 擬微分作用素が一般に非可換であるため

• 基本となる２つの作用素を$\hat{p} \equiv - {\rm i} \mu \partial_x, \hat{x} \equiv x$ と定義したとき, $\hat{x}\hat{p} \neq \hat{p}\hat{x}$.

WKB近似では$\mu \ll 1$を仮定. 最も粗い近似では, 作用素を可換とみなす. 作用素の非可換性によって摂動が生じる. このことがRossby波の発生に関与する. ($f$と$\partial_y$が非可換)

擬微分作用素の積に対応する表象の演算: スター積

擬微分作用素の定義式(p15)に立ち返って考察すると

が得られる. 被積分関数の$f$と$g$を$(x, p)$の周りでテイラー展開し, 積分を実行すると,

となる. 結果として, パラメータ$\mu$のベキ級数展開で表現される. 以降は $\mu \ll 1$を仮定し, 高次の項を摂動としてみなして話を進める.

高次元系への拡張

一般の多変量関数 $\bm{\phi}(\bm{x}): \mathbb{R}^d \to \mathbb{C}^n$に対して定義される擬微分作用素:

ここで, $\bm{f}$は$\mathbb{R}^{2d}$で定義された$n \times n$の行列値関数. スター積は

($\bm{m}$と$\bm{n}$には多重指数記法を用いている)

以降は簡略化のため, 表象$\bm{f}(\bm{x}, \bm{p})$に対応する擬微分作用素を$\hat{\bm{F}} = \bm{f}(\hat{\bm{x}}, \hat{\bm{p}})$等と表記する.

Wigner変換の視覚的表現

• 核関数の相対位置についてのFourier変換.

• $\bm{p}$は準局所的に定義された波数.

Wigner変換のご利益

作用素の性質や微積演算を, 表象行列の性質や演算へ置きかえる.

• エルミート作用素 ↔ エルミート行列

• 作用素の分解 ↔ 表象行列の分解

• 逆作用素の計算 ↔ スター積についての逆要素の計算

• 作用素の対角化 ↔ 表象行列の対角化

表象行列の対角化により, 流体方程式が含む波動成分を分離できる.

逆作用素の計算

ある作用素$\hat{\bm{L}}$に対し, 以下の式を満たす逆作用素$\hat{\bm{L}}^{-1}$を求めたいとする:

($\hat{\bm{I}}$は単位行列). ここで $\hat{\bm{L}}$と$\hat{\bm{L}}^{-1}$をWigner変換し, それぞれの表象を$\bm{l}(\bm{x}, \bm{p})$, $\bm{l}^{I}(\bm{x}, \bm{p})$とすると, 解くべき方程式は $\bm{l}$ を既知として,

となる. この式を解くため, 表象を

と展開して代入する.

(逆作用素続き) 方程式の各項を$\mu$のオーダーで分けると,

となる (表象$\bm{a}$に対して$\overline{\bm{a}}^1$は$\mu$に比例した項を取り出すことを意味する). $\mu^0$の式は通常の逆行列の計算により, $\bm{l}^I_0 = \bm{l}^{-1}_0$と解くことができ, 高次の式は

といった形で順次求めることができる. こうして得られた$\bm{l}^I$を表象にもつ擬微分作用素が所期の $\bm{L}^{-1}$になる.

作用素の対角化

与えられた作用素行列$\bm{H} = \bm{h}(\hat{\bm{x}}, \hat{\bm{p}})$を, 対角形式に変換する問題を考える:

表象行列に変換すると,

となる. 通常の行列の固有値問題に似ているが, 積がスター積に置き換わっている.

(対角化続き1)

表象 $\bm{h}, \bm{\omega}, \bm{u}, \bm{u}^I$ を$\mu$について級数展開し, $\bm{u}$と$\bm{u}^I$を以下のように書き換える:

解くべき方程式から$\mu^0$と$\mu^1$に比例した項を取り出すと, それぞれ

となる. 1つめの式は通常の行列の固有値問題であり、これは解けると仮定する.

(対角化続き2)

$\mu^1$の式の性質は, $\mu^0$の方程式の結果によって異なる:

• (1) $\bm{\omega}_0$のすべての対角成分 $(\omega_{0,1}, \omega_{0,2}, \ldots)$ が異なる場合 (縮退がない)

• (2) $\bm{\omega}_0$の対角成分 $(\omega_{0,1}, \omega_{0,2}, \ldots)$ に重複がある場合 (縮退がある)

(1) 縮退のない場合: $\bm{\omega}$と$\bm{u}$の$\mu^1$の補正項が次のように求まる:

(括弧書きの下付き添字は行列内における要素の位置を表す)

(対角化続き3)

(2) 縮退がある場合:

• $\mu^1$の式は一般には計算できない.
  → 周波数が縮退した成分は分離することが出来ない.

特別な状況として縮退した固有値が0の場合 ($\omega_{0, i_1} = \omega_{0, i_2} = \ldots = 0$):

• $\mu^1$の式を計算することができ, 縮退が解ける.
  → 周波数0に縮退していた成分が空間の不均質性により, 固有の周波数を獲得する.

[具体的な計算手順については論文(Onuki 2020)を参照]

回転成層流体では, $\beta$効果により定常渦モードがRossby波として周波数を獲得することに相当する.

ここまでのまとめ:

• Wigner変換は, ベクトル値関数に対する一般の作用素を, 表象行列に変換する.

• 作用素の積は表象行列のスター積に対応する.

• 逆作用素の計算や作用素の対角化等の演算は, 対応する表象行列を級数展開することで, 代数的な計算アルゴリズムに帰着される.

変数係数偏微分方程式の厄介な計算を, Wigner変換を通して見通しよく行うことができる.

3. Wigner変換を用いた流体方程式の理論解析

• 流体方程式の非正準ハミルトン形式

• 線形化方程式の擬微分作用素表現

• 単独方程式への分離

• Wigner分布関数と輸送方程式

流体方程式の一般的表現

非粘性の流体運動を支配する方程式は, 一般に次の非正準ハミルトン形式で書かれる:

• $\bm{\upsilon}(\bm{x}, t)$は状態ベクトルであり, 流速や密度等を表す.

• $\mathcal{H}[\bm{v}]$はハミルトニアン汎関数で, エネルギーと適当なCasimirの和である.

• $\hat{\bm{J}}$は$\bm{\upsilon}$の関数として定義された歪エルミート作用素 ($\hat{\bm{J}}^\dag = -\hat{\bm{J}}$を満たす).

地球流体力学における非正準ハミルトン形式の解説はShepherd (1990) あるいは Salmon (1998) を参照.

方程式の線形化とエネルギーの定義

$\mathcal{H}[\bm{\upsilon}]$が$\bm{\upsilon} = \bm{\Upsilon}$で極小となると仮定, そこからの擾乱を$\bm{\upsilon} = \bm{\Upsilon} + \bm{\upsilon}'$と定義して方程式$(\star)$に代入, $\hat{\bm{A}} \equiv \delta^2 \mathcal{H} / \delta\bm{\upsilon}^2, \hat{\bm{B}} = {\rm i} \hat{\bm{A}} \hat{\bm{J}}(\bm{\Upsilon}) \hat{\bm{A}}$として方程式を線形化すると

が得られる. ここで$\hat{\bm{A}}$は正値エルミート作用素, $\hat{\bm{B}}$はエルミート作用素.

この式から(擬)エネルギーの保存則が導かれる: $\mathcal{E} = (1/2)\int \bm{\upsilon}'^\dag \hat{\bm{A}} \bm{\upsilon}' d\bm{x}$として,

が成立.

作用素の分解

前のページで定義した正値エルミート作用素$\hat{\bm{A}}$を, Wigner変換と表彰行列の演算(論文参照)によって, 以下のように分解する:

ここで, $\hat{\bm{L}}$は可逆作用素であり, $\hat{\bm{L}}^{-1}$が存在. 状態ベクトルを$\bm{\psi} = \hat{\bm{L}}^\dag \bm{\upsilon}' / \sqrt{2}$と再定義すると, 方程式$(\star \star)$は

と変形できる. ここで$\hat{\bm{H}} = \hat{{\bm{L}}}^{-1}\hat{\bm{B}}\hat{\bm{L}}^{\dag-1}$はエルミート作用素. 系のエネルギーは, $\bm{\psi}$のノルムとして, $\mathcal{E} = \int |\bm{\psi}|^2 d\bm{x}$と書かれる.

方程式の分離

Wigner変換によってエルミート作用素に対応する表象行列 $\bm{h}(\bm{x}, \bm{p})$ を定義, 固有値方程式$\bm{u} \star \bm{h} \star \bm{u}^I = \bm{\omega}$を解き対角行列 $\bm{\omega} = diag(\omega_1, \omega_2, \ldots)$ を求める. 状態変数を$\bm{\psi}' = \bm{u}(\hat{\bm{x}}, \hat{\bm{p}}) \bm{\psi}$と再定義すると, 方程式($\star\star\star$)は

と分離できる. ここで, $\omega_i(\hat{\bm{x}}, \hat{\bm{p}})$はそれぞれエルミート作用素であり, $\bm{u}(\hat{\bm{x}}, \hat{\bm{p}})$は, ユニタリ作用素である.

分散関係と偏光関係

• 前のページで得られた表象行列$\bm{\omega}(\bm{x}, \bm{p})$の対角成分$(\omega_1, \omega_2, \ldots)$はそれぞれの波動成分の分散関係を表す.

• 対角化に用いた表象行列$\bm{u}$の各列ベクトルが偏光関係を表す.

• いずれの表象も, 背景パラメータの空間勾配$(\nabla_x f, \nabla_x H)$と波束の大きさの有限性を考慮に入れて, 準局所的に定義された物理量と言える.

• ユニタリ作用素$\bm{u}(\hat{\bm{x}}, \hat{\bm{p}})$は, 状態ベクトル$\bm{\psi}$を固有空間へ写す射影作用素である.

以上の手順により, 空間に依存した系における波動分離を準局所的に行うことができる.

密度作用素とWigner分布関数

ある波動成分の状態ベクトル$\psi(\bm{x})$に対し, 次のように定まる作用素$\hat{W}$を定義する:

($\phi(\bm{x})$は任意の関数). このように定義された$\hat{W}$は, 量子力学では密度作用素と呼ばれる. さらに, 密度作用素をWigner変換して得られる表象$w(\bm{x}, \bm{p})$は, Wigner分布関数と呼ばれる:

Wigner分布関数は物理-波数空間内で定義されたエネルギー密度と解釈でき, 積分値は全エネルギーに一致する: $(1/2\pi\mu)^d\iint w(\bm{x}, \bm{p}) d\bm{x}d\bm{p} = \mathcal{E}$.

輸送方程式

状態ベクトルの時間発展が ${\rm i} \partial_t \psi = \omega (\hat{\bm{x}}, \hat{\bm{p}}) \psi$であるとき, Wigner分布関数の時間発展は

となる. スター積を$\mu$で展開し, $O(\mu^3)$の項を無視すると,

となる. この式は物理-波数空間内におけるエネルギー密度の伝達を記述するものであり, 輸送方程式(あるいは放射伝達方程式など)と呼ばれる.

• 輸送方程式の特性曲線は, $\dot{\bm{x}} = \nabla_p \omega, \ \ \dot{\bm{p}} = - \nabla_x \omega$ で定められる.
  これは通常のRay-tracing 方程式である.

エネルギーフラックス

輸送方程式を物理空間に射影(波数 $\bm{p}$ で積分)すると, エネルギー密度の時間発展が

と得られる. ここで,

である. エネルギーフラックス $=$ エネルギー密度 $\times$ 群速度という形式になっている.

ここまでのまとめ

• 一般の流体方程式に対し, Wigner変換を利用して波動成分を準局所的に分離する手順を構築.

• Wigner分布関数を導入することで, 物理-波数空間内でのエネルギー伝達を記述.

• 群速度に比例した形式のフラックスを計算できることを示した.

4. シミュレーションデータ分析への利用

• 浅水方程式の擬微分作用素表現

• 数値モデルデータから波動成分を分離してエネルギーフラックスを描画

浅水方程式の擬微分作用素表現

流速を$(u, v)$, 水深を$h$, 重力加速度を$g$, コリオリパラメータを$f(y)$, 水底の地形を$D(x, y)$とする. ハミルトニアン(全エネルギー)を

と定義して, 浅水方程式は非正準ハミルトン形式で次のように書ける:

(続き) 定常解として静止状態

を考え, 前述の手順(p31)に沿って擾乱について線形化された方程式を導く:

状態ベクトルを$\bm{\psi} = ( \sqrt{H/2}u', \sqrt{H/2}v', \sqrt{g/2}\eta')$と定義すると, 上の式は

と書ける. (表象行列$\bm{h}$の具体的な形は次のページ)

(続き) 浅水方程式を特徴づけるエルミート作用素の表象行列$\bm{h}(\bm{x}, \bm{p})$:

ここで, $\bm{p} = (p_x, p_y)$, $c = \sqrt{gH}$, $\nabla_x c = (c_x, c_y)$とした.

浅水方程式の対角化

定められた手順(p24-p27)に沿い, 方程式$\bm{u} \star \bm{h}(\bm{x}, \bm{p}) \star \bm{u}^I = \bm{\omega}$を解いて作用素を対角化する. 偏光関係を表す行列$\bm{u}$は

となる. $\bigl($ $p = \sqrt{p_x^2 + p_y^2}$, $\theta = \sqrt{f^2 + c^2 p^2}$とした$\bigr)$

アルゴリズムの実行にはMaximaを利用.

(続き) 分散関係を表す $\bm{\omega} = diag(\omega_{G+}, \omega_{G-}, \omega_R)$は,

となる($\beta \equiv \partial_y f$とした). それぞれ慣性重力波とRossby波の周波数に対応する.

• パラメータ$f, H$の空間勾配を摂動として生じるRossby波の分散関係を正確に導くことができた.

• $\beta$効果により, 慣性重力波は波数ベクトルの方向によって分散関係が変化する.
  (このことは最近 Perez et al. でも議論されている.)

数値モデルデータの分析

• 浅水系の数値シミュレーション.

• 二重周期境界条件の正方形領域で, 水底の形状を右上の図に設定.

• 初期の水面に変位を与えて, 波束の運動を観察する.

境界条件: 二重周期
格子点数: $128 \times 128$
格子幅$\sim$変形半径/$2$

基本的な出力データ

左から: $u, v, \eta$

慣性重力波の速い伝搬と, 地形性Rossby波の遅い伝搬が見られる.

波動分離解析

射影作用素 $\bm{u}(\hat{\bm{x}}, \hat{\bm{p}})$を用いて波動成分を分離する.

• 左: 慣性重力波の信号 $\psi_G$
• 右: Rossby波の信号 $\psi_R$
  (上が実部, 下が虚部)

(Rossby波の信号を解析する際, $\omega_R$が正の部分のみを取り出すフィルタ操作を行っている. 詳細は論文参照)

エネルギー解析

• Wigner分布関数$w(\bm{x}, \bm{p})$を計算, 物理空間に射影してエネルギー密度$E(\bm{x})$とフラックス$\bm{F}(\bm{x})$を得る.　[左:慣性重力波, 右:Rossby波]

エネルギー解析

• ある瞬間のスナップショット.
  それぞれの波の群速度に比例したフラックスを描画することができた.

全体まとめ

• 流体波動を記述する偏微分方程式に対し, Wigner変換を用いることで, 見通しよく方程式の変形操作を行うことができる.

• Wigner変換により得られた表象行列の対角化により, 連立偏微分方程式に含まれる独立した波動成分が準局所的に分離される.

• それぞれの波動の信号からWigner分布関数を定義し, 物理空間へ射影することで, 群速度に比例したフラックスが得られる.

(補足)

一般的な振動モードの分離

• 閉じた領域での固有値解析

  • 離散モードとして得られる.

  • 各モードは領域全体に広がった構造.
    (右: Rossby-Haurwitz 波)

• 今回のテーマ

  • 無限に広い領域.

  • 空間に局在した波束.

(補足) Maximaを用いた表象行列の対角化計算

(続き)