【伝播計算:第2回】角スペクトル法で見るガウシアン・エルミートガウシアンの伝播

はじめに
角スペクトル法の計算条件
ガウシアンの伝播計算
1. ガウシアンビーム伝播の計算結果
エルミートガウシアンの伝播計算
ガウシアン＋エルミートガウシアンの伝播
まとめ
参考文献

はじめに

　前回記事にてレーザー光の伝播を角スペクトル法で計算できるという話をしました。
本記事では、レンズで集光した際のレーザー光の伝播を角スペクトル法で計算します。ガウシアンビームだけでなく、高次のエルミートガウシアン(Hermite-Gaussian)が伝播してもプロファイルを維持する様子や、複数モードを重ね合わせた場合に強度分布がどのように変化するかを確認します。

前回記事

角スペクトル法の計算条件

波長 : 532nm
ビーム径 : 3mm　(waist:1.5mm)
焦点距離 f = 1 m のレンズで集光して測定
1 ± 0.2 mの範囲におけるビームを角スペクトル法で計算

ガウシアンの伝播計算

まずは一番単純なガウシアンの伝播を計算してみます。
このコードで可視化までいけます。

import numpy as np
import matplotlib.pyplot as plt


# パラメータ設定
wavelength = 532e-9        # 波長 [m]
k = 2 * np.pi / wavelength
z = np.linspace(0.8,1.2,10)# 伝播距離 [m]
beam_radius = 1.5e-3       # 半径1.5mmのビーム
N = 1024                   # グリッド点数（大きめが安定）
L = 10e-3                  # グリッド物理サイズ（10mm）
focal_length = 1.0 #レンズの焦点距離
dx = L / N
x = np.linspace(-L/2, L/2, N)
X, Y = np.meshgrid(x, x)

# 初期ビーム
r2 = X**2 + Y**2
U0_amp = np.exp(-r2 / (beam_radius**2)) 
#レンズ集光による位相
focus = np.exp(-1j * k * r2 / (2*focal_length))

kx = 2 * np.pi * np.fft.fftfreq(N, dx)
ky = kx.copy()
KX, KY = np.meshgrid(kx, ky)
KZ = np.sqrt(k**2 - KX**2 - KY**2 + 0j)

#電場分布
U0 = U0_amp *focus# 振幅 × 位相
U0_ft = np.fft.fft2(U0)

#SA法の計算を随時実施
Pz_stack=[]
for z_prop in z:
    Uz_ft = U0_ft * np.exp(1j * KZ * z_prop)
    Uz = np.fft.ifft2(Uz_ft)
    Pz = np.abs(Uz) ** 2

    Pz_stack.append(Pz.copy())

#グラフ作成
fig, axes = plt.subplots(2,5,figsize=(9, 4), tight_layout=True)
axes = axes.ravel() # 現在の描画対象のサブプロット（アックス）
start=N//2-100
end =N//2+100

for i, idx in enumerate(Pz_stack):
    axes[i].imshow(idx[start:end,start:end])
    axes[i].set_title(f"z = {float(z[i]):.2f}m")
    axes[i].set_aspect("equal")
    axes[i].axis("off")
plt.subplots_adjust(
    left=0.01, right=0.99,
    bottom=0.02, top=0.90,
    wspace=0.002, hspace=0.50
)
plt.show()

import numpy as np
import matplotlib.pyplot as plt


# パラメータ設定
wavelength = 532e-9        # 波長 [m]
k = 2 * np.pi / wavelength
z = np.linspace(0.8,1.2,10)# 伝播距離 [m]
beam_radius = 1.5e-3       # 半径1.5mmのビーム
N = 1024                   # グリッド点数（大きめが安定）
L = 10e-3                  # グリッド物理サイズ（10mm）
focal_length = 1.0 #レンズの焦点距離
dx = L / N
x = np.linspace(-L/2, L/2, N)
X, Y = np.meshgrid(x, x)

# 初期ビーム
r2 = X**2 + Y**2
U0_amp = np.exp(-r2 / (beam_radius**2)) 
#レンズ集光による位相
focus = np.exp(-1j * k * r2 / (2*focal_length))

kx = 2 * np.pi * np.fft.fftfreq(N, dx)
ky = kx.copy()
KX, KY = np.meshgrid(kx, ky)
KZ = np.sqrt(k**2 - KX**2 - KY**2 + 0j)

#電場分布
U0 = U0_amp *focus# 振幅 × 位相
U0_ft = np.fft.fft2(U0)

#SA法の計算を随時実施
Pz_stack=[]
for z_prop in z:
    Uz_ft = U0_ft * np.exp(1j * KZ * z_prop)
    Uz = np.fft.ifft2(Uz_ft)
    Pz = np.abs(Uz) ** 2

    Pz_stack.append(Pz.copy())

#グラフ作成
fig, axes = plt.subplots(2,5,figsize=(9, 4), tight_layout=True)
axes = axes.ravel() # 現在の描画対象のサブプロット（アックス）
start=N//2-100
end =N//2+100

for i, idx in enumerate(Pz_stack):
    axes[i].imshow(idx[start:end,start:end])
    axes[i].set_title(f"z = {float(z[i]):.2f}m")
    axes[i].set_aspect("equal")
    axes[i].axis("off")
plt.subplots_adjust(
    left=0.01, right=0.99,
    bottom=0.02, top=0.90,
    wspace=0.002, hspace=0.50
)
plt.show()

ガウシアンビーム伝播の計算結果

各距離でのビームプロファイルです。画像の1辺が1.95mmです。
焦点距離1mのレンズを使っていて、確かにz = 1mの辺りで一番絞れていることが分かります。
また各z位置でビーム径は変化していますが、強度分布はガウシアンを維持していることが分かります。

エルミートガウシアンの伝播計算

エルミートガウシアンの伝播を計算してみます。
レーザーだとガウシアンビームはTEM00、1次のエルミートガウシアンはTEM01と呼ばれたりします。

コードの変更は簡単でU0_ampを1行変えるだけです。

#エルミートガウシアン
U0_amp = (2 * X / beam_radius) * np.exp(-r2 / beam_radius**2)

#エルミートガウシアン
U0_amp = (2 * X / beam_radius) * np.exp(-r2 / beam_radius**2)

エルミートガウシアンで表されるビームプロファイルを確認してみます。
x方向に2山の形であることが分かります。

#初期ビームプロファイル確認
plt.imshow(np.abs(U0_amp)**2)

#初期ビームプロファイル確認
plt.imshow(np.abs(U0_amp)**2)

エルミートガウシアンの伝播計算は以下のようになります。先ほどと同じく1辺が1.95mmです。
2山のプロファイルにも関わらず、1m付近で集光されていますし2山のプロファイルやピーク強度位置も維持していることが分かります。

ガウシアン＋エルミートガウシアンの伝播

　最後にガウシアンとエルミートガウシアンを足し合わせた時の強度分布と伝播を計算してみます。
X方向とY方向それぞれ2山あるエルミートガウシアン(TEM11)と通常のガウシアンを足し合わせた分布をみてみます。
コードは以下で書きかえればよいです。

#エルミートガウシアン、TEM11+ガウシアン
U0_amp = (4 * X * Y / beam_radius**2) * np.exp(-r2 / beam_radius**2) + np.exp(-r2 / (beam_radius**2))

#エルミートガウシアン、TEM11+ガウシアン
U0_amp = (4 * X * Y / beam_radius**2) * np.exp(-r2 / beam_radius**2) + np.exp(-r2 / (beam_radius**2))

伝播計算結果は下図のようになります。z=1mの集光されるのは通常と同じです。
十字のような形状は維持されていますが、集光点前後で強度が強くなる場所が移り変わっていることが分かります。

まとめ

　角スペクトル法を用いてガウシアンとエルミートガウシアンの伝播を計算しました。
・単一のガウシアン、エルミートガウシアンは強度分布を持ったまま伝播する
・複数のモードを重ね合わせると、伝播に伴い強度分布が変化する。
　
　レーザーの教科書を読んでいると、ガウシアンビームの説明はページ数を割いて説明されるのですが高次モードの説明はすぐ終わってしまいます。私自身も今までは「エルミートガウシアンはマクスウェル方程式の解で基底をなす、高次モードで山の数が増える」程度の理解でした。複数山のあるエルミートガウシアンが分布変えずに伝播していく様は意外でした。
　また任意のビームプロファイルはエルミートガウシアンの和で書けるため、伝播にともない強度分布が変わることが実感持って理解できました。