はじめに
近年では、自動運転やモバイル向けなどで測距センサーが注目されている。その測距センサーとして、dToFであったりFMCWというのがよく話題に上がることが多い。
本記事ではFMCWの解説記事の続きとして、書いているのでできればこちらも一緒に見てほしい。(Piggotの解説論文でFMCW法を学ぶ① -SNRの計算-)
本記事では、前回同様Piggottの論文“Understanding the physics of coherent LiDAR”[Piggott2022]をもとにFMCWの物理の紹介をしていきたいと思う。
| 記事番号 | タイトル概要 | 主な内容 |
|---|---|---|
| 1回目 | Piggottの解説論文でFMCW法を学ぶ① -SNRの計算- | FMCWの基本構造（三角チャープによる距離・速度の同時測定）と、平衡光ヘテロダイン受信機の動作を定式化。ショットノイズ制限下でのSNRが信号光電子数\(N\)に一致することを導出し、太陽光耐性がコヒーレンス時間の短さに起因することを整理 |
| 2回目 (本記事) |
Piggottの解説論文でFMCW法を学ぶ② -集光効率とアンテナ定理- | コヒーレントLiDARが単一モードしか検出できないことを重なり積分で確認（HG/LGモードのクロステーブルをJuliaで可視化）。Siegmanのアンテナ定理 \(A_R\Omega_R\simeq\lambda^2\) を導出し、ローレンツの相反定理から導波路への結合振幅、一様照明下での集光パワーまでを追う |
Understanding the physics of coherent LiDAR (4章)
coherent LiDARシステムはとても基本的な制限がある。それは単一モードの光しか検出できないという点である。このことについて、Piggottの論文ではあまり詳しく解説されていなかったのでSiegmanの論文”The Antenna Properties of Optical Heterodyne Receivers”[siegman1966]を使って解説する。
この単一モードの光しか検出できないというのは、dToFなどのAmplitude modulated LiDARと大きな違いを生む。それは、検出器の開口部の取り扱いである。Amplitude modulated LiDARでは開口部を大きくすることで集光効率を高めることができる。しかし、coherent LiDARでは開口部を大きくすることは必ずしも集光効率を向上させることには直結しない。
理想的なコヒーレントLiDARが拡散性ターゲットから収集する平均光出力\( \langle P\rangle \)は
\begin{equation}\label{eq:1}
\langle P\rangle = \frac{\lambda^2 I}{2\pi} \tag{1}
\end{equation}
ここで、\( I\)はターゲット表面で散乱して帰ってきた光の強度であり、\( \lambda\)は光の波長である。光の波長は\(10^{-6}\)オーダーなので、LiDARシステムによって集光される光は非常に小さい。
多くのcoherent LiDARは、送信光と受光が同じ光路を共有するモノスタティック構成をとなっている。理想的なモノスタティック構成のLiDARの場合、ターゲットから収集する平均光出力\( \langle P\rangle \) は次のようになる。
\begin{equation}\label{eq:2}
\langle P\rangle = \frac{\lambda^2 P_s}{2\pi A_{eff}} \tag{2}
\end{equation}
ここで\(P_s\)は散乱光の合計パワー、\(A_{eff}\)はビームの有効エリアである。この式から読み取れることとして、信号強度はビーム面積に反比例するため、集光することが非常に重要である。
一般的なビーム形状と有効エリアの関係をここに示す。
| Type | Intensity \(I(x,y)\) | Effective area \(A_{\mathrm{eff}}\) |
|---|---|---|
| Flat-top | \(\displaystyle \begin{cases} 1, & x^2+y^2 < r^2 \\ 0, & \text{otherwise} \end{cases}\) | \(\pi r^2\) |
| Gaussian | \(\displaystyle \exp\!\left(-\frac{2(x^2+y^2)}{w^2}\right)\) | \(\pi w^2\) |
Siegmanのアンテナ定理[siegman1966]
前の記事(Piggotの解説論文でFMCW法を学ぶ① -SNRの計算-)で書いた通り、FMCWはヘテロダイン検波を用いた測距システムである。このヘテロダイン検波について解説した論文に”The Antenna Properties of Optical Heterodyne Receivers”[siegman1966]がある。
この節ではPiggottの論文から離れて、Siegmanのアンテナ定理について深堀りしていく。
図1. 光ヘテロダイン受信機の構成
図1に示されるヘテロダイン受信機を考える。光源から分岐された参照用のLO(:Local Oscillator)光電場\( \tilde{u}_0(x,y) e^{j\omega_0 t}\)と測距対象物にあたって帰ってきた信号光電場\( \tilde{u}_1(x,y) e^{j\omega_1 t}\)が存在する。
簡単のため、基準面\(z=0\)を信号波とLO波が合波されている光検出器の近傍にある面をとる。また、信号波とLO波は同じ偏光と仮定する。(後述する式から、同一ではない偏光はキャンセルされることがわかる。)
基準面上の微小面積 \( dA \)に入射する複素スカラー振幅 \( \tilde{u}(x, y, t)\) によって光検出素子の出力端子に生じる微小光電流を、次式のように規格化する。
\begin{equation}\label{eq:3}
I(t) = \int \int \eta (x,y) | \tilde{u}(x, y, t) | ^2 dA \tag{3}
\end{equation}
ここで、\( \eta (x,y)\)は入射光に対する量子効率である。また、電場の合計は
\begin{equation}\label{eq:4}
\tilde{u}(x,y,t) = \tilde{u}_0(x,y) e^{j\omega_0 t} + \tilde{u}_1(x,y) e^{j\omega_1 t} \tag{4}
\end{equation}
また、検出器に実際に流れる全光電流\(I(t)\)は
\begin{equation}\label{eq:5}
I(t) = I_0 + I_1 + \frac{1}{2}[\tilde{I}_{10}e^{j(\omega_1-\omega_0)t} + \tilde{I}_{10}^*e^{-j(\omega_1-\omega_0)t} ] \tag{5}
\end{equation}
ここで、\(I_0\)と\(I_1\)はLO波と信号波によって与えられる光電流であり、
\begin{equation}\label{eq:6}
I_0(t) = \int \int \eta (x,y) | \tilde{u}_0(x, y) | ^2 dA \tag{6}
\end{equation}
\begin{equation}\label{eq:7}
I_1(t) = \int \int \eta (x,y) | \tilde{u}_1(x, y) | ^2 dA \tag{7}
\end{equation}
そして、一般的に差周波光電流の複素位相振幅\(\tilde{I}_{10} \)は
\begin{equation}\label{eq:8}
\frac{1}{2} \tilde{I}_{10} = \int \int \eta (x,y) \tilde{u}_1(x, y)\tilde{u}_0^*(x, y) dA \tag{8}
\end{equation}
\(\tilde{u}_0(x,y) \)と \(\tilde{u}_1(x,y) \) が基準面上で（振幅と位相の両方において）同じ空間変化を持つ場合、これは最適な光混合の状況に帰着し、
\begin{equation}\label{eq:9}
| \tilde{I}_{10}|^2 = 4I_0I_1 \tag{9}
\end{equation}
となる。これは最大感度を与える。ここで入射信号光は\(k_1 = (k_x, k_y, k_z) \)の波数ベクトルを持つ正規化平面波として
\begin{equation}\label{eq:10}
\tilde{u}_1(x,y,z) = \tilde{e}_1 e^{-jk_1r} \tag{10}
\end{equation}
参照面においてこれは
\begin{equation}\label{eq:11}
\tilde{u}_1(x,y,z=0) = \tilde{e}_1 e^{-j(k_xx+k_yy)} \tag{11}
\end{equation}
と表すことができる。もし信号光の波数ベクトルの強度が\( k_1 = |k_1| = \omega_1/c \)、また図1において極座標系をとるとき(\(\theta=0\)はz軸方向、\(\phi\)はxy平面内でのx軸からの角度とする)のとき、次のようになる。
\begin{equation}
k_x = k_1 \sin \theta \cos \phi, k_y = k_1 \sin \theta \sin \phi
\end{equation}
今回の場合の差周波光電流\(\tilde{I}_{10-plane} \)は以下のようになる。(式\eqref{eq:8}とはことなり、信号波が一様な平面波と仮定したのが今回の表現となる。)
\begin{equation}\label{eq:12}
\frac{1}{2} \tilde{I}_{10-plane} = \tilde{e}_1 \int \int \eta (x,y)\tilde{u}_0^*(x, y) e^{-j(k_xx+k_yy)} dxdy \tag{12}
\end{equation}
ここで平均量子効率\( \tilde{ \eta}\)は
\begin{equation}\label{eq:13}
\tilde{\eta} \equiv \frac{\int \int \eta(x,y) | \tilde{u}_0(x, y) | ^2 dxdy }{\int \int | \tilde{u}_0(x, y) | ^2 dxdy } \tag{13}
\end{equation}
で与えられる。この定義と式\eqref{eq:12}を用いて差周波光電流の大きさの2乗は次のように表せる。
\begin{equation}\label{eq:14}
|\tilde{I}_{10}|^2 = 4\tilde{\eta}|\tilde{e}_1|^2A_RI_0 \tag{14}
\end{equation}
この式における \(\tilde{\eta}|\tilde{e}_1|^2\) という積は、直流の光電流密度である。すなわち、複素振幅 \(\tilde{e}_1 \) を持つ平面波が量子効率\(\tilde{\eta} \)の光検出器に入射したときに生じる直流光電流を表す。
この量 \( A_R\) は面積の次元を持ち、実際、光ヘテロダイン受信機がLO波と混合して平面波信号を受信する際の「実効受信面積」すなわち「受信アパーチャ」として解釈できる。
この実効受信アパーチャ \( A_R\) は、通常の電波用アンテナにおける実効アパーチャと非常によく似た概念である。これは信号波の到来方向を表すベクトル\( \Omega\)の関数であり、前の3つの式\eqref{eq:12}\eqref{eq:13}\eqref{eq:14}から次のように与えられる。
\begin{equation}\label{eq:15}
A_R(k_x,k_y) = A_R(\Omega)\\
= \frac{\left|\int \int \eta (x,y)\tilde{u}_0^*(x, y) e^{-j(k_xx+k_yy)} dxdy\right|^2}{\tilde{\eta}^2\int \int | \tilde{u}_0(x, y) | ^2 dxdy }\tag{15}
\end{equation}
式\eqref{eq:15}からは重要なことがわかる。実効的な受信アパーチャは「検出器の感度\(\eta (x,y)\)」と「LO光の強さ・位相\(\tilde{u}_0^*(x, y) \)」を掛け合わせたものをフーリエ変換\(e^{-j(k_xx+k_yy)} \)させたものの2乗となっている。つまり、LO光の強さ・位相が受信アパーチャに影響をするといったことである。
次に、全方向から来た光: \(A_R(\Omega)\) を積分した値を評価したい。立体角の微小角は次のように書ける。
\begin{equation}
dΩ = \sinθ dθ dφ
\end{equation}
しかし、前述の関係式から、次のようにも書けることに注意する。
\begin{equation}
dk_xdk_y = k_1^2 \cos θ \sin θ dθdφ = k_1^2 \cos θ dΩ
\end{equation}
前述したように、光ヘテロダイン受信機はθ=0(z軸方向)を向いた、1つかつ十分に狭い受信角度しか持たないと仮定する。したがって、\(A_R(Ω)\)が無視できない大きさを持つ角度範囲においては、\( \cos θ≈1\) 、\(dΩ≈k_1^{-2}dk_xdk_y\) という近似を用いることができる。これにより、次のように書け直せる。
\begin{equation}\label{eq:16}
∬_{\mathrm{all\_ solid\_ angle}} A_R(Ω)dΩ \simeq (1/k₁)^2 \int \int_{-\infty}^\infty A_R(k_x,k_y)dk_xdk_y \tag{16}
\end{equation}
しかし、\(A_R(k_x,k_y)\) を空間周波数のパワー密度として
\begin{equation}\label{eq:17}
(1/k₁)^2 \int \int_{-\infty}^\infty A_R(k_x,k_y)dk_xdk_y = \left( \frac{2\pi}{k_1}\right)^2 \frac{\int \int \eta^2(x,y) | \tilde{u}_0(x, y) | ^2 dxdy }{\tilde{\eta}^2 \int \int | \tilde{u}_0(x, y) | ^2 dxdy } \tag{17}
\end{equation}
また、
\begin{equation}\label{eq:18}
∬_{\mathrm{all\_ solid\_ angle}} A_R(Ω)dΩ =\frac{\tilde{\eta^2}}{\tilde{\eta}^2}\lambda^2 \tag{18}
\end{equation}
もし量子効率が光検出器全体にわたって一様であれば、\(\frac{\tilde{\eta^2}}{\tilde{\eta}^2} \) は1になる。そして、実際の場合でもこの比は1からそう大きくは離れないと考えられる。したがって、光ヘテロダイン受信機の実効受信アパーチャは次のように近似できる。
\begin{equation}\label{eq:19}
∬ A_R(Ω)dΩ \simeq \lambda^2 \tag{19}
\end{equation}
もしアンテナが、\( Ω_R\)の立体角の視野をもつ単一の主受信ローブを持ち、この視野内では実効アパーチャが\(A_R\)、視野外ではゼロであるとすれば、次のようになる
\begin{equation}\label{eq:20}
A_RΩ_R \simeq \lambda^2 \tag{20}
\end{equation}
この関係がよく言われているアンテナ定理と呼ばれるもので、光の周波数においても同様のことが言える。これはPiggottの論文内での式\eqref{eq:1}と一致する。
単一モード検出についての確認
前述した”単一モードしか検出できない”というのを式で表しているのが式\eqref{eq:8}にあたる。実際に、過去に記事にした高次ガウシアンモードの関数を使ってこれを確認する。(レーザービーム品質(M²)の解説と高次ガウシアンモードのM²)
# %% setup
using PyPlot
using PyCall
using Printf
pyimport("japanize_matplotlib") # 日本語フォント
mcolors = pyimport("matplotlib.colors")
mcm = pyimport("matplotlib.cm")
rc("figure", dpi=130); rc("font", size=11)
n = 300; half = 3.5; w = 1.0
x = range(-half, half, length=n)
X = [xj for xi in x, xj in x] # X[i,j] = x[j]
Y = [xi for xi in x, xj in x] # Y[i,j] = x[i]
R = hypot.(X, Y)
Φ = atan.(Y, X)
# %% modes & overlap
function herm(m, ξ) # 物理学者版エルミート H_m(ξ)
m == 0 && return ones(eltype(ξ), size(ξ))
h0 = ones(eltype(ξ), size(ξ)); h1 = 2 .* ξ
for k in 1:m-1
h0, h1 = h1, 2 .* ξ .* h1 .- 2k .* h0
end
return h1
end
function lague(p, α, x) # 一般化ラゲール L_p^α(x)
p == 0 && return ones(eltype(x), size(x))
l0 = ones(eltype(x), size(x)); l1 = @. 1 + α - x
for k in 1:p-1
l0, l1 = l1, @. ((2k + 1 + α - x) * l1 - (k + α) * l0) / (k + 1)
end
return l1
end
hg(m, k) = herm(m, √2 .* X ./ w) .* herm(k, √2 .* Y ./ w) .* exp.(-(X.^2 .+ Y.^2) ./ w^2)
function lg(p, l) # index=(p,l), l≠0 は複素 (位相 e^{ilφ})
al = abs(l)
rad = (√2 .* R ./ w).^al .* lague(p, al, 2 .* R.^2 ./ w^2) .* exp.(-(R.^2) ./ w^2)
return rad .* exp.(im .* l .* Φ)
end
eta(a, b) = abs2(sum(a .* conj.(b))) / (sum(abs2, a) * sum(abs2, b))
function domain_color(field) # 明度=振幅, 色相=位相 (実配列なら 0/π→シアン/赤)
amp = abs.(field); m = maximum(amp); amp = m == 0 ? amp : amp ./ m
h = (angle.(complex.(field)) .+ π) ./ (2π)
return mcolors.hsv_to_rgb(cat(h, ones(size(h)), amp; dims=3))
end
# %% 5x5 cross-table : 内側=被積分関数, 対角=符号一定→残る, 非対角=正負相殺→∬=0
function cross_table(idx, efun, kind, title)
N = length(idx); M = N + 1
E = [efun(i...) for i in idx]
en = [e ./ sqrt(sum(abs2, e)) for e in E] # 単位ノルム規格化
fig = plt.figure(figsize=(2.15*M, 2.35*M))
for r in 0:N, c in 0:N
ax = fig.add_subplot(M, M, r*M + c + 1)
ax.set_xticks([]); ax.set_yticks([])
if r == 0 && c == 0
ax.axis("off"); continue
elseif r == 0
a, b = idx[c]
ax.imshow(domain_color(E[c]), origin="lower")
ax.set_title("\$E_{$a$b}=\\mathrm{$kind}_{$a$b}\$", fontsize=12)
elseif c == 0
a, b = idx[r]
ax.imshow(domain_color(E[r]), origin="lower")
ax.set_ylabel("\$E_{$a$b}\$", fontsize=13, rotation=0, labelpad=16, va="center")
else
integ = real.(en[r] .* conj.(en[c]))
ax.imshow(domain_color(integ), origin="lower")
ipos = sum(integ[integ .> 0]); ineg = sum(integ[integ .< 0])
diag = r == c
for k in ("top","bottom","left","right")
ax.spines[k].set_color(diag ? "#0a7d00" : "#555")
ax.spines[k].set_linewidth(diag ? 2.2 : 1.0)
end
if diag
ax.text(0.5, -0.10, "\$\\eta=1.0,\\ \\iint=$(round(ipos+ineg, digits=2))\$", transform=ax.transAxes, ha="center", va="top", color="#0a7d00", fontsize=11, fontweight="bold")
ax.text(0.5, -0.35, "符号一定→残る", transform=ax.transAxes, ha="center", va="top", color="#0a7d00", fontsize=9)
else
ax.text(0.5, -0.10, "\$\\iint=0\$ (相殺)", transform=ax.transAxes, ha="center", va="top", color="#333", fontsize=11, fontweight="bold")
ax.text(0.5, -0.3, "\$\\iint_+=$(@sprintf("%+.2f",ipos)),\\ \\iint_-=$(@sprintf("%+.2f",ineg))\$", transform=ax.transAxes, ha="center", va="top", color="#666", fontsize=9)
end
end
end
fig.suptitle(title, fontsize=14, y=0.995)
fig.text(0.5, 0.952, "内側=被積分関数 \$\\mathrm{Re}(E_{mn}E_{lk}^{*})\$（単位ノルム規格化）. 対角=符号一定→\$\\iint\$残る, 非対角=正負が相殺→\$\\iint=0\$", ha="center", fontsize=10, color="#444")
fig.tight_layout(rect=[0, 0.075, 1, 0.93])
sm = mcm.ScalarMappable(cmap="hsv", norm=mcolors.Normalize(-π, π))
cax = fig.add_axes([0.30, 0.032, 0.40, 0.016])
cb = fig.colorbar(sm, cax=cax, orientation="horizontal", ticks=[-π, -π/2, 0, π/2, π])
cb.ax.set_xticklabels(["\$-\\pi\$", "\$-\\pi/2\$", "\$0\$", "\$\\pi/2\$", "\$\\pi\$"])
cb.set_label("位相マップ（明度=振幅）: \$0\\to\$ シアン(正), \$\\pm\\pi\\to\$ 赤(負)", fontsize=10)
display(gcf())
end
# %% Hermite / Laguerre : 0 ≤ m,n,l,k ≤ 1
hg_idx = [(0,0),(0,1),(1,0),(1,1)]
lg_idx = [(0,0),(0,1),(1,0),(1,1)]
cross_table(hg_idx, hg, "HG", "エルミート・ガウス: 重なり積分クロステーブル (\$0\\leq m,n\\leq 1\$)")
cross_table(lg_idx, lg, "LG", "ラゲール・ガウス: 重なり積分クロステーブル (\$0\\leq p,l\\leq 1\$)")
図を確認するとわかるように同じモード(一つ目の電場が\( E_{mn}\), 二つの電場が\( E_{lk}\)としたときに\(m = l\)かつ\( n = k\) )のときは積分値は1になる。逆に別のモード( \(m ≠ l\)または\( n ≠ k\) ))のときは積分値は0になる(逆位相同士が打ち消しあうため)。しかし、この打消しあうというのは単一の受光機(例えば図2では入出力端子)のときであり、画像そのままを受光する場合はこの限りでない。
これに関しては従来のCMOSイメージセンサではGHzのビート信号は受け取れないので、シリコンフォトニクスによる受光とPDアレイのセットで考えられていた。しかし、2025年にAppleの公開した特許ではSPADを用いてこれを回避している(US20250035760A1)。SPADではsub nsオーダーの時間応答を追えるため、SPADならではの回避方法と考えられる。詳細についてはまた新しい記事に乗せようと思う。
平衡光ヘテロダイン受信機におけるsiegman
図1における光ヘテロダイン受信機は、厳密には以前紹介したFMCWにおける平衡光ヘテロダイン受信機とはことなる。
以前のシステムでは、このようにカプラーを使って合波を考えており、その後各PDで受光＆光電流に変換という流れになる(上図参考)。この図の解釈を行う。
①シングルモードファイバーでの構成
…この場合は結論はシンプルでシングルモードファイバーに光をカップリングする時点でモード選択はされており、受信機での光モードについては考える必要はない。
②シングルモードファイバーではなく自由空間またはマルチモードファイバーでの構成
…この場合は様々な光モードが許容される。しかし、上下に存在するPDに光が入射されるときに、式\eqref{eq:4}と同じように電場が加算される形で合波する。その後はsiegmanの論文の議論同様の形になるのでお互い同じモードしか許容されないことがわかる。
構成
図2. 典型的なcoherent LiDARの集光光学系
図2に今回想定する光学系の構成を示す。siegmanの議論とはことなり、ファイバーでの構成となっている。我々のゴールは入力導波路に再入射する光パワーを決めることにある。
入射導波路内のパワー
図3. 典型的なcoherent LiDARの集光光学系の閉曲面について
受信信号を重なり積分の形で表す。レンズと導波路カプラを内包する閉曲面\(S\)を考える。寄与が無視できないのは前面 \(S_1\) と背面 \( S_2\) のみとなるように \(S\) を選ぶ
また、散乱電場と散乱磁場はそれぞれ\(\boldsymbol{E}\)と\(\boldsymbol{H}\)として、「モード場（modal fields）」である\(\boldsymbol{E_m}\)および\(\boldsymbol{H_m} \)を導入する。これらは、望みの導波路モードに対して逆向きに注入した際に生じる場である。
「散乱場」\(\boldsymbol{E}, \boldsymbol{H}\) は、ターゲットに当たって散乱され、LiDARの受光系（レンズ→導波路カプラ）に戻ってくる実際の光の電磁場である。これは我々が受け取りたい「信号」そのものとなっている。
「モード場（modal fields）」\(\boldsymbol{E_m}, \boldsymbol{H_m} \)は、受け取りたい導波路モードそのものの場である。ただし定義の仕方が少し特殊で、論文では「知りたい導波路モードを逆向き（出力から遠ざかる \( −z\) 方向）に注入したときに系全体に生じる場」として定義している。
\begin{equation}\label{eq:21}
\boldsymbol{E_m}(x,y,z)=\mathcal{E}_m(x,y)e^{jβz} \tag{21}
\end{equation}
\begin{equation}\label{eq:22}
\boldsymbol{H_m}(x,y,z)=\mathcal{H}_m(x,y)e^{jβz} \tag{22}
\end{equation}
ここで \(\beta\) は導波路の伝搬定数である。導波路は相反・無損失で、\( \hat{z}\)垂直なある面に対して反射対称と仮定する。横方向のモード場\(\boldsymbol{E_m}\)および\(\boldsymbol{H_m} \)は純実数(本論文appendixB参照)になる。単位パワーに規格化したうえで規格化条件が以下の形になる。
\begin{equation}\label{eq:23}
1 = \frac{1}{2} \Re \left\{ −\int \int_{z=0} \mathcal{E_m}×\mathcal{H_m}^∗⋅\hat{z} dxdy\right\} \tag{23}
\end{equation}
次が重要で散乱場\(\boldsymbol{E}, \boldsymbol{H}\)とが導波路モードに結合してできる振幅 \(a\)を、\(z=0\) 面上のモード重なり積分として
\begin{equation}\label{eq:24}
a = \frac{1}{N} \int \int_{z=0} \left( \boldsymbol{E_m}×\boldsymbol{H} – \boldsymbol{E}×\boldsymbol{H_m}\right)⋅\hat{z} dxdy \tag{24}
\end{equation}
規格化定数\(N\)は
\begin{equation}\label{eq:25}
N = -2 \int \int_{z=0} (\mathcal{E_m}×\mathcal{H_m})⋅\hat{z} dxdy \tag{25}
\end{equation}
とされる。式\eqref{eq:25}に式\eqref{eq:23}を適用することで、\(N=4\)とシンプルに定まる。
規格化定数\(N\)について少し考えてみよう。\(N\)は、式\eqref{eq:24}で定義した重なり積分の値 \(a\)が「導波路モードの振幅」として正しいスケールになるように調整する係数だというのが結論である。
重なり積分\(\int \int_{z=0} \left( \boldsymbol{E_m}×\boldsymbol{H} – \boldsymbol{E}×\boldsymbol{H_m}\right)⋅\hat{z} dxdy \)自体は、モード場の大きさをどう規格化するかで値が変わってしまう。そこで「散乱場 \(\boldsymbol{E}, \boldsymbol{H}\) がちょうどモード場そのものだったら \(a=1 \)になる」ように割り算する分母が \( N\) である。式\eqref{eq:24}の\(N = -2 \int \int_{z=0} (\mathcal{E_m}×\mathcal{H_m})⋅\hat{z} dxdy\) は、モード場を自分自身と重ねた「自己重なり(セルフオーバーラップ)」で、これで割ることで規格化することができる。
\(z=0 \)面が \(S_1\) と一致し \(\hat{z}dxdy=d\boldsymbol{A}\) なので、出力振幅\(a\)は
\begin{equation}\label{eq:26}
a = \frac{1}{4} \int \int_{S_1} \left( \boldsymbol{E_m}×\boldsymbol{H} – \boldsymbol{E}×\boldsymbol{H_m}\right)⋅d\boldsymbol{A} \tag{26}
\end{equation}
となる。
ローレンツの相反定理
式\eqref{eq:26}は、出力振幅 \(a\) を面 \(S_1\)（レンズと導波路カプラを通り抜けた後の面）上の重なり積分で表している。ところがこの積分を実際に評価するには、散乱光がレンズや導波路カプラをどう伝搬するかをモデル化しなければならず、扱いにくい。そこで論文は、散乱面のすぐ手前の面 \(S2\)上の積分に移し替えたい、と考えている。\(S_2\) 上なら余計な光学素子が入っておらず散乱光のコヒーレンス特性以外は無視できるためである。この移し替えを可能にするのがローレンツの相反定理となっている。
Source(今回の場合は光源:Light source)を含まない閉曲面\(S\)について、ローレンツの相反定理は次のようになる。
\begin{equation}\label{eq:27}
0 = \oint_{S} \left( \boldsymbol{E_m}×\boldsymbol{H} – \boldsymbol{E}×\boldsymbol{H_m}\right)⋅d\boldsymbol{S} \tag{27}
\end{equation}
つまり、散乱場\(\boldsymbol{E}, \boldsymbol{H}\)とモード場\(\boldsymbol{E_m}, \boldsymbol{H_m} \)で作ったこの被積分量を閉曲面全体で積分するとゼロになる。ここで前提として前面/背面の場(\(S_1,S_2\))の寄与のみを考えていた。そこで式\eqref{eq:27}は次のように書き換わる。
\begin{equation}\label{eq:28}
\int \int_{S_1} \left( \boldsymbol{E_m}×\boldsymbol{H} – \boldsymbol{E}×\boldsymbol{H_m}\right)⋅d\boldsymbol{S} = -\int \int_{S_2} \left( \boldsymbol{E_m}×\boldsymbol{H} – \boldsymbol{E}×\boldsymbol{H_m}\right)⋅d\boldsymbol{S} \tag{28}
\end{equation}
したがって、式\eqref{eq:26}は扱いやすい散乱面直後の\(S_2\)上の積分で書き直すことができ、
\begin{equation}\label{eq:29}
a = \frac{1}{4} \int \int_{S_2} \left( \boldsymbol{E}×\boldsymbol{H_m} – \boldsymbol{E_m}×\boldsymbol{H}\right)⋅d\boldsymbol{S} \tag{29}
\end{equation}
となる。最後に、この時点で光学系の効率 \(η\) を「モードパワーのうちターゲット面 \(S_2\) に到達する割合」として定義して
\begin{equation}\label{eq:30}
\eta = \frac{1}{2} \Re \left\{ \int \int_{S_2} \boldsymbol{E_m}×\boldsymbol{H_m}^∗⋅d\boldsymbol{S} \right\} \tag{30}
\end{equation}
となる。
平面波近似
図4. 積分面\(S_2\)と散乱面の詳細図
積分面 \(S_2\) を波長スケールで局所的に平面とみなす。局所座標を \(z=0\) が \(S_2\)に一致し、面法線 \(\hat{z}\)が散乱面を向くようにとる。さらに、モード場 \(E_m, H_m\) はこの局所パッチ上でほぼ平行光（コリメート）で、面 \(S_2\) に対して入射角 \(\theta\)で当たっていると仮定する。目標は式\eqref{eq:29}を局所平面領域で書いた式\eqref{eq:29′}を簡単化することとなる。
\begin{equation}\label{eq:29′}
a = \frac{1}{4} \int \int_{z=0} \left( \boldsymbol{E}×\boldsymbol{H_m} – \boldsymbol{E_m}×\boldsymbol{H}\right)⋅\hat{z} dxdy \tag{29′}
\end{equation}
散乱場もモード場も平面波展開を行う。平面波展開すると
\begin{equation}\label{eq:31}
\boldsymbol{E} (\boldsymbol{r}) = \int_{-\infty}^\infty\hat{\boldsymbol{E}}(k_x,k_y) e^{-j\boldsymbol{k} \cdot \boldsymbol{r}}dk_xdk_y \tag{31}
\end{equation}
\begin{equation}\label{eq:32}
\boldsymbol{H} (\boldsymbol{r}) = \frac{1}{Z_0} \int_{-\infty}^\infty \hat{k}×\hat{\boldsymbol{E}}(k_x,k_y) e^{-j\boldsymbol{k} \cdot \boldsymbol{r}}dk_xdk_y \tag{32}
\end{equation}
\begin{equation}\label{eq:33}
\boldsymbol{E_m} (\boldsymbol{r}) = \int_{-\infty}^\infty\hat{\boldsymbol{E}}_m (k_x,k_y) e^{-j\boldsymbol{k}_m \cdot \boldsymbol{r}}dk_xdk_y \tag{33}
\end{equation}
\begin{equation}\label{eq:34}
\boldsymbol{H_m} (\boldsymbol{r}) = \frac{1}{Z_0} \int_{-\infty}^\infty \hat{k}_m×\hat{\boldsymbol{E}}_m(k_x,k_y) e^{-j\boldsymbol{k}_m \cdot \boldsymbol{r}}dk_xdk_y \tag{34}
\end{equation}
ここで\(\boldsymbol{r}=(x,y,z)\)であり、波数ベクトル\(\boldsymbol{k},\boldsymbol{k}_m\)は\((k_x,k_y,\sqrt{4\pi^2/\lambda^2 – k_x^2 -k_y^2})\)である。また、\(\hat{k}\)は\(\boldsymbol{k}/|\boldsymbol{k}|\)で与えられる進行方向の単位ベクトル、\(Z_0\)は自由空間インピーダンスである。ここで注意しなければいけないのは、散乱波は\(-z\)方向に進み、モード場は\(z\)方向に進むことである。
ここでプランシュレルの定理(参考:フーリエ変換で積分を解く。)から
\begin{equation}\label{eq:35}
\int_{-\infty}^\infty f(x)g(x) dx = \int_{-\infty}^\infty \hat{f}(-k)\hat{g}(k) dk \tag{35}
\end{equation}
となり。これを使って、式\eqref{eq:29′}の第一項を平面波展開すると
\begin{equation}\label{eq:36}
\frac{1}{4} \int \int_{S_2} \left( \boldsymbol{E}×\boldsymbol{H_m}\right)⋅\hat{z} dxdy \\
=\frac{1}{4Z_0} \int \int_{-\infty}^\infty \left[ \hat{\boldsymbol{E}}(-k_x,-k_y)×\left( \hat{k}_m(k_x,k_y)×\hat{\boldsymbol{E}}_m(k_x,k_y) \right) \right]⋅\hat{z} dxdy \tag{36}
\end{equation}
また、\(\hat{k}_m(k_x,k_y) = -\hat{k}(-k_x,-k_y)\)という関係があり、かつ平面波振幅\(\hat{E}\)は伝搬方向に垂直（横波）であること。さらに、モードビームが入射角\(\theta\)でほぼ平行光という仮定から、近似\(\hat{k}_m(k_x,k_y)\cdot\hat{z} \simeq \cos\theta\)を使える。これで空間周波数領域の
\begin{equation}\label{eq:37}
\frac{1}{4} \int \int_{z=0} \left( \boldsymbol{E}×\boldsymbol{H_m}\right)⋅\hat{z} dxdy \\
=\frac{1}{4Z_0} \int \int_{-\infty}^\infty \left( \hat{\boldsymbol{E}}(-k_x,-k_y)\cdot \hat{\boldsymbol{E}}_m(k_x,k_y) \right) \cos\theta dk_xdk_y \tag{37}
\end{equation}
が得られ、\(\cos\theta\)が積分の外へ括り出せる形になる。(ここの式変換はかなり端折ったので元論文を確認してください。)これを実空間に戻すと
\begin{equation}\label{eq:38}
\frac{1}{4} \int \int_{z=0} \left( \boldsymbol{E}×\boldsymbol{H_m}\right)⋅\hat{z} dxdy = \frac{1}{4Z_0} \int \int_{-\infty}^\infty \boldsymbol{E}\cdot \boldsymbol{E_m} \cos\theta dxdy \tag{38}
\end{equation}
となり、式\eqref{eq:29′}の第二項も同様に計算すると
\begin{equation}\label{eq:39}
-\frac{1}{4} \int \int_{z=0} \left( \boldsymbol{E_m}×\boldsymbol{H}\right)⋅\hat{z} dxdy = \frac{1}{4Z_0} \int \int_{-\infty}^\infty \boldsymbol{E}\cdot \boldsymbol{E_m} \cos\theta dxdy \tag{39}
\end{equation}
よって、式\eqref{eq:29′}から出力振幅 \(a\) は\(S_2\)上におけるシンプルな積分となる。
\begin{equation}\label{eq:40}
a = \frac{1}{2Z_0} \int \int_{S_2}\boldsymbol{E} \cdot \boldsymbol{E_m} \cos\theta dxdy \tag{40}
\end{equation}
一様照明のもとでの集光パワー
前節の式\eqref{eq:40}で得られた振幅 \(a\) は、あくまで「ある一回の測定における値」である。しかし実際にLiDARが受け取るのは、散乱面のランダムな凹凸に由来するスペックル的なゆらぎを平均した期待パワー \(\langle P \rangle\) である。この節では、その統計的な平均化を行い、集光パワーの一般形を導く。
まず話を簡単にするため、次の仮定を置く。散乱面\(S_2\)は完全に平らであり、モード場はほぼ平行光で入射角 \(\theta\) がほぼ一定とする。モード電場は固定偏光 \(\hat{p}\) を持ち、ゆっくり変化する振幅 \(\psi(x,y)\) で
\begin{equation}\label{eq:41}
\boldsymbol{E_m}(x,y,0) = \psi(x,y)\,e^{-j(k_x x + k_y y)}\,\hat{p} \tag{41}
\end{equation}
と書けるものとする。ここで \(\psi(x,y)\) は集光モードの形状（＝ビームの形）を表す。さらに、散乱面は一様に照らされており、散乱電場 \(\boldsymbol{E}\) の自己相関が位置の差 \(\boldsymbol{r}_1-\boldsymbol{r}_2\) だけに依存すると仮定する。
\begin{equation}\label{eq:42} \langle \boldsymbol{E}(x_1,y_1,0)\,\boldsymbol{E}^\dagger(x_2,y_2,0)\rangle = \mathcal{W}(x_1-x_2,\,y_1-y_2) \tag{42} \end{equation}
ここで \(\langle\cdot\rangle\) は期待値、\(\dagger\) は共役転置であり、\(\mathcal{W}\) は一般には \(3\times3\) のテンソルである。
集光パワーは \(P=|a|^2\) で与えられるので、その期待値 \(\langle P \rangle\) を計算するには式\eqref{eq:40}の \(a\) に含まれる散乱電場 \(\boldsymbol{E}\) を上の自己相関で置き換えればよい。整理のため、スカラー自己相関関数 \(\rho(x,y)\) を
\begin{equation}\label{eq:43}
\rho(x,y) = \frac{\cos\theta}{2Z_0}\,\bar{p}^\dagger\,\mathcal{W}(x,y)\,\bar{p}\;e^{-j(k_x x + k_y y)} \tag{43}
\end{equation}
と定義すると、集光パワーの期待値は
\begin{equation}\label{eq:44}
\langle P \rangle = \frac{\cos\theta}{2Z_0} \int\int\int\int \rho(x_1-x_2,\,y_1-y_2)\,\psi(x_1,y_1)\,\psi^*(x_2,y_2)\,dx_1\,dx_2\,dy_1\,dy_2 \tag{44}
\end{equation}
となる。ここで重要なのが、\(\rho(x,y)\) と \(\psi(x,y)\) の変化の速さの違いである。スカラー自己相関関数 \(\rho(x,y)\) は \(x,y\) がある長さ（コヒーレンス幅）を超えるとほぼゼロになる。一方でモード形状 \(\psi(x,y)\) はビームサイズ程度でしか変化しないので、\(\rho\) に比べてずっと緩やかである。したがって、\(\psi\) が実質的に一定とみなせる長さスケール \(L\) に対して \(\sqrt{x^2+y^2}>L\) で \(\rho\approx0\) とみなせるとき、\(\rho\) をデルタ関数で近似してよい。
\begin{equation}\label{eq:45}
\rho(x,y) \approx P_0\,\delta(x)\,\delta(y) \tag{45}
\end{equation}
ここで \(P_0\) は1つのスペックルあたりの特性パワーであり、\(\rho\) を全空間で積分した値、すなわち自己相関テンソルの空間フーリエ変換 \(\hat{\mathcal{W}}(k_x,k_y)\) を用いて
\begin{equation}\label{eq:46}
P_0 = \int\int \rho(x,y)\,dxdy = \frac{\cos\theta}{2Z_0}\,\bar{p}^\dagger\,\hat{\mathcal{W}}(k_x,k_y)\,\bar{p} \tag{46}
\end{equation}
と書ける。このデルタ関数近似を式\eqref{eq:44}に代入すると、二重の積分がきれいに畳み込まれて、最終的に集光パワーの期待値は非常にシンプルな形にまとまる。
\begin{equation}\label{eq:47}
\langle P \rangle = P_0\,\frac{\cos\theta}{2Z_0} \int\int |\psi(x,y)|^2\,dxdy \tag{47}
\end{equation}
この式は、集光パワーが「1スペックルあたりの特性パワー \(P_0\)」と「ビーム形状の積分 \(\int\int|\psi(x,y)|^2\,dxdy\)」の単純な積に分解された点である。前者はターゲット（散乱面）の性質、後者は光学系・ビームの性質に対応していることがわかる。
本記事では、Siegmanのアンテナ定理 \(A_R\Omega_R\simeq\lambda^2\) と、ローレンツの相反定理を用いた導波路結合振幅の計算から、集光パワーの一般形（式\eqref{eq:47}）まで到達した。集光パワーが「1スペックルあたりの特性パワー \(P_0\)」と「ビーム形状の積分 \(\iint|\psi|^2,dxdy\)」の積に分解される、というのが4章の到達点である。
冒頭の式(1)(2)を再現するには、さらに \(\iint|\psi|^2\) から有効エリア \(A_{\mathrm{eff}}\) を定義し、\(P_0\) をLambertian散乱のモデルから求める必要がある。これは本論文に譲る。
コメント