Radiative Transfer方程式の拡散近似 -

これはレイトレ Advent Calendar 2019(https://qiita.com/advent-calendar/2019/raytracing)の記事です.

はじめに

BSSRDFモデルや subsurface scattering の理論的背景には, Radiative Transfer方程式の拡散近似があります[JMLH01]. ただし[JMLH01]ではこのことを軽く触れているのみで, 詳細は[Ish78]へ委ねられています. 一方で, Radiation Hydrodynamicsの本[Cas04]を見ると, この話題はモーメント定式化[Tho81]を用いてわかりやすく整理されています. そこで, この記事では [Tho81, Cas04]を参考にしてRadiative Transfer方程式のモーメント定式化を行い, [JMLH01]にある拡散方程式を求めます.

Radiative Transfer方程式

位置 ${\mathbf x}$ , 時刻 $t$ で方向 ${\mathbf \omega}$ へ向かう光線の放射輝度を, $L({\mathbf x},\,{\mathbf \omega},t)$ $[\text{W}\cdot\text{m}^{-2}\cdot\text{sr}^{-1}]$ と書きます. このとき, 媒質内を進行する光線の放射輝度の変化は, Radiative Transfer方程式
${ \begin{align} \frac{1}{c}\frac{\partial L({\mathbf x},{\mathbf \omega},t)}{\partial t} + {\mathbf \omega}\cdot\nabla L({\mathbf x},{\mathbf \omega},t) = -\mu_e({\mathbf x}) L({\mathbf x},{\mathbf \omega},t) +\mu_s({\mathbf x}) L_s({\mathbf x},{\mathbf \omega},t) +S({\mathbf x},{\mathbf \omega},t), \label{eq:radiative_transfer}\tag{1} \end{align} }$
${ \begin{align*} \qquad\qquad\qquad\qquad\qquad L_s({\mathbf x},{\mathbf \omega},t) := \int_{S^2} d{\mathbf \omega}'\,p({\mathbf \omega},\,{\mathbf \omega}')L({\mathbf x},{\mathbf \omega}',t) \end{align*} }$

によって与えられます. ここで $\int_{S^2}d{\mathbf \omega}=\int^{2\pi}_{\varphi=0}d\varphi\int^{\pi}_{\theta=0}d\theta\sin\theta$ は単位球面 $S^2$ 上の積分を表します. そして $c$ $[\text{m}\cdot\text{s}^{-1}$ ] は光速, $\mu_s$ $[\text{m}^{-1}$ ], $\mu_a$ $[\text{m}^{-1}$ ], $\mu_e = \mu_s + \mu_a$ は媒質の散乱係数, 吸収係数, 消散係数です. $p({\mathbf \omega},{\mathbf \omega}')$ は散乱の位相関数で,
${ \begin{align} \int_{S^2}d{\mathbf \omega}\,p({\mathbf \omega}\cdot{\mathbf \omega}')=1 \label{eq:phase_function_normalization_condition}\tag{2} \end{align} }$
と規格化されているとします. また散乱は軸対称 $p({\mathbf \omega},{\mathbf \omega}')=p({\mathbf \omega}\cdot{\mathbf \omega}')$ の場合に限るとして, その異方性因子を
${ \begin{align} g := \int_{S^2}d{\mathbf \omega}({\mathbf \omega}\cdot{\mathbf \omega}')p({\mathbf \omega}\cdot{\mathbf \omega}') \label{eq:asymmetry_factor}\tag{3} \end{align} }$
と定義しておきます[JMLH01].
なお, [JMLH01]にある Radiative Transfer方程式は放射輝度の時間変化を無視しているのですが, この記事では残して導出を進めます. ほかの分野と比較したくなった時に便利だからです.
Radiative Transfer方程式(\ref{eq:radiative_transfer})は, 放射輝度 $L$ に関する積分微分方程式となっていることに加えて, 位置 ${\mathbf x}\in\mathbb{R}^3$ , 時間 $t\in\mathbb{R}$ , 方向 ${\mathbf \omega}\in S^2$ の6つの独立変数を持っています. 従って簡単には解けません. そこで次節では, Radiative Transfer方程式を方向 ${\mathbf \omega}$ に依存しない方程式群へ展開する, モーメント定式化を導入します.

モーメント定式化

モーメント展開

単位球面上に値をとる関数 $f({\mathbf \omega})$ はモーメント展開によって,

${ \begin{align} f({\mathbf \omega}) = \frac{1}{4\pi}\int d{\mathbf \omega}f({\mathbf \omega}) + \frac{3}{4\pi}{\mathbf \omega}\cdot\int d{\mathbf \omega}{\mathbf \omega}f({\mathbf \omega}) + \cdots \end{align} }$

と展開できます[Tho81]. この式の第1項, 第2項, $\cdots$ は球面調和関数展開の 0次, 1次, $\cdots$ に対応しています. 展開の各項に現れる量 $\int_{S^2} d{\mathbf \omega}\,f({\mathbf \omega})$ , $\int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}f({\mathbf \omega})$ を一般的に表すために, 関数 $f({\mathbf \omega})$ の $l$ 次のモーメントを
${ \begin{align} \mu_l\left[f({\mathbf \omega})\right] := \int_{S^2}d{\mathbf \omega}\,\underbrace{{\mathbf \omega}\otimes{\mathbf \omega}\otimes\cdots\otimes{\mathbf \omega}}_{l\text{ factors}}f({\mathbf \omega}) \label{eq:define_moment}\tag{4} \end{align} }$
と導入します[JAM+10]. $\otimes$ はテンソル積です.
これを使うと, 放射輝度 $L({\mathbf x},{\mathbf \omega},t)$ は
${ \begin{align} L({\mathbf x},{\mathbf \omega},t) = \frac{1}{4\pi}E({\mathbf x},t) + \frac{3}{4\pi}{\mathbf \omega}\cdot{\mathbf E}_1({\mathbf x},t) + \cdots \label{eq:moment_expansion}\tag{5} \end{align} }$
と展開されます. ここで $L({\mathbf x},{\mathbf \omega},t)$ の 0次のモーメントを scalar irradiance $E({\mathbf x},t)$ , 1次のモーメントを vector irradiance ${\mathbf E}_1({\mathbf x},t)$ と書きました[JMLH01]:
${ \begin{align} E({\mathbf x},t) &:= \mu_0\left[L({\mathbf x},{\mathbf \omega},t)\right] = \int_{S^2}d{\mathbf \omega}\,L({\mathbf x},{\mathbf \omega},t)\qquad[\text{W}\cdot\text{m}^{-2}] \label{eq:scalar_irradiance}\tag{6}\\ {\mathbf E}_1({\mathbf x},t) &:= \mu_1\left[L({\mathbf x},{\mathbf \omega},t)\right] = \int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\,L({\mathbf x},{\mathbf \omega},t)\qquad[\text{W}\cdot\text{m}^{-2}] \label{eq:vector_irradiance}\tag{7} \end{align} }$
後のために, この記事では, 2次のモーメントも tensor irradiance ${\mathbf E}_2({\mathbf x},t)$
${ \begin{align} {\mathbf E}_2({\mathbf x},t) &:= \mu_2\left[L({\mathbf x},{\mathbf \omega},t)\right] = \int_{S^2}d\omega\,{\mathbf \omega}\otimes{\mathbf \omega}\,L({\mathbf x},{\mathbf \omega},t)\qquad[\text{W}\cdot\text{m}^{-2}] \label{eq:tensor_irradiance}\tag{8} \end{align} }$
と書いておきます. また, 自己放射 $S({\mathbf x},{\mathbf \omega},t)$ の0次と1次のモーメントも
${ \begin{align} Q_0({\mathbf x},t) &:= \mu_0\left[S({\mathbf x},{\mathbf \omega},t)\right] = \int_{S^2}d{\mathbf \omega}\,S({\mathbf x},{\mathbf \omega},t) \label{eq:scalar_emission}\tag{9}\\ {\mathbf Q}_1({\mathbf x},t) &:= \mu_1\left[S({\mathbf x},{\mathbf \omega},t)\right] = \int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\,S({\mathbf x},{\mathbf \omega},t) \label{eq:vector_emission}\tag{10} \end{align} }$
と書いておきます.

Radiative Transfer方程式のモーメント定式化

それでは, (\ref{eq:define_moment})を使って Radiative Transfer方程式の0次と1次のモーメント式を求めていきます.

0次のモーメント式

Radiative Transfer方程式(\ref{eq:radiative_transfer})の0次のモーメントをとると,

${ \begin{align*} \frac{1}{c}\frac{\partial}{\partial t}\underbrace{\int_{S^2}d{\mathbf \omega}\,L}_{ = E\text{ (\ref{eq:scalar_irradiance})} } + \underbrace{\int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\cdot\nabla L}_{ = \nabla\cdot{\mathbf E}_1 \text{ (\ref{eq:nablaL_0th_moment})} } = -\mu_t({\mathbf x})\underbrace{\int_{S^2}d{\mathbf \omega}\,L}_{=E\text{ (\ref{eq:scalar_irradiance})}} + \mu_s({\mathbf x})\underbrace{\int_{S^2}d{\mathbf \omega}\,L_s}_{=E\text{ (\ref{eq:pL_0th_moment})}} + \underbrace{\int_{S^2}d{\mathbf \omega}\,S}_{=Q_0\text{ (\ref{eq:scalar_emission})}}. \end{align*} }$

そして下かっこで示した式の整理を行うと,
${ \begin{align} \frac{1}{c}\frac{\partial E({\mathbf x},t)}{\partial t} + \nabla\cdot{\mathbf E}_1({\mathbf x},t) = -\mu_a({\mathbf x})E({\mathbf x},t) + Q_0({\mathbf x},t) \label{eq:rte_0th_moment}\tag{11} \end{align} }$
となります.

1次のモーメント式

同様に, Radiative Transfer方程式(\ref{eq:radiative_transfer})の1次のモーメントをとると,

${ \begin{align*} \frac{1}{c}\frac{\partial}{\partial t} \underbrace{\int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}L}_{ = {\mathbf E}_1\text{ (\ref{eq:vector_irradiance})} } + \underbrace{\int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\left({\mathbf \omega}\cdot\nabla L\right)}_{ = \nabla\cdot{\mathbf E}_2 \text{ (\ref{eq:nablaL_1st_moment})} } = -\mu_t({\mathbf x}) \underbrace{\int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}L}_{ = {\mathbf E}_1\text{ (\ref{eq:vector_irradiance})} } + \mu_s({\mathbf x})\underbrace{\int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}L_s}_{=g{\mathbf E}_1\text{ (\ref{eq:pL_1st_moment})}} + \underbrace{\int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}S}_{=Q_1\text{ (\ref{eq:vector_emission})}} \end{align*} }$

そして下かっこで示した式の整理を行うと,
${ \begin{align} \frac{1}{c}\frac{\partial{\mathbf E}_1({\mathbf x},t)}{\partial t} + \nabla\cdot{\mathbf E}_2({\mathbf x},t) = -(\mu_a({\mathbf x}) + \mu'_s({\mathbf x})){\mathbf E}_1({\mathbf x},t) + {\mathbf Q}_1({\mathbf x},t) \label{eq:rte_1st_moment}\tag{12} \end{align} }$
となります.ここで $\mu'_s:=\mu_s(1-g)$ としました.

Radiative Transfer方程式の拡散近似

前節で得た 0次のモーメント式には1次の量 ${\mathbf E}_1$ が, 1次のモーメント式には2次の量 ${\mathbf E}_2$ が表れます. このままではモーメント式の階層が無限に続いてしまうので, 放射輝度のモーメント展開(\ref{eq:moment_expansion})を1次で打ち切る近似[JMLH01]を行います
${ \begin{align} L({\mathbf x},{\mathbf \omega},t) = \frac{1}{4\pi}E({\mathbf x},t) + \frac{3}{4\pi}{\mathbf \omega}\cdot{\mathbf E}_1({\mathbf x},t) \label{eq:moment_expansion_1st_approximation}\tag{13} \end{align} }$
この近似は, 拡散近似と呼ばれています[Cas04]. これを, 1次のモーメント式(\ref{eq:rte_1st_moment})の第2項へ適用すると, 下記が得られます.

${ \begin{align*} \nabla\cdot{\mathbf E}_2({\mathbf x},t) = \frac{1}{3}\nabla E({\mathbf x},t) \end{align*} }$

詳細

${ \begin{align*} \nabla\cdot{\mathbf E}_2({\mathbf x},t) &= \int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\left[{\mathbf \omega}\cdot\nabla L({\mathbf x},{\mathbf \omega},t)\right]\\ &= \frac{1}{4\pi}\underbrace{\int_{S^2}d{\mathbf \omega}{\mathbf \omega}\left({\mathbf \omega}\cdot\nabla E({\mathbf x},t)\right)}_{ = \frac{4\pi}{3}\nabla E({\mathbf x},t)\text{ (\ref{eq:linear_function_1st_moment})} } +\frac{3}{4\pi}\int_{S^2}d{\mathbf \omega}{\mathbf \omega}\left[ {\mathbf \omega}\cdot\nabla\left({\mathbf \omega}\cdot{\mathbf E}_1({\mathbf x},t)\right) \right]\\ &= \frac{1}{3}\nabla E({\mathbf x},t) +\frac{3}{4\pi}\int_{S^2}d{\mathbf \omega}{\mathbf \omega} [{\mathbf \omega}^T\underbrace{\left(\nabla{\mathbf \omega}\right)}_{=0}{\mathbf E}_1({\mathbf x},t)] +\frac{3}{4\pi} \underbrace{ \int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\left[{\mathbf \omega}^T\left(\nabla{\mathbf E}_1({\mathbf x},t)\right){\mathbf \omega}\right] }_{=0\text{ (\ref{eq:quadratic_form_1st_moment})}}\\ &= \frac{1}{3}\nabla E({\mathbf x},t) \end{align*} }$

$\Box$

さらに, 平均自由行程 $1/(\mu_a + \mu'_s)$ $[\text{m}]$ のあいだで ${\mathbf E}_1$ の時間変化が無視できるくらい小さい $\left|\frac{1}{c}\frac{\partial{\mathbf E}_1}{\partial t}\right| \ll \left|(\mu_a + \mu'_s){\mathbf E}_1\right|$ とすれば, 1次のモーメント式 (\ref{eq:rte_1st_moment}) は

${ \begin{align*} \frac{1}{3}\nabla E({\mathbf x},t) = -(\mu_a({\mathbf x}) + \mu'_s({\mathbf x})){\mathbf E}_1({\mathbf x},t) + {\mathbf Q}_1({\mathbf x},t) \end{align*} }$

となります. これを変形して

${ \begin{align*} {\mathbf E}_1({\mathbf x},t) = 3D({\mathbf x}){\mathbf Q}_1({\mathbf x},t) - D({\mathbf x})\nabla E({\mathbf x},t) ,\quad D({\mathbf x}):=\frac{1}{3(\mu_a({\mathbf x})+\mu'_s({\mathbf x}))} \end{align*} }$

と書き, 0次のモーメント式 (\ref{eq:rte_0th_moment}) へ代入すれば

${ \begin{align*} \frac{1}{c}\frac{\partial E({\mathbf x},t)}{\partial t} - \nabla\cdot\left[D({\mathbf x})\nabla E({\mathbf x},t)\right] = -\mu_a({\mathbf x})E({\mathbf x},t) + Q_0({\mathbf x},t) - 3\nabla\cdot\left[D({\mathbf x}){\mathbf Q}_1({\mathbf x},t)\right] \end{align*} }$

が得られます. 同様にして $E$ の時間変化も無視でき, 媒質が一様である $\mu_a({\mathbf x})\rightarrow\mu_a$ , $\mu_s({\mathbf x})\rightarrow\mu_s$ と仮定すれば,
${ \begin{align} D\nabla^2 E({\mathbf x},t) = \mu_aE({\mathbf x},t) - Q_0 + 3D\nabla\cdot{\mathbf Q}_1({\mathbf x},t) \label{eq:isotropic_diffusion_equation}\tag{14} \end{align} }$
となります. これで, 文献[Ish78, JMLH01, JAM+10]の拡散方程式が求まりました.

まとめ

この記事では, 放射輝度の低次近似(\ref{eq:moment_expansion_1st_approximation})を取ることで, Radiative Transfer方程式(\ref{eq:radiative_transfer})から拡散方程式(\ref{eq:isotropic_diffusion_equation})が得られることを確認しました. ここで書いた内容が, Radiative Transfer方程式の拡散近似を用いた文献[JMLH01, JAM+10]などを読む際の助けになれば幸いです.

付録

この記事内で使った式を挙げます.

内積 ${\mathbf a}\cdot{\mathbf \omega}$ と 2次形式 ${\mathbf \omega}^T{\mathbf A}{\mathbf \omega}$ の 0次と1次のモーメント

[JAM+10] から, これらの関係を引用します.
${ \begin{align} \int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\cdot{\mathbf a}=0 \label{eq:linear_function_0th_moment}\tag{A.1} \end{align} }$
${ \begin{align} \int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}({\mathbf \omega}\cdot{\mathbf a})=\frac{4\pi}{3}{\mathbf a} \label{eq:linear_function_1st_moment}\tag{A.2} \end{align} }$
${ \begin{align} \int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}^T{\mathbf A}{\mathbf \omega} = \frac{4\pi}{3}\text{Tr}({\mathbf A}) \label{eq:quadratic_form_0th_moment}\tag{A.3} \end{align} }$
${ \begin{align} \int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\left({\mathbf \omega}^T{\mathbf A}{\mathbf \omega}\right) = 0 \label{eq:quadratic_form_1st_moment}\tag{A.4} \end{align} }$

位相関数 $p({\mathbf \omega}\cdot{\mathbf \omega}')$ の1次のモーメント

${ \begin{align} \int_{S^2}d{\mathbf \omega}{\mathbf \omega}p({\mathbf \omega}\cdot{\mathbf \omega}') = {\mathbf \omega}'g \label{eq:phase_function_1st_moment}\tag{A.5} \end{align} }$
$p({\mathbf \omega}\cdot{\mathbf \omega}')$ の1次のモーメントをとると,

${ \begin{align*} \mu_1\left[p({\mathbf \omega},{\mathbf \omega}')\right] &= \int_{S^2}d{\mathbf \omega}{\mathbf \omega}p({\mathbf \omega}\cdot{\mathbf \omega}')\\ &= \int_{S^2}d{\mathbf \omega}\left[ {\mathbf \omega}'({\mathbf \omega}\cdot{\mathbf \omega}') + {\mathbf \omega}'\times({\mathbf \omega}\times{\mathbf \omega}') \right]p({\mathbf \omega}\cdot{\mathbf \omega}')\\ &= {\mathbf \omega}'\underbrace{ \int_{S^2}d{\mathbf \omega}({\mathbf \omega}\cdot{\mathbf \omega}')p({\mathbf \omega}\cdot{\mathbf \omega}') }_{=g \text{ (\ref{eq:asymmetry_factor})}} + {\mathbf \omega}'\times\left[ \int_{S^2}d{\mathbf \omega}{\mathbf \omega}p({\mathbf \omega}\cdot{\mathbf \omega}') \right]\times{\mathbf \omega}' \end{align*} }$

初めの等号でベクトル3重積 ${\mathbf \omega}={\mathbf \omega}'({\mathbf \omega}\cdot{\mathbf \omega}') + {\mathbf \omega}'\times({\mathbf \omega}\times{\mathbf \omega}')$ を使いました. また,

${ \begin{align*} {\mathbf \omega}'\times\int_{S^2}d{\mathbf \omega}{\mathbf \omega}p({\mathbf \omega}\cdot{\mathbf \omega}') \underbrace{\sim}_{\text{(\ref{eq:linear_function_1st_moment})}} {\mathbf \omega}'\times{\mathbf \omega}' = 0 \end{align*} }$

を使うと,

${ \begin{align*} \int_{S^2}d{\mathbf \omega}{\mathbf \omega}p({\mathbf \omega}\cdot{\mathbf \omega}') = {\mathbf \omega}'g \end{align*} }$

となります. $\Box$

放射輝度 $\nabla L({\mathbf x},{\mathbf \omega})$ の0次と1次のモーメント

${ \begin{align} \int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\cdot\nabla L({\mathbf x},{\mathbf \omega}) = \nabla\cdot\int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}L \label{eq:nablaL_0th_moment}\tag{A.6} \end{align} }$
${ \begin{align} \int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\left[{\mathbf \omega}\cdot\nabla L({\mathbf x},{\mathbf \omega})\right] = \nabla\cdot\int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\otimes{\mathbf \omega}L({\mathbf x},{\mathbf \omega}) \label{eq:nablaL_1st_moment}\tag{A.7} \end{align} }$
これらは, $\nabla$ が ${\mathbf \omega}$ に作用しないことから得られます.

$L_s({\mathbf x},{\mathbf \omega}) = \int_{S^2}d{\mathbf \omega}'p({\mathbf \omega}\cdot{\mathbf \omega}')L({\mathbf x},{\mathbf \omega}')$ の0次と1次のモーメント

${ \begin{align} \int_{S^2}d{\mathbf \omega}L_s({\mathbf x},{\mathbf \omega}) = E({\mathbf x}) \label{eq:pL_0th_moment}\tag{A.8} \end{align} }$
${ \begin{align} \int_{S^2}d{\mathbf \omega}{\mathbf \omega}L_s({\mathbf x},{\mathbf \omega}) = g{\mathbf E}_1({\mathbf x}) \label{eq:pL_1st_moment}\tag{A.9} \end{align} }$
(\ref{eq:pL_0th_moment}): $L_s$ の0次のモーメントを取ると

${ \begin{align*} \mu_1\left[L_s({\mathbf x},{\mathbf \omega})\right] &= \int_{S^2}d{\mathbf \omega}'L({\mathbf x},{\mathbf \omega}') \underbrace{ \int_{S^2}d{\mathbf \omega}p({\mathbf \omega}\cdot{\mathbf \omega}') }_{=1 \text (\ref{eq:phase_function_normalization_condition})}\\ &= E({\mathbf x}) \end{align*} }$

となります. $\Box$
(\ref{eq:pL_1st_moment}): $L_s$ の1次のモーメントを取ると

${ \begin{align*} \mu_1\left[L_s({\mathbf x},{\mathbf \omega})\right] &= \int_{S^2}d{\mathbf \omega}'L({\mathbf x},{\mathbf \omega}') \underbrace{\int_{S^2}d{\mathbf \omega}{\mathbf \omega}p({\mathbf \omega}\cdot{\mathbf \omega}')}_{ =g{\mathbf \omega}'\text{ (\ref{eq:phase_function_1st_moment})} }\\ &= g\int_{S^2}d{\mathbf \omega}'{\mathbf \omega}'L({\mathbf x},{\mathbf \omega}')\\ &= g{\mathbf E}_1({\mathbf x}) \end{align*} }$

となります. $\Box$

参考文献

[Cas04] John I. Castor. Radiation Hydrodynamics. Cambridge University Press (2004).

[Ish78] Akira Ishimaru. Wave Propagation and Scattering in Random Media, Chapter 9 - diffusion approximation. Academic Press (1978).

[JAM+10] Wenzel Jakob, Adam Arbree, Jonathan T. Moon, Kavita Bala, and Steve Marschner. A radiative transfer framework for rendering materials with anisotropic structure. In SIGGRAPH 2010, ACM, New York, NY, USA, 53:1–53:13 (2010).

[JMLH01] Henrik Wann Jensen, Stephen R. Marschner, Marc Levoy, and Pat Han-rahan. A practical model for subsurface light transport. SIGGRAPH ’01, pages 511–518 (2001).

[Tho81]Kip S. Thorne. Relativistic radiative transfer: moment formalisms. Monthly Notices of the Royal Astronomical Society, 194(2):439–473, 02 (1981).