1.   引　言

高能超短激光脉冲在高能粒子加速、强场物理、惯性约束聚变“快点火”等领域具有广泛的应用^[1-3]. 光参量啁啾脉冲放大(OPCPA)技术在实现高功率、高能量、超短脉冲方面有着很大潜力^[4-6], 研发超短脉冲激光系统数值模拟仿真平台, 对于优化激光装置设计、提升激光质量十分重要. 其中, 仿真平台中对超短脉冲聚焦后的光学特性进行模拟计算时, 通常采用菲涅尔远场衍射传输方程^[7-9]. 首先需要计算超短脉冲中不同波长分量的光场传输分布状态, 然后利用傅里叶逆变换得到脉冲时空分布特性. 但是对于不同的波长分量, 远场算法得到的各自聚焦光场网格大小不相等, 导致后续傅里叶逆变换无法直接进行, 需要进行重新采样, 增大了计算量; 此外由于在脉冲展宽压缩、参量放大环节光场计算范围受到限制, 还存在聚焦光场分辨率差的问题.

有关聚焦光场的衍射传输计算问题, 较为经典的是Talanov 算法^[10], 它能根据传输距离自动变换计算网格以适应聚焦光束实际范围, 但不适用于聚焦场瑞利区; Feigenbaum等^[11]对Talanov 算法进行了修正, 通过对输入场再采样, 微调相位弯曲因子, 避免了Talanov 算法不适用于瑞利区的缺陷, 但再采样使算法变得复杂. Kozacki等 ^[12]在计算高数值孔径的衍射场时, 通过多次计算网格点数较小的傅里叶变换, 使光场分辨率得到提升; 杨美慧等^[13]采用两步坐标扩展变换方法, 按比例缩小输出窗口范围, 计算焦点处光场分布; Hu等^[14]利用啁啾Z变换调整空间频率计算范围, 进而得到细致的聚焦光场分布. 文献[15]中的两步衍射法假设除了输入场和输出场以外存在一个虚拟衍射场, 利用虚拟衍射场和输入、输出场之间的对应关系, 通过改变虚拟场的位置来调整输出场计算范围.

不同于前述文献方法, 本文直接从菲涅尔远场衍射公式出发, 在离散计算过程中, 对输出场独立采样, 以满足输出场计算范围灵活可调的要求, 在此基础上推导给出了具备快速傅里叶算法的衍射传输数学表达式. 然后进一步分析讨论了使用该算法时需要满足的采样条件. 最后利用该算法对超短脉冲聚焦特性进行模拟仿真计算. 由于该算法在提高聚焦光场空间分辨率的同时, 保持了各波长分量计算网格的一致性, 计算效率大大提高. 在该算法的支持下, 模拟发现超短脉冲聚焦场暗环区域表现出强烈的时空耦合特征.

2.   理论分析

2.1   输出场计算范围可调的理论公式推导

设输入和输出光场复振幅分别为${U_{{\text{in}}}}(x, y)$和${U_{{\text{out}}}}(x', y')$, 衍射传输距离为$z$. 对聚焦过程进行模拟计算, 通常采用菲涅尔衍射方程的远场计算形式^[9], 输出场的复振幅表达式为

$$\begin{split} & {U_{{\text{out}}}}\left( {x',y'} \right) = \frac{1}{{{\text{i}}\lambda z}} \exp\Big\{ {\text{j}} \frac{k}{{2z}} \left[(x')^2 +(y')^2\right] \Big\}\\ & ~~~~\times \iint_{-\infty }^\infty \Big\{ U_{\rm in} (x,y) \exp\Big[ {\text{j}} \frac{k}{{2z}} (x^2 + y^2) \Big] \Big\} \\ &~~~~\times \exp \Big[ - {\text{j}} \dfrac{{2{\text{π }}}}{{\lambda z}}( {xx' + yy'} )\Big] {\text{d}}x{\text{d}}y, \end{split}$$

(1)

式中, k表示波矢($k = {{2{\text{π}}}/\lambda}$), $\lambda$是光波长, 其中的积分是输入场复振幅与一个二次相位因子乘积的傅里叶变换, 空间频率和坐标的映射关系为${f_x} = {x'}/({{\lambda z}})$, $f_y = y'/(\lambda z)$. 显然, 输出场计算范围依赖于波长, 因此利用(1)式数值计算时, 不同的波长分量对应聚焦光场网格间距不相等, 使得后续直接进行傅里叶逆变换变得困难; 另一方面, 数值仿真平台中, 为了兼顾参量放大过程, 输入场计算范围通常只需包含脉冲光束有效能量范围即可, 此时可能会使聚焦光场分辨率过低.

为了解决上述问题, 有必要发展输出光场计算网格能灵活调整的传输算法, 既能使不同波长分量输出光场的计算网格大小保持一致, 也能根据分辨率需求合理选择输出场范围大小. 在离散数值计算过程中, 设输入和输出光场计算范围分别为${L_x} \times {L_y}$和${L'_x} \times {L'_y}$, 抽样点数均为$N$, 则采样间隔分别是$\Delta x{\text{ = }}{{{L_x}} \mathord{\left/ {\vphantom {{{L_x}} N}} \right. } N}$, $\Delta y{\text{ = }}{{{L_y}} \mathord{\left/ {\vphantom {{{L_y}} N}} \right. } N}$, $\Delta x'{\text{ = }}{{{L'_x}} /N}$, $\Delta y'{\text{ = }} {{{L'_y}} / N}$, 离散采样点坐标可表示为

$$\left\{ \begin{aligned} &x' = {x'_0} + p\Delta x' ,\;&p = 0,1,2, \cdot \cdot \cdot ,N - 1, \\[2mm] &y' = {y'_0} + q\Delta y' ,\;&q = 0,1,2, \cdot \cdot \cdot ,N - 1, \\[2mm] &x = {x_0} + m\Delta x, \;&m =0,1,2, \cdot \cdot \cdot ,N - 1, \\[2mm] &y = {y_0} + n\Delta y ,\;& n = 0,1,2, \cdot \cdot \cdot ,N - 1, \end{aligned} \right.$$

(2)

其中$\left( {{x_0}, {y_0}} \right)$, $\left( {{x'_0}, {y'_0}} \right)$分别表示输入场和输出场的起始坐标. 用数值求和替代积分后, (1)式的离散形式为

$$\begin{split} & U_\text{out} \big( x'_0 + p\Delta x',y_0' + q\Delta y'\big) \\ =\;& \frac{{\Delta x\Delta y}} {{{\text{i}}\lambda z}} \text{e}^{\text{j} \tfrac{k}{2z} \left[ (x'_0 + p\Delta x')^2 + (y'_0 + q\Delta y' )^2 \right]} \sum\limits_{n = 0}^{N - 1} \sum\limits_{m = 0}^{N - 1} \Big\{ U_\text{in} \left( x_0 + m\Delta x ,y_0 + n\Delta y \right) \\ &\times \text{e}^{\text{j} \tfrac{k}{2z}\left[ {{{\left( {{x_0} + m\Delta x } \right)}^2} + {{\left( {{y_0} + n\Delta y} \right)}^2}} \right]} {\rm e}^{-{\text{j}} \tfrac{{2{\text{π }}}} {\lambda z}\left[ {\left( {{x_0} + m\Delta x} \right)\left( {{x'_0} + p\Delta x'} \right) + \left( {{y_0} + n\Delta y} \right)\left( {{y'_0} + q\Delta y' } \right)} \right]} \Big\}. \end{split}$$

(3)

此处, 不能直接对(3)式进行快速离散傅里叶变换计算, 因为它虽然形式上和(1)式类似, 但是输出场独立取样以后, 其采样间隔不一定满足傅里叶变换要求的空间频率和坐标映射关系. (3)式中${{\text{e}}^{ - {\text{j}}\tfrac{{2{\text{π }}}}{{\lambda z}}\left[ {\left( {{x_0} + m\Delta x} \right)\left( {{x'_0} + p\Delta x'} \right){\kern 1 pt} + \left( {{y_0} + n\Delta y} \right)\left( {{y'_0} + q\Delta y'{\kern 1 pt} } \right)} \right]}}$项展开可以得到

$$\begin{split} &{U_{{\text{out}}}} \big( {{x'_0} + p\Delta x',{y_0'} + q\Delta y' } \big) \\ =\;& \frac{{\Delta x\Delta y}}{{{\text{i}}\lambda z}} {{\text{e}}^{{\text{j}}\tfrac{k}{{2z}}\left[ {{{\left( {{x'_0} + p\Delta x'} \right)}^2} + {{\left( {{y'_0} + q\Delta y' } \right)}^2}} \right]}}{{\text{e}}^{ - {\text{j}}\tfrac{{2{\text{π }}}}{{\lambda z}}\left[ {{x_0}p\Delta x' + {y_0}q\Delta y'} \right]}} \sum\limits_{n = 0}^{N - 1} \sum\limits_{m = 0}^{N - 1} \Big\{ {U_{{\text{in}}}}\left( {{x_0} + m\Delta x ,{y_0} + n\Delta y} \right) \\ &\times {{\text{e}}^{{\text{j}}\tfrac{k}{{2z}}\left[ {{{\left( {{x_0} + m\Delta x } \right)}^2} + {{\left( {{y_0} + n\Delta y} \right)}^2} + \left( {{x_0} + m\Delta x} \right){x'_0} + \left( {{y_0} + n\Delta y} \right){y'_0}} \right]}} {{{\text{e}}^{ - {\text{j}}\tfrac{{2{\text{π }}}}{{\lambda z}}\left[ {mp\Delta x\Delta x' + nq\Delta y\Delta y'} \right]}}} \Big\}, \end{split}$$

(4)

为使式(4)求和部分能够进行快速卷积计算, 对其中的$p$和$q$变量进行以下变换处理:

$$p=\frac{{m}^{2}+{p}^{2}-{(m-p)}^{2}}{2m},~~~~ q=\frac{{n}^{2}+{q}^{2}-{(n-q)}^{2}}{2n}.$$

(4)式成为

$$\begin{split} &{U_{{\text{out}}}}\left( {{x'_0} + p\Delta x',{y_0'} + q\Delta y' } \right) \\ =\;& \frac{{\Delta x\Delta y}}{{{\text{i}}\lambda z}} {{\text{e}}^{{\text{j}}\tfrac{k}{{2z}}\left[ {{{\left( {{x'_0} + p\Delta x'} \right)}^2} - 2{x_0}p\Delta x' - {p^2}\Delta x\Delta x' + {{\left( {{y'_0} + q\Delta y' } \right)}^2} - 2{y_0}q\Delta y' - {q^2}\Delta y\Delta y'} \right]}} \sum\limits_{n = 0}^{N - 1} \sum\limits_{m = 0}^{N - 1} \Big\{ {U_{{\text{in}}}} \left( {{x_0} + m\Delta x ,{y_0} + n\Delta y} \right) \\ &\times {{\text{e}}^{{\text{j}}\tfrac{k}{{2z}}\left[ {{{\left( {{x_0} + m\Delta x } \right)}^2} - 2\left( {{x_0} + m\Delta x} \right)x' - {m^2}\Delta x\Delta x' + {\left( {{y_0} + n\Delta y} \right)}^2_0 - 2\left( {{y_0} + n\Delta y} \right){y'_0} - {n^2}\Delta y\Delta y'} \right]}} \\ & \times {\rm{e}}^{ {{\text{j}}\tfrac{{\text{π }}}{{\lambda z}}\left[ {{{\left( {m - p} \right)}^2}\Delta x\Delta x' + {{\left( {n - q} \right)}^2}\Delta y\Delta y'} \right]}} \Big\},\\[-15pt] \end{split}$$

(5)

进一步运算后得到

$$\begin{split} &{U_{{\text{out}}}}\left( {{x'_0} + p\Delta x',{y_0'} + q\Delta y' } \right) \\ =\;& \dfrac{{\Delta x\Delta y}}{{{\text{i}}\lambda z}}{{\text{e}}^{{\text{j}}\tfrac{k}{{2z}}\left[ {\left( {1 - \tfrac{{\Delta x}}{{\Delta x'}}} \right){{\left( {p\Delta x' + \tfrac{{\Delta x'{x'_0} - \Delta x'{x_0}}}{{\Delta x' - \Delta x}}} \right)}^2} + \left( {1 - \tfrac{{\Delta y}}{{\Delta y'}}} \right){{\left( {q\Delta y' + \tfrac{{\Delta y'{y'_0} - \Delta y'{y_0}}}{{\Delta y' - \Delta y}} } \right)}^2}} \right]}} \\ &\times \sum\limits_{n = 0}^{N - 1} {\sum\limits_{m = 0}^{N - 1} {\bigg\{ {{U_{{\text{in}}}}\left( {{x_0} + m\Delta x ,{y_0} + n\Delta y} \right){\text{e}}} } } { ^{{\text{j}}\tfrac{k}{{2z}}\left[ {\left( {1 - \tfrac{{\Delta x'}}{{\Delta x}}} \right){{\left( {m\Delta x - \tfrac{{\Delta x{x'_0} - \Delta x{x_0}}}{{\Delta x - \Delta x'}} } \right)}^2} + \left( {1 - \tfrac{{\Delta y'}}{{\Delta y}}} \right){{\left( {n\Delta y - \tfrac{{\Delta y{y'_0} - \Delta y{y_0}}}{{\Delta y - \Delta y'}}} \right)}^2}} \right]}} \\ & \times {\rm{e}}^{\text{j}\tfrac{k}{{2z}} \left[ {{{\left( {m - p} \right)}^2}\Delta x\Delta x' + {{\left( {n - q} \right)}^2}\Delta y\Delta y'} \right]} \bigg\}, \\[-15pt] \end{split}$$

(6)

为使表达式更加简洁, 引入以下变量代换:

$$\begin{gathered} u = 1 - \dfrac{{\Delta x'}}{{\Delta x}}, \; {u'} = 1{{ - }}\dfrac{{\Delta x }}{{\Delta x'}} ~~ (\Delta x' \ne \Delta x ) , \qquad v = 1 - \dfrac{{\Delta y'}}{{\Delta y}}, \; {v'} = 1{{ - }}\dfrac{{\Delta y }}{{\Delta y'}} ~~ (\Delta y' \ne \Delta y), \end{gathered}$$

化简后得到

$$\begin{split} &{U_{{\text{out}}}}\left( {{x'_0}{\text{ + }}p\Delta x',{y'_0} + q\Delta y'} \right) \\ =\;& \dfrac{{\Delta x\Delta y}}{{{\text{i}}\lambda z}}\exp \left\{ {\dfrac{{{\text{i\pi }}}}{{\lambda z}}\left[ {u'{{\left( {p\Delta x' + \dfrac{{{x'_0} - {x_0}}}{{u'}}} \right)}^2} + {v'}{{\left( {q\Delta y' + \dfrac{{{y'_0} - {y_0}}}{{v'}}} \right)}^2}} \right]} \right\}\\ &\times \sum\limits_{n = 0}^{N - 1} {\sum\limits_{m = 0}^{N - 1} {\left\{ {{U_{{\text{in}}}}\left( {{x_0} + m\Delta x,{y_0} + n\Delta y} \right)} \right.} } \exp \left\{ {\left. {\dfrac{{{\text{i\pi }}}}{{\lambda z}}\left[ {u{{\left( {m\Delta x - \dfrac{{{x'_0} - {x_0}}}{u}} \right)}^2} + v{{\left( {n\Delta y - \dfrac{{{y'_0} - {y_0}}}{v}} \right)}^2}} \right]} \right\}} \right. \\ &\left. { \times \exp \left[ {\dfrac{{{\text{i\pi }}\Delta x'}}{{\lambda z\Delta x}}{{\left( {p\Delta x - m\Delta x} \right)}^2} + \dfrac{{{\text{i\pi }}\Delta y'}}{{\lambda z\Delta y}}{{\left( {q\Delta y - n\Delta y} \right)}^2}} \right]} \right\} ,\end{split}$$

(7)

(7)式中的求和运算表示输入场${U_{{\text{in}}}}\left( {{x_0} + m\Delta x, {y_0} + n\Delta y} \right)$乘上一相位因子后, 再和序列$\exp \bigg[ \dfrac{{{\text{i\pi }}\Delta x'}}{{\lambda z\Delta x}}{{\left( {m\Delta x} \right)}^2} + \dfrac{{{\text{i\pi }}\Delta y'}}{{\lambda z\Delta y}}{{\left( {n\Delta y} \right)}^2} \bigg]$进行卷积, 由于二者参与卷积计算的步长均为$\Delta x$($y$方向同理), 已经可以使用快速傅里叶算法进行卷积运算, 其中卷积核的傅里叶变换解析表达式为

$$\begin{gathered} \Im \left\{ {\exp \left[ {\dfrac{{{\text{i\pi }}\Delta x'}}{{\lambda z\Delta x}}{{\left( {m\Delta x} \right)}^2} + \dfrac{{{\text{i\pi }}\Delta y'}}{{\lambda z\Delta y}}{{\left( {n\Delta y} \right)}^2}} \right]} \right\}{\text{ = }} {\text{i}}\lambda z\sqrt {\dfrac{{\text{1}}}{{\Delta x\Delta x'}}} \sqrt {\dfrac{{\text{1}}}{{\Delta y\Delta y'}}} \exp \left[ { - {\text{i\pi }}\lambda z\left( {\dfrac{{\Delta x}}{{\Delta x'}}f_x^2 + \dfrac{{\Delta y}}{{\Delta y'}}f_y^2} \right)} \right], \end{gathered}$$

(8)

符号$\Im$表示傅里叶变换, 利用(8)式将(9)式进一步改写为以下直观的快速计算形式:

$$\begin{split} & {U}_{\text{out}}\left(x'_{0}+p\Delta x',y'_{0}+q\Delta y'\right) \\ =\;& \sqrt{\dfrac{\Delta x}{\Delta x'}}\sqrt{\dfrac{\Delta y}{\Delta y'}}\cdot {\rm{exp}}\left\{\dfrac{\text{iπ}}{\lambda z}\left[{u}'{\left(p\Delta x'+\dfrac{x'_{0}-x_{0}}{{u}'}\right)}^{2}+{v}{'}{\left(q\Delta y'+\dfrac{y'_{0}-y_{0}}{{v}'}\right)}^{2}\right]\right\}\\ &\times {\Im }^{-1}\left\{\Im \left\{{U}_{\text{in}}\left(x_{0}+m\Delta x,y_{0}+n\Delta y\right){\rm{exp}}\left[u\dfrac{\text{iπ}}{\lambda z}{\left(m\Delta x-\dfrac{x'_{0}-x_{0}}{u}\right)}^{2}+v\dfrac{\text{iπ}}{\lambda z}{\left(n\Delta y-\dfrac{y'_{0}-y_{0}}{v}\right)}^{2}\right]\right\}\right.\\ & \left. \times {\rm{exp}}\left[-\text{iπ}\lambda z\left(\dfrac{\Delta x}{\Delta x'}{f}_x^{2}+\dfrac{\Delta y}{\Delta y'}{f}_y^{2}\right)\right]\right\},\\[-15pt] \end{split}$$

(9)

式中, 符号${\Im ^{ - 1}}$表示逆傅里叶变换. 当输入场和输出场中心重合且计算范围相同时, 式(9)简化为

$$\begin{split} &{U_{{\text{out}}}}({x'_0} + p\Delta x',{y'_0} + q\Delta y') \\ =\;& \dfrac{{1 }}{{{\text{i}}\lambda z}} \cdot {\Im ^{ - 1}}\big\{ \Im \left[ {{U_{{\text{in}}}}({x_0} + m\Delta x,{y_0} + n\Delta y)} \right] \\ & \times \exp \left[ { - {\text{i\pi }}\lambda z\left( {f_x^2 + f_y^2} \right)} \right] \big\}, \end{split}$$

(10)

该式回归到文献[8, 9]中菲涅尔衍射算法的近场计算形式.

2.2   聚焦场采样条件讨论

在OPCPA飞秒激光装置中, 压缩后的脉冲通常利用离轴抛面镜进行聚焦^[4]. 输入光场通过焦距为$f$的离轴抛面镜后, 其离散化复振幅可表示为

$$\begin{split} &~~~{U_{{\text{in}}}}\left( {{x_0} + m\Delta x,{y_0} + n\Delta y} \right)\\ &\times \exp\Big\{ { - {\text{j}}\dfrac{{\text{π}}}{{\lambda f}}\big[ {{{( {{x_0} + m\Delta x} )}^2} + {{( {{y_0} + n\Delta y} )}^2}} \big]} \Big\}, \end{split}$$

(11)

将其代入(9)式后, (9) 式中傅里叶变换项的指数部分变为

$$\exp \left[ {{\text{i\pi }}\left( {\dfrac{u}{{\lambda z}} - \dfrac{1}{{\lambda f}}} \right){{\left( {{x_0} + m\Delta x} \right)}^2}} \right].$$

(12)

为简化分析, 令输入场和输出场中心重合, 根据奈奎斯特采样定理, 任意两个相邻采样值之间的相位差应不大于π, 以$x$方向为例(以下均同), 存在以下不等式:

$${\left| {{\text{π}}\Big( {\dfrac{u}{{\lambda z}} - \dfrac{1}{{\lambda f}}} \Big) ( { - N\Delta x + 2m\Delta x} )\Delta x} \right|_{m = N}} \leqslant {\text{π}},$$

(13)

其中, $m = N$对应输入场边缘位置所在点, 代入$u$的表达式展开得到

$$\Delta x - \dfrac{{z\Delta x}}{f} - \dfrac{{\lambda z}}{{N\Delta x}} \leqslant \Delta x' \leqslant \Delta x - \dfrac{{z\Delta x}}{f}{\text{ + }}\dfrac{{\lambda z}}{{N\Delta x}}.$$

(14)

通常, 输入场计算范围大于光束能量范围, $m$的取值可以小于$N$, 所以(14)式属于充分条件. 同理, (9)式中与傅里叶变换项相乘的另一指数项也需要满足采样定理, 则应有以下不等式:

$$\left| {2{\text{π }}\lambda z\left( {\dfrac{{\Delta x}}{{\Delta x'}}{f_x} \cdot \dfrac{1}{{N\Delta x}}} \right)} \right| \leqslant {\text{π }}{\text{. }}$$

(15)

理论上, 上述不等式中${f_x}$可以取到离散计算时的最高频率${1 \mathord{\left/ {\vphantom {1 {\left( {2\Delta x} \right)}}} \right. } {\left( {2\Delta x} \right)}}$处, 但实际上输入光场空间频率带宽范围是有限的, 输入光场和相位变换(12)式乘积共同决定了空间频率带宽, 即(9)式中傅里叶变换项对应的空间频率带宽存在上限${B_{\max }}$, 则(15)式中${f_x}$的最大取值应不大于${B_{\max }}$, 由此得到

$$\Delta x' \geqslant {{2\lambda z{B_{\max }}}}/{N} ,$$

(16)

在离散计算时, 应保证${1}/({{2\Delta x}} )\geqslant {B_{\max }}$; 并且当衍射传播距离$z$和焦距$f$大致相等时, (14)式中最左边项为负数, 结合(14)式和(16)式可得到输出场采样约束条件:

$$\dfrac{{2\lambda z{B_{\max }}}}{N} \leqslant \Delta x' \leqslant \Delta x{\text{ + }}\dfrac{{\lambda z}}{{N\Delta x}} - \dfrac{{z\Delta x}}{f} ,$$

(17)

实际计算时, 通常使$\Delta x'$接近 ${{2\lambda z{B_{\max }}}}/{N}$, 尽可能提高输出光场分辨率.

3.   超短脉冲聚焦模拟计算

输入超短光脉冲的复振幅由下式给出, 其空间为$n$阶圆形超高斯分布, 时间为高斯分布:

$$\begin{split} & {E_{{\text{in}}}}(x,y,t) = \exp \bigg[ { - {{\left( {\dfrac{{{x^2} + {y^2}}}{{{R^2}}}} \right)}^{2n}}}\; \bigg]\\ &~~\times\exp\left( { - 2\ln 2\dfrac{{{t^2}}}{{\tau _{{\text{FWHM}}}^2}}} \right)\exp \left( {{\text{i}}\dfrac{{2{\text{π }}c}}{{{\lambda _0}}}t} \right), \end{split}$$

(18)

式中, $R$为光束半径, $n$为超高斯光束阶数, ${\tau _{{\text{FWHM}}}}$为输入脉冲强度半峰全宽, ${\lambda _0}$为超短脉冲的中心波长, $c$为光速.

计算超短脉冲聚焦时, 首先通过时频变换得到超短脉冲中的光谱分布, 即:

$${U_{{\text{in}}}}\left( {x,y,\lambda } \right) = {\Im _t}\left[ {{E_{{\text{in}}}}(x,y,t)} \right].$$

(19)

然后, 利用(9)式计算出每一光谱分量聚焦场${U_{{\text{out}}}}\left( {x', y', \lambda } \right)$, 再作时频逆变换得到最终的超短脉冲聚焦后的时空分布:

$${E_{{\text{out}}}}(x',y',t) = \Im _t^{ - 1}\left\{ {{U_{{\text{out}}}}\left( {x',y',\lambda } \right)} \right\}.$$

(20)

参考OPCPA实验^[4]中的光脉冲情况, 模拟计算采用的具体参数由表1列出.

对于20 fs的超短脉冲, 计算出的光谱如图1所示, 其波长范围很宽, 从0.74—0.86 μm均有能量分布.

图2和图3分别给出(1)式传统衍射算法和本文算法的不同波长分量聚焦场计算结果. 如图2所示, 传统衍射算法结果中, 大量采样点光场强度几乎为零, 光场分辨率低, 特别是当输入光场波前存在畸变的情况下, 这样的计算结果不能清楚地反映波前畸变对聚焦光斑的影响; 另外, 输出光场范围依赖于波长, 各光谱分量网格大小不一样, 导致(20)式的傅里叶逆变换难以直接进行, 而此时若是对输出光场进行重采样将会大大增加整个系统计算量.

从图3可看出, 使用本文算法后, 可以灵活调整输出光场计算范围, 只针对能量不为零的区域进行计算. 这里将各波长分量的输出场范围统一为边长$102.4{\kern 1 pt} {\kern 1 pt} {\kern 1 pt} {\text{μm}}$的正方形, 不但提高分辨率, 而且有利于后续直接进行时频变换来分析脉冲聚焦后的时空特性. 另一方面可以由(16)式计算出聚焦场采样间隔${\text{δ}}x' = 0.05{\kern 1 pt} {\kern 1 pt} {\kern 1 pt} {\kern 1 pt} {\text{μm}}$对应的输入场有效带宽, 在这个带宽范围内, 各波长光束能量占比均超过总能量的99.8%, 由此保证了计算精度.

图4给出了超短脉冲聚焦后的时空分布图, 图4(a)为整体脉冲光的空间强度分布, 其$x'$轴上的一维强度分布如图4(b)中实线所示, 虚线对应中心波长分量情况, 此处, 两条曲线相对自身最大强度分别进行归一化处理. 通过对比可知, 在中心亮斑和各阶亮环附近, 整体脉冲光和中心波长分量分布比较接近、无明显区别, 而各谱分量暗环出现位置略有差异, 导致整体脉冲光的暗环相对强度和中心波长分量相比有明显的提高. 图4(c)和图5(d)分别为$y' = 0{\kern 1 pt} {\kern 1 pt} {\kern 1 pt} {\text{μm}}$和$y' = 8.0{\kern 1 pt} {\kern 1 pt} {\kern 1 pt} {\kern 1 pt} {\text{μm}}$时, $x'$轴上不同位置的时间分布图, 可以看出, 在暗环附近超短脉冲聚焦后会产生畸变, 原因在于暗环附近的光谱分布和传输之前相比发生很大变化. 为了更清楚地反映暗环附近的时空耦合演变规律, 图5给出了$x'$轴上暗环附近不同坐标点的归一化时间强度分布图, 其中实线为超短脉冲聚焦后的时间分布, 作为对比, 虚线是输入脉冲的时间分布. 从图5可以看出, 在暗环处, 脉冲时域波形分裂成两个波峰, 而且在零时刻强度最弱, 随着逐渐偏离暗环区域, 两个波峰逐渐演变合成为单一短脉冲. 值得说明的是, 尽管聚焦脉冲在暗环区域会发生严重畸变, 但由于脉冲暗环区域功率占总功率比例很小, 对光脉冲的整体时间分布影响却不大.

4.   总　结

本文提出了一种针对超短脉冲聚焦场的衍射传输算法, 并推导给出了具备快速傅里叶算法的数学表达式. 该方法优点表现在, 可以根据实际需求灵活调整网格大小, 提高聚焦场分辨率, 克服了传统衍射算法各谱分量输出场计算范围不一致的缺点, 为后续直接进行时频变换带来便利. 同时本文还讨论了使用该算法时需要满足的采样条件, 为数值计算时采用合适的窗口大小和采样间隔提供依据, 避免计算失真. 利用该算法对超短脉冲聚焦进行模拟计算, 结果发现超短脉冲聚焦场暗环区域表现出强烈的时空耦合特征. 本算法已成功应用于OPCPA数值仿真平台研发中, 可以用来分析超短脉冲聚焦过程发生的时空耦合效应.