Occultation diffraction profile with poppy

小惑星による掩蔽とは, 手前にある小惑星が背景の恒星を覆い隠すことによって一時的に恒星が見えなくなる (あるいは暗くなる) 現象を指します. 小惑星そのものの大きさを画像から測ることは大口径の望遠鏡をもちいても困難ですが, 掩蔽観測では小惑星の大きさや形状が減光の光度変化プロファイルに反映されるため, 小口径の望遠鏡であっても, 小惑星の大きさや形状を明らかにすることができます.

ここでは簡単のために恒星は大きさのない点源であり, 十分遠方にあるため小惑星の位置では平面波で近似できると仮定します. 小惑星が平面波の一部を遮蔽するため, 小惑星の後ろには影が落ちます. 観測者がこの影の中に入ったときに掩蔽が観測されます. 掩蔽が観測できる時間は, 影の大きさや形状と影に対する観測者の運動によって変わります. ここで多くのケースでは, 後者は小惑星の位置天文観測から精確に推定することができます. 掩蔽を多くの観測点で捉えることで, 影の大きさと形状, すなわち小惑星の大きさと形状を推定することができます.

Introduction

小惑星が光を遮ったあとの光の伝播が幾何光学で近似できる場合には, 影は明瞭な境界を持つため, 影の形状は小惑星のシルエットに直接対応します. しかし, 小惑星と観測者との距離は ∼1 au のオーダで大きく離れているため, 影の形状は光の回折による影響を受けて変化する可能性があります.

ここでは掩蔽による光の伝播を波動光学的に定式化してみましょう. 光源は \(z\) 軸の負の無限大に位置しているとします. 小惑星は \(z\) 軸方向に厚みのないペラペラの板で代用して, 小惑星の存在する位置を \(z = 0\) と定義します (小惑星面). \(xy\) 平面で小惑星が占めている領域を集合 \(\phi_\text{ast}\) で定義して, 小惑星の形状を表します.

観測者は \((x,y,z) = (0, 0, D)\) に位置しているとします (観測者面). 無限遠からやってきた平面波の一部を波動光学に従って観測者の位置まで伝播させることで, 観測者が見る星の明るさを計算できます. 小惑星面を通過した直後の電場強度を \(E_t(x,y,0)\) とすると, 観測者の位置での電場強度 \(E(0, 0, D)\) は以下のように書けます.

\[ E(0, 0, D) = \frac{1}{i\lambda} \iint E_t(x, y, 0) \frac{\mathrm{e}^{ikr}}{r}\frac{D}{r} \left(1 + \frac{i}{kr} \right) \mathrm{d}x\,\mathrm{d}y, \quad r = \sqrt{x^2 + y^2 + D^2}. \]

\(\lambda\) は観測波長, \(k\) は波数 (\(2\pi/\lambda\)) です. \(E_t(x, y, 0)\) は \((x, y) \not\in \phi_{\text{ast}, t}\) のときにのみ値を持つ関数です. 時刻 \(t\) における小惑星のシルエット \(\phi_{\text{ast}, t}\) が小惑星面を移動することで, 観測者の位置での電場強度 \(E(0, 0, D)\) が変化します. これが掩蔽による減光の光度プロファイルに対応します.

実用的には電場そのものではなく掩蔽が生じていない状態からの比が重要になります. 天体は遥か遠方に存在しているので, 掩蔽が生じていないときの電場強度は小惑星面でも観測者面でも定数 \(E_0\) になるとみなせます.

\[ E_0(0, 0, D) = \frac{1}{i\lambda} \iint E_0 \frac{\mathrm{e}^{ikr}}{r}\frac{D}{r} \left(1 + \frac{i}{kr} \right) \mathrm{d}x\,\mathrm{d}y = E_0 \]

\(E_t(x, y, 0)\) を \(E_0 (1 - A_t(x, y))\) と置き換えて, 電場の相対振幅強度 \(a(0, 0, D) = E(0, 0, D)/E_0\) を計算します.

\[ a(0, 0, D) = 1 - \frac{1}{i\lambda} \iint A_t(x, y) \frac{\mathrm{e}^{ikr}}{r}\frac{D}{r} \left(1 + \frac{i}{kr} \right) \mathrm{d}x\,\mathrm{d}y, \quad r = \sqrt{x^2 + y^2 + D^2}. \]

積分の中では光を透過する項が遮蔽する項へと反転しています. 円形の遮蔽物が作るパターンと円形の開口が作るパターンは明暗が逆にはなりますが基本的には同じ式で記述されることがわかります.

ここで, 観測者から見ると小惑星が移動しているわけですが, 小惑星面と観測者面は平行なので, 相対関係を保ったまま平行移動させても問題ありません. 小惑星の形状変化 (自転) を無視できるのであれば, 小惑星のシルエットを固定して, 以下のように書き直すこともできます.

\[ a(X_t, Y_t, D) = 1 - \frac{1}{i\lambda} \iint A(x, y) \frac{\mathrm{e}^{ikr}}{r}\frac{D}{r} \left(1 + \frac{i}{kr} \right) \mathrm{d}x\,\mathrm{d}y, \quad r = \sqrt{(x - X_t)^2 + (y - Y_t)^2 + D^2}. \]

こちらのケースでは小惑星が観測者面に作る影のパターンの中を, 観測者 \((X_t, Y_t)\) が移動するというイメージになります.

上記の式を直接計算することは大変なため, 多くのケースでは近似してから計算をします. 特に \(\rho^2 = (x - X_t)^2 + (y - Y_t)^2\) が \(D^2\) よりも十分小さい場合には, \(r\) を以下のように展開できます.

\[ r = \sqrt{\rho^2 + D^2} = D\sqrt{1 + \frac{\rho^2}{D^2}} = D \left[ 1 + \frac{1}{2}\frac{\rho^2}{D^2} - \frac{1}{8}\frac{\rho^4}{D^4} + \cdots \right] \]

Fresnel Approximation

Fresnel 近似とは, おおよその項では \(r \approx D\) としつつも, exp の中身では第 2 項まで採用したケースです. この近似を採用することで, 観測者平面における電場振幅を以下のようにかなり簡単な式に整理することができます.

\[ a(X_t, Y_t, D) = 1 - \frac{\mathrm{e}^{ikD}}{i\lambda D} \iint A(x, y) \exp\left[ \frac{ik}{2D}\left( (x-X_t)^2 + (y-Y_t)^2 \right) \right] \mathrm{d}x\,\mathrm{d}y. \]

積分領域を細かく分割して \(A(x, y)\) が定数だとみなせるような小領域の集合を考えると, 以下のように書き換えることもできます. ここで, \(F(t)\) は複素 Fresnel 積分と呼ばれる関数で定義されます.

\[ \begin{aligned} a(X_t, Y_t, D) & = 1 - \frac{\mathrm{e}^{ikD}}{i\lambda D} \sum_{(x_i, y_i) \in \phi_\text{ast}} \int^{x_{1,i}}_{x_{0,i}} \exp\left[ \frac{ik}{2D} (x_i-X_t)^2 \right] \mathrm{d}x \int^{y_{1,i}}_{y_{0,i}} \exp\left[ \frac{ik}{2D} (y_i-Y_t)^2 \right] \mathrm{d}y \\ & = 1 - \frac{\mathrm{e}^{ikD}}{i\lambda D} \sum_{(x_i, y_i) \in \phi_\text{ast}} \left[F(u_{1,i}) - F(u_{0,i})\right] \left[F(v_{1,i}) - F(v_{0,i})\right] \end{aligned} \]

Fresnel 近似が成立するのは \(\rho\) が小さい範囲, 定性的には小惑星面から観測者面に向けた光線が \(z\) 軸となす角が小さいことが求められます. こうした性質から Fresnel 近似は paraxial 近似 (近軸近似)と呼ばれます. 定量的には \(r\) を展開した第 3 項が波長と比較して十分に小さいことが求められます. 小惑星の典型的なサイズを \(l\) とすると \(\rho^2 \simeq 2l^2\) であることが期待されるので, Fresnel 近似が成立するための条件は以下のように整理できます.

\[ \frac{1}{8}\frac{\rho^4}{D^3} \ll \lambda ~~\longrightarrow~~ \frac{1}{2}\frac{l^4}{\lambda D^3} \ll 1 \]

説明の流れから, 一般的な文献とは定数倍程度の違いがあるので注意してください.

Fraunhofer Approximation

さらに距離 \(z\) が大きい領域では小惑星の大きさが距離 \(D\) に対してほぼ無視できるという近似を考えることができます. これを Fraunhofer 近似と呼びます. Fraunhofer 近似では exp の中の \(r\) を以下のように置き換えます. 距離 \(D\) が十分に大きい領域で成り立つ近似であることから, Fraunhofer 近似のことを Far-field 近似 (遠方界近似) とも呼ばれます. 定性的には, 小惑星面から観測者面のある点に向かう光線はすべて並行光線だとみなすことに相当します.

\[ r ~~ \longrightarrow ~~ \underset{\text{Fresnel approx.}} {\underline{ z\left( 1 + \frac{1}{2}\frac{\rho^2}{z^2} \right) }} ~~ \longrightarrow ~~ \underset{\text{Fraunhofer approx.}} {\underline{ z\left( 1 - \frac{xX_t + yY_t}{z^2} + \frac{1}{2}\frac{X^2 + Y^2}{z^2} \right) }} \]

Fraunhofer 近似で得られる電場は以下のとおりです. 観測されるプロファイルを実質的に決めているのは積分の部分ですが, 小惑星の形状 \(A(x,y)\) の Fourier 変換になっていることがわかります.

\[ a(X_t, Y_t, D) = 1 - \frac{\mathrm{e}^{ikD}}{i\lambda D} \mathrm{e}^{ik\frac{X^2+Y^2}{D}} \iint A(x, y) \exp\left[ -\frac{ik}{D}\left( xX_t + yY_t \right) \right] \mathrm{d}x\,\mathrm{d}y. \]

Fourier 変換で定義されているために計算は容易ですが, ひとつ問題があります. Fourier 変換によって影の形状は一意に決まってしまいます. 距離 \(D\) は影の大きさをスケールする因子になるため, 小惑星の大きさと強く縮退してしまいます. この縮退を解くためには外部情報として小惑星までの距離, あるいは影と観測者の相対速度が正確にわかっている必要があります.

Fraunhofer 近似が成立するためには, 無視した項が波長に比べて十分に小さいことが求められます.

\[ \frac{x^2 + y^2}{D} \ll \lambda ~~\longrightarrow~~ 2\frac{l^2}{\lambda D} \ll 1 \]

こちらも一般的な文献とは定数倍程度の違いがあるので注意してください.

Fresnel scale

逆に回折の効果がほとんど効かないケースはどのような場合があるでしょうか. Fresnel 近似で導出した式の一部を整理すると以下のように記述できます.

\[ \exp\left[ \frac{ik}{D}(x - X_t)^2 \right] = \exp\left[ \frac{i\pi}{2}\left(\frac{x - X_t}{f}\right)^2 \right], \quad f = \sqrt{\frac{\lambda D}{2}} \]

\(f\) は Fresnel scale と呼ばれる量で, 回折の効果が現れる典型的な長さに対応しています. Fresnel scale が小惑星の大きさに対してほぼ無視できるような大きさであれば, 小惑星の大きさを推定するときにほぼ無視できます. ここでは暫定的に小惑星の典型的な長さの 1/100 になる場合を基準としてみます.

\[ f \ll \frac{l}{100} ~~\longrightarrow~~ 5{\times}10^3\frac{\lambda D}{l^2} \ll 1 \]

Validity of Approximations

上記の近似が成立する範囲を確認します. 観測波長は 0.55 μm として, 観測者から 0.1, 1, 10 au の距離にある天体のサイズと近似の基準となる値の関係をプロットしました. また Fresnel scale が小惑星の大きさの 1/100 になるケースを想定して, 幾何光学で近似してしまえる基準についても示しました.

小天体のサイズがおおよそ 10 km を超えると幾何光学で近似できる領域に入ってきます. こうした天体では回折の効果を考慮に入れなくても, 小天体のサイズを掩蔽の光度変化プロファイルから精度良く推定することができます.

また, Fresnel 近似は 1,000 km の天体 (∼冥王星) が 0.1 au 以内にあるようなケースを除いてほぼ成立するとみなすことができそうです. 太陽系小天体による掩蔽を考えている限り, Fresnel 近似よりも高い精度での計算が必要になるケースはなさそうです. 距離 10 au で天体のサイズが 1 km を切ってくると Fraunhofer 近似が成立する領域に入ってきます. 距離 1 au でも 100 m 級の小天体であれば, 掩蔽の光度プロファイルは Fraunhofer 近似で表現できそうです.

ここまでをまとめると以下のようになります.

太陽系小天体の掩蔽プロファイルは Fresnel 近似で十分計算できる.
Fresnel 回折の効果が重要になるのはあるサイズ帯 (距離 1 au では ∼0.1–1 m) に限られる.
あるサイズ以下 (距離 1 au では ∼0.1 km) では Fraunhofer 近似で計算できる.

Fresnel diffraction profile

ここでは poppy を使用して掩蔽の減光プロファイルを計算してみます. poppy とは光学系の point-spread function を簡単に計算することができる Python パッケージです.

通常の用途では Fraunhofer 回折による波面伝播を計算しますが, poppy には Fresnel 回折によって波面伝播を計算することもできます. そこで poppy の機能を使用して小惑星による掩蔽の減光プロファイルを計算してみます.

まずは必要なモジュールを読み込みます.

import matplotlib.pyplot as plt
import numpy as np
import poppy
import astropy.units as u

poppy で Fresnel 回折を計算するためには FresnelWavefront で波面を定義して, 中心に遮蔽物を起きます. ここでは CircularAperture で円形開口を定義した後に透過率を反転させて遮蔽物としています. シミュレーションの領域は 4 km × 4 km として, 観測波長は 0.55 μm としました. また, 中心に設置した小惑星の半径は 1 km です.¹

npix = 1024
wf = poppy.FresnelWavefront(
  4000 * u.m, wavelength=0.55 * u.um, npix=npix, oversample=1)
wf *= poppy.InverseTransmission(poppy.CircularAperture(radius=1000))

小惑星面での波面振幅と位相を可視化すると以下のとおりです.

plt.figure(figsize=(10,4))
wf.display('both',colorbar=True, showpadding=True)
plt.suptitle("Asteroid plane (entrance pupil)", fontsize=14)
plt.tight_layout()
plt.show()

Fresnel 回折で波面を観測者面まで伝播させます. ここでは距離 1 au だけ波面を進ませます.

z = 1.0 * u.au
wf.propagate_fresnel(z)

先ほどと同様に波面振幅と位相を可視化します. 回折の効果によって影の縁では波面の振幅が波打っていることがわかります. 波面振幅の画像をある方向に切ったプロファイルを見ることで, 観測される掩蔽の減光プロファイルを得ることができます.

plt.figure(figsize=(10,4))
wf.display('both',colorbar=True)
plt.suptitle(f'Observer plane (at {z})', fontsize=14)
plt.tight_layout()
plt.show()

影の中央を通る線で切ったプロファイルを表示してみます.

fig, ax = plt.subplots(figsize=(10,5))

y, x = wf.coordinates()
ax.plot(x[wf.intensity.shape[1]//2,:], wf.intensity[wf.intensity.shape[1]//2,:])

ax.set_title(f'Cut through profile at z={z:0.1f}')
ax.set_xlabel('Position [m]', fontsize=12)
ax.set_ylabel('Intensity relative to entrance pupil', fontsize=12)
fig.tight_layout()
plt.show()

影の中心から少し離れた位置で観測した場合を考えてみましょう. 観測者は影の中心から \(y\) 方向に 400 m の位置にいるとします. 影が観測者の上を \(x\) 方向に 1500 の速度で通り過ぎたとします. 時刻 \(t\) において観測者がいる位置は以下の式で与えられます.

\[ (x, y) = (1500\, t, -400), \quad t \in [-2.2, +2.2] \]

このラインに沿ってプロファイルを切り出します. ここでは 2 次元マップで隣接するデータから補完するために RectBivariateSpline を使用しました.

実際の観測データはある時間間隔で積分したものになります. ここでは積分時間を \(\Delta{t} = 0.1\) と設定して積分された減光プロファイルを計算しています.

from scipy.interpolate import RectBivariateSpline
from scipy.integrate import cumulative_trapezoid

fit = RectBivariateSpline(x[0, :], y[:, 0], wf.intensity)

t = np.linspace(-2.2, 2.2, 1001)
tx = -2.1 + 0.1 * np.arange(44)
xx = 1500 * t
yy = 0 * t - 400
zz = fit(xx, yy, grid=False)
zt = np.array([0, *cumulative_trapezoid(zz, t)])
z0 = np.interp(tx[0:-1], t, zt)
z1 = np.interp(tx[1:], t, zt)

左右のパネルに観測者面で見た観測者の軌跡と観測される減光プロファイルを図示します.

from astropy.visualization import simple_norm

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

ax = axes[0]
extent = [np.min(x), np.max(x), np.min(y), np.max(y)]
norm = simple_norm(wf.intensity, 'asinh')
ax.imshow(
  wf.intensity, cmap='gray', origin='lower',
  norm=norm, extent=extent)
ax.plot(xx, yy)
ax.set_xlabel('Observer plane X (m)', fontsize=14)
ax.set_ylabel('Observer plane Y (m)', fontsize=14)

ax = axes[1]
ax.plot(t, zz, linewidth=0.5)
ax.scatter(
  tx[:-1] + 0.05, (z1-z0) / 0.1, marker='x', color='C1')
ax.stairs(
  (z1-z0) / 0.1, edges=tx, color='C1')
ax.set_xlabel('Time', fontsize=14)
ax.set_ylabel('Relative Intensity', fontsize=14)

fig.tight_layout()
plt.show()

右パネルの青実線は減光プロファイルをそのまま図示したものです. オレンジ色の実線はある時間で積分した減光プロファイルです. 青実線のプロファイルは激しく変動していますが, 積分した後のプロファイルでは変動がかなりなまされています. 回折プロファイルを測定するためには, 適切な積分時間を選択することが重要になります.

FresnelWavefront のパラメタ oversample は 1 に設定します. 実行すると Warning が出ますが気にせず進めます. ↩