球形粒子による光の散乱

液滴や気泡などの透明球形粒子による光の散乱は，Lorents-Mie散乱理論として知られている．
下記の式は厳密解（解析解）のため，いわゆるMie散乱領域（$I\propto d^2$）だけではなく，より小径な Rayleigh散乱領域を含めて，任意の粒径，（複素数を含む）任意の屈折率に対して計算ができる．

散乱光強度分布の計算

散乱光強度の角度依存性の計算に際しては，粒子直径 $d$ の代わりに，真空中の光波長 $\lambda_{vac}$ に対して規格化した粒径パラメータ $x$ を用いる．

\begin{align} x = \frac{\pi d}{\lambda_{vac}/m_2} \label{eqn:x_def} \end{align}

$m_2$ を球体外部の周囲媒体の屈折率とする．
$m_1$ を球体内の屈折率とし，球体内部で光の吸収がある場合は，$m_1$ は複素数となる．

以下の式では，相対屈折率 $m = m_1/m_2$ を用いる． これは複素数になる可能性がある．$m\in \mathbb C$

散乱角 $\theta$ に対して，散乱光強度は以下の2種類を考慮する必要がある．
入射光，散乱光の光軸をともに含む面を，散乱断面 "plane of reference"と定義する．

1. plane of referenceに対して，垂直な偏光による散乱光強度は \begin{align} I = \frac{i_1}{k^2r^2}I_0, \ \ i_1 = |S_1(\theta)|^2 \notag \end{align}
2. plane of referenceに対して，平行な偏光による散乱光強度は \begin{align} I = \frac{i_2}{k^2r^2}I_0, \ \ i_2 = |S_2(\theta)|^2 \notag \end{align}
3. 円偏光の光強度は \begin{align} I = \frac{\frac{1}{2}(i_1+i_2)}{k^2r^2}I_0 \notag \end{align}

と表される．(van de Hulst p.35)

$S(\theta)$ は，Amplitude function と呼ばれ，
\begin{align} S_1(\theta) &= \sum_{n=1}^{\infty} \frac{2n+1}{n(n+1)} \left( a_n\pi_n(\cos\theta) + b_n\tau_n(\cos\theta) \right) \label{eqn:S_1_def}\\ S_2(\theta) &= \sum_{n=1}^{\infty} \frac{2n+1}{n(n+1)} \left( b_n\pi_n(\cos\theta) + a_n\tau_n(\cos\theta) \right) \label{eqn:S_2_def} \end{align}

である．(同，p.125)

ここで，係数 $a_n, b_n$ は，Riccati Bessel functions（リカッチ・ベッセル関数）, $\psi, \zeta$ と，その微分 $\psi', \zeta'$ を用いて \begin{align} a_n &= \frac{\psi'_n(mx) \psi_n(x) - m\psi_n(mx) \psi'_n(x)} {\psi_n'(mx) \zeta_n(x) - m\psi_n(mx) \zeta'_n(x) } \label{eqn:a_n_def}\\ b_n &= \frac{m\psi'_n(mx) \psi_n(x) - \psi_n(mx) \psi'_n(x)} {m\psi_n'(mx) \zeta_n(x) - \psi_n(mx) \zeta'_n(x) } \label{eqn:b_n_def} \end{align} とあらわされる．（van de Hulst, p.123の式では， $mx=y$ と表記）

また，$\pi, \tau$ は，Angular function と呼ばれ，Associated Legendre Polynomial（ルジャンドル陪関数）, $P_l^m$ を用いて

\begin{align} \pi_n (\cos\theta) &= \frac{1}{\sin\theta} P_n^1 (\cos\theta) \label{eqn:pi_n_def}\\ \tau_n(\cos\theta) &= \frac{d}{d\theta} P_n^1 (\cos\theta) \label{eqn:tau_n_def} \end{align}

と表される．

様々な関数が登場するが，これらのうち，$\psi_n, \psi'_n$，は複素数の関数となる．

• H.C. van de Hulst, "Light Scattering by Small Particles", Wiley New York (1957)

具体的な計算方法

PC などを用いて実際に $S_1, S_2$ を求めるには，以下の項の計算が必要．

• 係数 $a_n, b_n$ の計算のために，Riccati Bessel functions, $\psi, \zeta$ と，その微分 $\psi', \zeta'$
• Angular functions, $\pi_n, \tau_n$ の計算のために，Associated Legendre Polynomial, $P_l^m, m=1$ とその微分

計算式にはいくつかの特殊関数が出てくるが，漸化式を用いる計算法が数多く提案されている．
また，ここ最近の STL(C++), Boost, Pythonのパッケージでは，特殊関数などを計算できる（C言語の）関数などが利用でき，計算の難易度がぐっと下がった．
ただし，計算条件によっては複素数の計算に対応している関数を用いる必要がある．

Coefficient， $a_n, b_n$

Riccati Bessel functions

$a_n, b_n$に出てくる Riccati-Bessel functions, $\psi_n, \zeta_n$ は，Spherical Bessel functions（球ベッセル関数）, Bessel functions（ベッセル関数）を用いて，以下のように定義される． Wiki(fr) Wolfram MathWorld
（一般には，記号 $S, C$ が用いられるが，ここでは p.123，Debye,1909の表記法を採用．）
$\psi_n$ については，屈折率が複素数になる可能性があるので計算時に注意する必要がある．

以下，$z$は複素数とする．
（$a_n, b_n$ の計算に必要なのは $\psi_n, \zeta_n$ のみだが，全て記載）

\begin{eqnarray} \psi_n(z) &=& S_n(z) &\equiv& z j_n(z) &=& \sqrt{\frac{\pi z}{2}} J_{n+1/2}(z) \label{eqn:RB1}\\ \chi_n(z) &=& C_n(z) &\equiv& -z n_n(z) &=& -\sqrt{\frac{\pi z}{2}} N_{n+1/2}(z) \label{eqn:RB2}\\ \xi_n(z) && &\equiv& zh^{(1)}_{n}(z) &=& \sqrt{\frac{\pi z}{2}}H^{(1)}_{n+1/2}(z) &=& \psi_n(z) - i \chi_n(z) \label{eqn:RB3}\\ \zeta_n(z) && &\equiv& zh^{(2)}_{n}(z) &=& \sqrt{\frac{\pi z}{2}}H^{(2)}_{n+1/2}(z) &=& \psi_n(z) + i \chi_n(z) \label{eqn:RB4}\\ \end{eqnarray}

大文字の $J_n, N_n, H_n$ はそれぞれ，Bessel関数，Neumann関数，Hankel関数．（ = 第1, 2, 3種 Bessel関数）
さらに，$H^{(1)}_n,H^{(2)}_n$ は，第1, 2種 Hankel関数

小文字の $j_n, n_n, h_n$ は，それぞれ，球Bessel関数，球Neumann関数，球Hankel関数．（ = 第1, 2, 3種 球Bessel関数）
さらに，$h^{(1)}_n, h^{(2)}_n$ は，第1, 2種 球Hankel関数．

漸化式による計算法

Spherical Bessel 関数，Riccati-Bessel 関数の漸化式の一般形 \begin{eqnarray} f_{n-1}(z) + f_{n+1}(z) &=& \frac{2n+1}{z}f_n(z) \end{eqnarray} より， \begin{eqnarray} f_{n+1}(z) &=& \frac{2n+1}{z}f_{n }(z) - f_{n-1}(z) \ \ \ n=0,1,2, ... \notag\\ \end{eqnarray} または \begin{eqnarray} f_{n }(z) &=& \frac{2n-1}{z}f_{n-1}(z) - f_{n-2}(z) \ \ \ n=1,2,3, ... \notag \end{eqnarray}

微分については \begin{eqnarray} f'_{n+1}(z) &=& f_{n }(z) - \frac{n+1}{z} f_{n+1}(z) \notag\\ f'_{n}(z) &=& f_{n-1}(z) - \frac{n }{z} f_{n }(z) \notag\\ \end{eqnarray} と表される．

結論

まとめると，以下の漸化式のようになる．

$S_n(z) = \psi_n(z)$

\begin{eqnarray} \psi_{n+1}(z) &=& \frac{2n+1}{z}\psi_{n }(z) - \psi_{n-1}(z) \ \ \ n=0, 1, 2, ... \notag\\ \psi_{n }(z) &=& \frac{2n-1}{z}\psi_{n-1}(z) - \psi_{n-2}(z) \ \ \ n=1, 2, ... \notag\\ \psi_{-1}(z) &=& \cos z \notag\\ \psi_{ 0}(z) &=& \sin z \notag\\ \end{eqnarray}

$C_n(z) = \chi_n$

漸化式は $\psi_n(z)$ と同様．
初項は， \begin{eqnarray} \chi_{-1}(z) &=& -\sin z \notag\\ \chi_{ 0}(z) &=& \cos z \notag \end{eqnarray}

$\zeta_n(z)$

漸化式は $\psi_n(z)$ と同様．
初項は， \begin{eqnarray} \zeta_{-1}(z) &=& \psi_{-1}(z) + i \chi_{-1}(z) = \cos z - i \sin z = \cos(-z) + i \sin(-z) = \exp(-iz) \notag\\ \zeta_0(z) &=& \psi_0(z) + i \chi_0(z) = \sin z + i \cos z = i(\cos(-z) + i \sin(-z)) = i\exp(-iz) \notag\\ \zeta_1(z) &=& \psi_1(z) + i \chi_1(z) = \frac{i(\cos(-z) + i \sin(-z))}{z} - (\cos(-z) + i\sin(-z)) = \frac{i\exp(-iz)}{z} - \exp(-iz)\notag \end{eqnarray}

$\psi'_{n+1}(z)$

微分形については， \begin{eqnarray} \psi'_{n+1}(z) &=& \psi_{n }(z) - \frac{n+1}{z}\psi_{n+1}(z) \ \ \ n=0,1,2, ... \notag\\ \psi'_{n }(z) &=& \psi_{n-1}(z) - \frac{n }{z}\psi_{n }(z) \ \ \ n=1,2,3, ... \notag\\ \notag\\ \psi'_0(z) &=& \cos z \notag\\ \end{eqnarray}

\begin{eqnarray} \zeta'_{n+1}(z) &=& \zeta_{n }(z) - \frac{n+1}{z}\zeta_{n+1}(z) \ \ \ n=0,1,2, ... \notag\\ \zeta'_{n }(z) &=& \zeta_{n-1}(z) - \frac{n }{z}\zeta_{n }(z) \ \ \ n=1,2,3, ... \notag\\ \notag\\ \zeta'_0(z) &=& \cos z - i \sin z = \exp(-iz) \notag\\ \end{eqnarray}

漸化式を用いて $a_n, b_n$ を計算するソースコード
an_bn.cpp

参考文献

• Wikipedia:Bessel functions，ベッセル関数
• Wikipedia Riccati-Bessel function(fr)
• NIST Digital Library of Mathematical Functions:Chapter 10 Bessel Functions
• C.F. Bohren and D.R.Huffman, Absorption and Scattering of Light by Small Particles, Wiley-VCH, 1983, 2012

Angular functionの計算

Angular functionのうち，$\pi_n$ は， \begin{eqnarray} \pi_n (\cos\theta) &= \frac{dP_n(\cos\theta)}{d\cos\theta} \label{eqn:pi_n_def2} \end{eqnarray} また，式\eqref{eqn:pi_n_def} より $P_n^1 (\cos\theta)= \sin\theta \cdot \pi_n (\cos\theta)$を用いて \begin{eqnarray} \tau_n(\cos\theta) &=& \frac{d}{d\theta}P_n^1 (\cos\theta) \notag = \cos\theta\cdot\pi_n(\cos\theta) - \sin^2\theta \frac{d\pi_n(\cos\theta)}{d\cos\theta} \label{eqn:tau_n_def2} \end{eqnarray} と表される．

漸化式による計算

ルジャンドル陪関数より，まともに計算することもできるが，漸化式から計算するほうが簡単．
多様な漸化式が提案されているが，コンピュータで数値的に安定に計算できるものを使用する．

$x=\cos\theta$ として，$n=2,3,4,...$ にたいして \begin{eqnarray} \pi_n(x) &=& \frac{2n-1}{n-1} x \pi_{n-1}(x) - \frac{n}{n-1}\pi_{n-2}(x) \notag\\ \tau_n(x) &=& n x \pi_n(x)-(n+1)\pi_{n-1}(x) \notag\\ \notag\\ \pi_0(x) = 0\notag\\ \pi_1(x)=1 \notag \end{eqnarray}

• 再帰関数を用いてAngular functionを計算するソースコード
  pi_tau.cpp
  
• 漸化式を用いてAngular functionを計算するソースコード
  pi_tau_recurr.cpp

計算結果

ルジャンドル陪関数の計算については，こちら

計算例

BHMIEでは，対数微分の計算を行ったり，漸化式の後方計算を行うなどの工夫がされているが，ここでは漸化式を定義通り計算している．

• Mie散乱の計算を行うソースコード
  https://github.com/tatsuya-kawaguchi/ktmie

計算結果の一例を示す．

色の順番が逆になった虹が二本見えるでしょうか．

中項目

中項目

中項目