レーザー光の基本ガウシアンモードの導出とGuoy phase shift等の特性

2020年4月18日

URLをコピーしました！

前回の記事で、電磁気学の基礎方程式であるMaxwell方程式から電磁波が導かれることを示しました。

方針

Maxwell 方程式から導いた電磁波の方程式にたいして、近軸近似を行い簡単な電磁波方程式にする。
これにより大きな分布を持った面のような波ではなく、まさにレーザーのような線上の電磁波となる。
上記の近似電磁波方程式をフーリエ変換して、解く。
さらに逆フーリエ変換して、電場(磁場)の分布を明らかにする。
結果、これがGaussian分布となるはずで、これがまさにGaussian beamを意味している。

電磁波方程式の復習

真空中のMaxwell方程式は下記です。

$$
\begin{align}
\nabla\cdot\vec{E}(t, \vec{x}) &= 0 \\
\nabla\cdot\vec{B}(t, \vec{x}) &= 0 \\
\nabla\times\vec{E}(t, \vec{x}) &= -\frac{\partial \vec{B}(t, \vec{x})}{\partial t} \\
\nabla\times\vec{B}(t, \vec{x}) &= \epsilon_0\mu_0\frac{\partial \vec{E}(t, \vec{x})}{\partial t}
\end{align}
$$

これを変形していくと次のような波動方程式が得られます。

$$
\left( \nabla^2 – \epsilon_0\mu_0\frac{\partial^2}{\partial t^2}\right) \vec{E} = 0 \\
\left( \nabla^2 – \epsilon_0\mu_0\frac{\partial^2}{\partial t^2}\right) \vec{B} = 0
$$

導出の詳細は前回の記事をご参照ください。

近軸近似

単一偏光で、かつ等位相面の法線がほぼ単一方向に揃っている場合には、電場をスカラー関数 $ u(t, x, y, z) $ で近似することができます。

さらに波が$ z $ 軸方向に進んでいると仮定して次のようにおきます。

$$
u(t,x,y,z) = \psi (x,y,z) \exp (i\omega _0 t – ik_0z) \ \ \ \ (1)
$$

ただし、波の一般的な性質として、次のような式が成り立つ。$ c $ を波の速度とした。

$$
c = \frac{\omega_0}{k_0}
$$

式(1)のように置いた理由としては、
1. 波のように時間変化しそうだから
2. 波のように位置変化しそうだから
です。

これを波動方程式に代入すると、

$$
\left( \frac{\partial ^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial ^2}{\partial z^2}\right)\psi (x,y,z) e ^{i\omega _0t – ik_0z} =
\frac{1}{c^2}\psi e^{-ik_0z} \frac{\partial ^2 e^{i\omega_0 t}}{\partial t^2}
$$

これを計算していきます。

まず右辺。

$$
\frac{1}{c^2}\psi e^{-ik_0z} \frac{\partial ^2 e^{i\omega_0
t}}{\partial t^2} = – \frac{\omega _0 ^2}{c^2} \psi e^{i\omega_0 t – ik_0z}
$$

次に左辺。

$$
\begin{align}
\frac{\partial ^2}{\partial x^2}\psi e^{i\omega _0 t – ik_0z} &=
e^{i\omega _0 t – ik_0z} \frac{\partial ^2 \psi}{\partial x^2} \\
\frac{\partial ^2}{\partial y^2}\psi e^{i\omega _0 t – ik_0z} &=
e^{i\omega _0 t – ik_0z} \frac{\partial ^2 \psi}{\partial y^2} \\
\frac{\partial ^2}{\partial z^2}\psi e^{i\omega _0 t – ik_0z} &=
e^{i\omega _0 t – ik_0z} \left( \frac{\partial ^2 \psi}{\partial z^2}
– 2k_0i\frac{\partial \psi}{\partial z} – k_0^2\psi \right)
\end{align}
$$

以上から、波動方程式は次のようになります。

$$
\left( \frac{\partial ^2}{\partial x^2} + \frac{\partial ^2}{\partial
y^2} + \frac{\partial ^2}{\partial z^2} – 2ik_0\frac{\partial}{\partial z} -k_0^2 \right) \psi (x,y,z) = -\frac{\omega_0^2}{c^2}\psi (x,y,z)
$$

ここで

$$
k_0 = \frac{\omega_0}{c}
$$

より、

$$
\left( \frac{\partial ^2}{\partial x^2} + \frac{\partial ^2}{\partial
y^2} + \frac{\partial ^2}{\partial z^2} – 2ik_0\frac{\partial}{\partial z} \right) \psi (x,y,z) = 0
$$

となります。

ここで次のような近似を使います。

$$
\left| \frac{\partial \psi}{\partial z} \right| \ll k_0\psi
$$

これはつまり、z軸方向の電場の変化が十分に穏やかであるという仮定です。

これを仮定すると、伝搬方向(z軸方向)の2階微分の項を無視でき、次のようになります。

$$
\left( \frac{\partial ^2}{\partial x^2} + \frac{\partial ^2}{\partial
y^2} – 2ik_0\frac{\partial}{\partial z} \right) \psi (x,y,z) = 0 \ \ \ \ (2)
$$

この式を近軸スカラー方程式というそうです。この方程式はちょうど 2 次元自由空間のシュレディンガー方程式の時間座標を空間座標zで置き換えた方程式と同じ形をしている点も面白いですね。

以上では Maxwell方程式から導いた電磁波の方程式に対して、近軸近似を行い、近軸スカラー方程式を導きました。
下記ではこの近軸スカラー方程式を解いていき、電場分布を求めます。

基本ガウシアンモードの導出

上記の近軸スカラー方程式をフーリエ変換を使用して解いていきます。

フーリエ変換、及びその逆変換を次のように定義します。

$$
\begin{align}
\tilde{\psi}(\nu_x, \nu _y, z) &= \int_{-\infty}^{\infty}{\rm d}x
\int_{-\infty}^{\infty}{\rm d}y \psi (x,y,z) \exp (-i\nu_xx-i\nu_yy)
\\
\psi (x,y,z) &= \frac{1}{(2\pi)^2}\int_{-\infty}^{\infty}{\rm d}\nu_x
\int_{-\infty}^{\infty}{\rm d}\nu_y \tilde{\psi}(\nu_x, \nu_y, z) \exp (i\nu_xx + i\nu_yy)
\end{align}
$$

これを使って、近軸スカラー方程式をフーリエ変換すると、例えば一項目が、次のようになります。

$$
\begin{align}
\int_{-\infty}^\infty {\rm d}x \int_{-\infty}^\infty {\rm d}y
\frac{\partial ^2 \psi}{\partial x^2} e^{-i\nu_xx -i\nu_yy}
&= \int_{-\infty}^\infty {\rm d}y e^{-i\nu_yy} \int_{-\infty}^\infty
\frac{\partial ^2 \psi}{\partial x^2} e^{-i\nu_xx} {\rm d}x \\
&= \int_{-\infty}^\infty {\rm d}y e^{-i\nu_yy} \left[\left(e^{-i\nu_xx}\frac{\partial \psi}{\partial x}\right)_{-\infty}^\infty -(-i\nu_x)\int_{-\infty}^{\infty}\frac{\partial \psi}{\partial x}e^{-i\nu_xx}{\rm d}x \right] \\
&= -\nu_x^2 \int_{-\infty}^{\infty}{\rm d}y
\int_{-\infty}^{\infty}{\rm d}x \psi e^{-i\nu_xx -i\nu_yy} \\
&= -\nu_x^2 \tilde{\psi}(\nu_x, \nu_y, z)
\end{align}
$$

なので、結局近軸スカラー方程式は次のようにフーリエ変換されます。

$$
\left( -\nu_x^2 -\nu_y^2 -2ik_0\frac{\partial}{\partial z}
\right)\tilde{\psi}(\nu_x, \nu_y, z) = 0 \ \ \ \ (3)
$$

この式は実はzに関しては簡単な微分方程式で、すぐに次のように解けます。

$$
\begin{align}
\tilde{\psi} &= A \exp \left(-\frac{\nu_x^2 + \nu_y^2}{2ik_0}z\right) \\
&= \sqrt{\frac{4\pi z_R}{k_0}} \exp \left[ -\frac{\nu_x^2 +
\nu_y^2}{2ik_0}(z-z_0+iz_R)
\right] \ \ \ \ \ (4)
\end{align}
$$

ここで積分定数Aは複素数であると考え、$ z_0, z_R $ という文字で表した。かなり意味深だが、悪いことはしていない。積分定数Aの範囲は、$ z_0, z_R $ を変化させることで、複素数の全領域にわたる。$ z_0, z_R $ はともに実数で、$ z_R $ は正でなくてはならない。

次に、$ \tilde{\psi} $ を逆フーリエ変換させ、$ \psi $ を求める。

$$
\begin{align}
\psi (x,y,z) &= \frac{1}{(2\pi)^2}\int_{-\infty}^{\infty}{\rm d}\nu_x
\int_{-\infty}^{\infty}{\rm d}\nu_y \tilde{\psi}(\nu_x, \nu_y, z) \exp
(i\nu_xx + i\nu_yy) \\
&= \frac{1}{(2\pi)^2}\int_{-\infty}^{\infty}{\rm d}\nu_x
\int_{-\infty}^{\infty}{\rm d}\nu_y \sqrt{\frac{4\pi z_R}{k_0}}\exp
\left[ -\frac{\nu_x^2 + \nu_y^2}{2ik_0}(z – z_0 + iz_R) \right] \exp
(i\nu_xx + i\nu_yy) \\
&= \frac{\sqrt{4\pi z_R}}{4\pi ^2 \sqrt{k_0}} \int_
{-\infty}^{\infty}{\rm d}\nu_y \exp (-\alpha \nu_y^2 + i\nu_yy) \int_{-\infty}^{\infty}{\rm d}\nu_x \exp (-\alpha \nu_x^2 + i\nu_xx)
\end{align}
$$

ここでは、

$$
\alpha = \frac{z-z_0+iz_R}{2ik_0}
$$

としました。

$ \nu_x $ の積分のみを行って見ると、平方完成とガウス積分を用いて、

$$
\begin{align}
\int_{-\infty}^{\infty}{\rm d}\nu_x e^{-\alpha \nu_x^2 + i\nu_xx} &=
e^{-\frac{x^2}{4\alpha}} \int_{-\infty}^{\infty}{\rm d}\nu_x
e^{-\alpha \nu_x^2} \\
&= e^{-\frac{x^2}{4\alpha}}\sqrt{\frac{\pi}{\alpha}}
\end{align}
$$

となります。

$ \nu_y $ での積分も同様にして計算して、それを合わせると次のように$ \psi $ が計算できます。

$$
\begin{align}
\psi ^{{\rm FG}} (x,y,z) &= \frac{\sqrt{4\pi z_R}}{4\pi^2
\sqrt{k_0}}\frac{\pi}{\alpha} e^{-\frac{1}{4\alpha}(x^2+y^2)} \\
&= \sqrt{\frac{k_0}{\pi z_R}} \frac{iz_R}{z-z_0+iz_R}\exp
\left(-\frac{ik_0}{z-z_0+iz_R}\frac{x^2 + y^2}{2}\right) \ \ \ \ (5)
\end{align}
$$

この解は基本ガウシアンモード(Fundamental Gaussian mode)とよばれ、近軸スカラー波動方程式の解となります。

ここまで、近軸スカラー関数式(2)の解を求めようと、式をフーリエ変換して式(3)に変形し、この解として式(4)を得ましたが、この式(4)は実は$ \mu_x, \mu_y $ を含むが$ z $ を含まない関数を掛けても式(3)をみたします。

つまり

$$
\tilde{\psi} = \sqrt{\frac{4\pi z_R}{k_0}} \exp \left[ -\frac{\nu_x^2 +
\nu_y^2}{2ik_0}(z-z_0+iz_R)
\right]f(\nu_x, \nu_y)
$$

これによっていくらでも複雑なbeam解を得ることができます。このようにして得られる解が高次モードとなります。

実際に高次モードの(レーザー)光をどのようにつくかは置いておいて、存在は許されることがわかった。
もちろんこれらはMaxwell方程式(と近軸近似)から導かれる結論です。

次に、以上の結果をさらに式変形してみて、その性質をみていきます。

基本ガウシアンモードの特徴

上記でもとまった近軸近似における電場の解(式(5))をもう一回書いてみる。

$$
\psi ^{{\rm FG}} (x,y,z) = \sqrt{\frac{k_0}{\pi z_R}} \frac{iz_R}{z-z_0+iz_R}\exp
\left(-\frac{ik_0}{z-z_0+iz_R}\frac{x^2 + y^2}{2}\right) \ \ \ \ (5)
$$

ただし、$ \nu_x, \nu_y $ の関数は掛けずに基本モードとした。

この式に対して、次の変数変換を行ってみます。

$$
\begin{align}
q(z) &= z + iz_R \label{eq:qnosiki} \\
z_R &= \frac{k_0w_0^2}{2}
\end{align}
$$

ここで$ k_0, \omega_0 $ は平面波$ \exp(i\omega_0t – ik_0z) $ のものと同じで、定数$ z_0 $ は0としました。

すると、式(5)は次のように書き換えられます。

$$
\psi ^{{\rm FG}} (x,y,z)
= \sqrt{\frac{k_0}{\pi z_R}} \frac{iz_R}{q(z)}\exp
\left(-\frac{ik_0}{q(z)}\frac{x^2 + y^2}{2}\right) \ \ \ \ (6)
$$

ここでは積分定数に値を代入しただけで、特殊な条件を課したわけでは無いです。

さらに長さの次元を持つ実数量$ R(z), S(z) $ を次のように定義しておきます。

$$
\frac{1}{q(z)} = \frac{1}{R(z)} + \frac{1}{iS(z)}
$$

これもいわば、変数を書き換えただけですね。

さらにさらに次のように$ \omega(z) $ を定義します。

$$
\frac{k_0}{2S(z)}=\frac{1}{\omega^2(z)}
$$

このように定義した$ R(z), \omega(z) $ を使って、式(6)を表すようにしてみる。そうすると式(6)のかなり物理的な意味が見えてきます。

まず$ q(z) $ の式と$ \frac{1}{q(z)} $ の式から$ R(z), S(z) $ を求めます。これは簡単に次のようになります。

$$
\begin{align}
R(z) &= z\left[ 1+\left( \frac{z_R}{z} \right)^2 \right] \\
S(z) &= z_R\left[ 1+\left( \frac{z}{z_R} \right)^2 \right]
\end{align}
$$

ここで$ S(z) $ が表せたので、$ \omega(z) $ も次のように書き下すことができます。

$$
\begin{align}
\omega &= \sqrt{\frac{2}{k_0}z_R\left[
1+\left(\frac{z}{z_R}
\right)^2 \right]} \\
&= \omega_0\sqrt{1+\left( \frac{z}{z_R} \right)^2}
\end{align}
$$

これらを用いて、式(6)つまり$ \psi^{FG} $ を変形していきます。

$$
\begin{align}
\psi^{FG}(x,y,z) &= \sqrt{\frac{k_0}{\pi
z_R}}\frac{iz_R}{q(z)}\exp\left[ -\frac{ik_0}{q(z)}\frac{x^2+y^2}{2}
\right] \\
&= \sqrt{\frac{k_0}{\pi z_R}}
\left[
\frac{iz_R}{z\left[ 1+\left( \frac{z_R}{z} \right)^2 \right]}
+ \frac{iz_R}{iz_R\left[ 1+\left( \frac{z}{z_R} \right)^2 \right]}
\right]
\exp\left[
-\left( \frac{ik_0}{R} + \frac{k_0}{S} \right)
\frac{x^2+y^2}{2}
\right] \\
&= \sqrt{\frac{k_0}{\pi z_R}}
\left[
\frac{izz_R}{z^2+z_R^2}
+ \frac{z_R^2}{z^2+z_R^2}
\right]
\exp\left[
-\left( \frac{ik_0}{2R} + \frac{k_0}{2S} \right)
(x^2+y^2)
\right] \\
&= \sqrt{\frac{2}{\pi}}
\sqrt{\frac{k_0z_R}{2}}
\frac{1}{z^2 + z_R^2}
(iz + z_R)
\exp\left[
-\left( \frac{ik_0}{2R} + \frac{1}{\omega^2} \right)
(x^2+y^2)
\right] \\
&= \sqrt{\frac{2}{\pi}}
\frac{k_0\omega_0}{2}
\frac{1}{z^2+z_R^2}
\sqrt{z^2+z_R^2}
e^{i\arctan\left( \frac{z_R}{z} \right)}
\exp\left[
-\left( \frac{ik_0}{2R} + \frac{1}{\omega^2} \right)
(x^2+y^2)
\right] \\
&= \sqrt{\frac{2}{\pi}}
\frac{k_0\omega_0^2}{2}
\frac{1}{\omega_0\sqrt{z^2+z_R^2}}
\exp\left[
-\left( \frac{ik_0}{2R} + \frac{1}{\omega^2} \right)
(x^2+y^2)
+ i\arctan\left( \frac{z_R}{z} \right)
\right] \\
&= \sqrt{\frac{2}{\pi}}
\frac{1}{\omega (z)}
\exp\left[
-\left( \frac{1}{\omega^2(z)}+\frac{ik_0}{2R(z)} \right)
(x^2+y^2)
+i\zeta(z)
\right] \\
&= \sqrt{\frac{2}{\pi}}
\frac{1}{\omega (z)}
\exp \left(-\frac{r^2}{\omega^2(z)}\right)
\exp\left[
-i\left(
\frac{k_0r^2}{2R(z)} – \zeta(z)
\right)
\right] \ \ \ \ (7)
\end{align}
$$

ただし、ここでは、

$$
\begin{align}
\zeta(z) &= \arctan\left( \frac{z}{z_R} \right) \\
r &= \sqrt{x^2+y^2}
\end{align}
$$

としました。

この$ \zeta(z) $ がGouy phaseと呼ばれる量となります。

ここで原点に戻って電場$ u(t,x,y,z) $ を求めると次のようになります。

$$
\begin{align}
u(t,x,y,z) &= \psi (x,y,z) \exp (i\omega _0 t – ik_0z) \ \ \ \ (1) \\
&= \sqrt{\frac{2}{\pi}}
\frac{1}{\omega (z)}
\exp\left( -\frac{r^2}{\omega^2(z)}\right)
\exp\left[
-i\left(
\frac{k_0r^2}{2R(z)} – \zeta(z)
\right)
\right] \exp (i\omega _0 t – ik_0z) \\
&= \sqrt{\frac{2}{\pi}}
\frac{1}{\omega (z)}
\exp\left(-\frac{r^2}{\omega^2(z)}\right)
\exp\left[
-i\left(
k_0z +\frac{k_0r^2}{2R(z)} – \zeta(z)
\right)
\right] \exp (i\omega _0 t) \ \ \ \ (8)
\end{align}
$$

だいぶ形が見えてきました。

ここで、$ z=0, r=0 $ での電場$ u $ を計算して見ると上記の式の形がさらに見えやすくなります。

まずは$ z=0 $ という条件によって次のようになります。

$$
\begin{align}
u(t, r, z=0) &= \sqrt{\frac{2}{\pi}}\frac{1}{\omega(0)}\exp\left( -\frac{r^2}{\omega^2(0)} \right)\exp\left( -\frac{ik_0r^2}{2R(0)} \right)\exp(i\omega_0t) \\
&= \sqrt{\frac{2}{\pi}}\frac{1}{\omega_0}\exp\left( -\frac{r^2}{\omega^2_0} \right)\exp(i\omega_0t)
\end{align}
$$

ここでは、$ \omega(z=0)=\omega_0, R(z=0)\rightarrow \infty, \zeta(z=0)=0 $ を利用しました。

次に$ r = 0 $ という条件で次のようになります。

$$
u(t, r=0, z=0) = \sqrt{\frac{2}{\pi}}\frac{1}{\omega_0}\exp(i\omega_0t) \equiv E_0(t)
$$

これはつまりz軸の原点(z=0)でかつ、中心(r=0)における電場なので$ E_0(t) $ という記号を付けました。
電場は波長に応じた固有振動数で時間的には揺れ動いているので時間tの依存性は残っています。

これを式(8)に代入すると次のようにはっきりと物理的な意味が見えてきます。

$$
u(t, r, z) = E_0(t)\frac{\omega_0}{\omega(z)}\exp\left( -\frac{r^2}{\omega^2(z)} \right)\exp\left[ -i\left( k_0z+\frac{k_0r^2}{2R(z)} – \zeta(z) \right) \right] \ \ \ \ (9)
$$

電場分布式に関する考察

上記で求められた電場の分布式である式(9)に関して、観察していく。

式を再記しておきます。

$$
u(t, r, z) = E_0(t)\frac{\omega_0}{\omega(z)}\exp\left( -\frac{r^2}{\omega^2(z)} \right)\exp\left[ -i\left( k_0z+\frac{k_0r^2}{2R(z)} – \zeta \right) \right] \ \ \ \ (9)
$$

ひとつずつ項(?)をみていきます。

$ E_0(t) $

この項、z=0, r=0の時に、唯一残る項(1でない項)です。

つまり原点での電場の時間変化を表しています。

$ \exp\left( -\frac{r^2}{\omega^2(z)} \right) $

この項は飛行光線軸から離れた位置での電場が$ \omega(z) $ を基準にして減少していくことを意味しています。

これがまさにGaussian beamを特徴付けている項ですね。

となると、項内の$ \omega(z) $ というのがいわばレーザー光の太さを意味していることになることがわかります。

つまり、$ r $ が$ \omega(z) $ に一致した位置では、電場振幅が$ \frac{1}{e} $ になるということを意味しています。

$ \frac{\omega_0}{\omega(z)} $

上記で、$ \omega(z) $ がレイザー光線の幅であるとわかったので、この項が意味することは、

レーザー光線の幅が太くなる部分では、太さに反比例して電場が弱くなる、

ということです。

$ \exp\left[ -i\left( k_0z+\frac{k_0r^2}{2R(z)} – \zeta(z) \right) \right] $

この項は$ \exp(i\theta) $ の形をしていますので、電場の空間的振動を表しています。

振動はsin波で表されると考えると、同時にこの項はsin波の位相を表しているということになります。

$ \exp(-ikz) $ 部分はz位置に従って波が変化することを表しています。

$ \exp\left( -i\frac{k_0r2}{2R(z)} \right) $ はある距離$ R(z) $ に特徴付けられて、光軸からの距離に応じて電場が振動していることを意味しています。

最後に$ \exp(-i\zeta(z)) $ はGouy phase shiftと呼ばれる効果を生み出している項でこの後に説明します。

Gouy phase $ \zeta(z) $ の物理的な意味

上記までで、$ \zeta(z) $ を除く項の物理的な理解ができました。

ここではこのGouy phase $ \zeta(z) $ の物理的な意味を考えていきます。

実はかなり簡単です。

平面波の電場分布$ E_p $について時間の成分を分けると、次のように表せます。

$$
E_p(t, z) = E_{p0}(t)\exp(-ik_0z)
$$

空間位相変化部分だけを比較してみます。

平面波	$ \exp(-ik_0z) $
Gaussian beam	$ \exp\left[ -i\left( k_0z+\frac{k_0r^2}{2R(z)} – \zeta(z) \right) \right] $

Gaussian beamの第一項目は平面波と一致していますね。

第二項目はちょっと複雑に見えるけど、レーザー光線の近傍に興味があるとすると、$ r $ は小さく、$ R(z) $ はそれと比べて大きな数値となりますので、この第二項は0に近似できます。

となると最後に第三項の意味が見えてきます。

そう、第三項のGouy phase項はGaussian beamと平面波との位相のズレを意味しているのです。

$$
\zeta(z) = \arctan\left( \frac{z}{z_R} \right)
$$

ですので、zが進むにつれて、位相のズレが生じてきます。

符号がマイナスなので、Gaussian beamの方が位相が遅れているということもわかります。
平面波とGaussian beamを競争させて発射した時、光面・光線が進むほどにGaussian beamの電場位相が遅れてくるということですね。

イメージ図を作ってみました。

縦軸が電場の強さで横軸が位置(z)です。
紫線が平面波の変化で、緑線がGaussian beamを表しています。$ z_R $ は適当に10としています。

この図から、ある適当な位置に注目すると、平面波を表す紫線の方が先に進んでいる、位相が進んでいることがわかります。
またzが進むにつれてその位相差が大きくなっていることが見えます。

まとめ

今回は、真空中のMaxwell方程式から得られる波動方程式

$$
\left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2} \right)\vec{E} = \frac{1}{c^2}\frac{\partial^2 \vec{E}}{\partial t^2}
$$

を近軸近似(飛行方向の電場は穏やかにしか変化しないという条件を課した)

$$
\left| \frac{\partial \psi}{\partial z} \right| \ll k_0\psi
$$

して、フーリエ変換を用いて解いた。

計算された電場分布が下記でした。

$$
u(t, r, z) = E_0(t)\frac{\omega_0}{\omega(z)}\exp\left( -\frac{r^2}{\omega^2(z)} \right)\exp\left[ -i\left( k_0z+\frac{k_0r^2}{2R(z)} – \zeta(z) \right) \right] \ \ \ \ (9)
$$

この式中に現れた$ \zeta(z) $ はGouy phaseと呼ばれ、平面波に対するGaussian beamの位相遅れを意味することがわかりました。