闪耀光栅 (DMD) 的衍射效应

摘要

这篇博客主要介绍DMD作为闪耀光栅的衍射效应。首先我们介绍单缝衍射、多缝衍射和光栅方程，然后讨论闪耀光栅，最后仿真给出不同入射角下的闪耀判据和光栅光谱的二维模拟。

DMD微镜翻转状态

现在wavefront shaping (WFS) 的实验多采用数字微镜阵列(DMD)，可用于快速波前整形。DMD是一种二元衍射光学元件，本质是个反射型闪耀光栅，可用于振幅的周期性调制。以DLP7000[1]芯片为例，其$1024\times 768$的微镜阵列相对flat state可翻转 $\pm 12^0$（对应“on-state”和“off-state”），刷新率可高达23KHz，见下图[1]：

单缝衍射

先简单提一下单缝衍射。考虑如下装置，单缝宽为a，在θ方向可以将单缝面元n等分，则相邻面元之间的光程差为$\Delta \delta =asin\theta /n$，相位差为$\Delta \varphi =\frac{2\pi asin\theta }{\lambda n}$，所有面元在$F$点产生的扰动振幅为$\Delta A\frac{sin(n\Delta \varphi /2)}{sin(\Delta \varphi /2)}$。由于$\Delta \varphi$很小，有$sin(\Delta \varphi )\approx \Delta \varphi$。因此$F$点振幅为$A\frac{sin(\alpha )}{\alpha }$，其中$A=n\Delta A$，$\alpha =\pi asin\theta /\lambda$.

多缝衍射 (光栅衍射)[2]

对于光栅衍射，多缝按缝距d将入射光波前分为N个部分，每个部分成为缝宽为a的单缝而发生夫琅禾费衍射。这样多个单缝衍射复振幅场发生缝间干涉而形成多缝衍射。

由单缝衍射理论可知，对于波长为$\lambda$、$\theta$角入射光，中心单缝衍射在透镜后焦面上一点P的复振幅为
$$\tilde{E}(P)=A(\frac{sin\alpha }{\alpha })$$其中振幅$A$为常数，$\alpha =\frac{\pi }{\lambda }asin\theta$。则相邻单缝在P点的衍射幅值与中心单缝的相同，只是产生一个相位差：
$$\delta =\frac{2\pi }{\lambda }dsin\theta$$由此多缝在P点的复振幅为N个振幅相同、相邻相位差恒定的多光束干涉结果。
E ~ ( P ) = A ( s i n α α ) { 1 + e x p ( i δ ) + . . . + e x p ( i ( N − 1 ) δ } = A ( s i n α α ) ( s i n N δ / 2 s i n δ / 2 ) e x p ( i ( N − 1 ) δ / 2 ) \widetilde{E}(P)=A(\frac{sin\alpha}{\alpha})\{1+exp(i\delta)+...+exp(i(N-1)\delta\}\\ =A(\frac{sin\alpha}{\alpha})(\frac{sinN\delta/2}{sin\delta/2})exp(i(N-1)\delta/2) E (P)=A(αsinα​){1+exp(iδ)+...+exp(i(N−1)δ}=A(αsinα​)(sinδ/2sinNδ/2​)exp(i(N−1)δ/2)

因此P点光强为$$I(P)=I_0(\frac{sin\alpha }{\alpha }{)}^2(\frac{sinN\delta /2}{sin\delta /2}{)}^2$$

其中$I_0={\left|A\right|}^2$是单缝在P点的光强，$(\frac{sin\alpha }{\alpha }{)}^2$为单缝衍射因子，$(\frac{sinN\delta /2}{sin\delta /2}{)}^2$为多光束干涉因子。

分析
每个单缝的衍射光强来自于各单缝的光振幅矢量$I_0$的大小，并随衍射角 $\theta$而变化。而多缝干涉主极大的光强为 $N^2{I}_0(\frac{sin\alpha }{\alpha }{)}^2$，因此光栅衍射图样是多缝干涉光强分布受单缝衍射光强分布调制的结果。

光栅方程

衍射光栅可对入射光波的振幅或相位进行空间周期性调制，按工作方式可分为透射式和反射式。光栅最重要的性能就是用作分光元件，值得一提的是近些年空间光调制器 (SLM, DMD & DM)用于波前整形有系列应用前景。一般讨论光栅的分光性能包括光栅方程、光栅的色散和色分辨本领，这里单讲光栅方程。

由光栅衍射的多光束干涉因子可知，当
$$\delta =\frac{2\pi }{\lambda }dsin\theta =2m\pi m=0,\pm 1,\pm 2,\cdots$$即$dsin\theta =m\lambda m=0,\pm 1,\pm 2,\cdots$，相邻光束相位差为$2\pi$整数倍时发生干涉叠加，形成主极大，该式子决定各级干涉主极大位置，称为光栅方程。

上式为垂直光栅入射的情形，下面分析更为一般的斜入射情形，导出光栅方程。

平行光以 $\alpha$ 角斜入射时，考察的衍射光与入射光分别位于光栅法线的两侧（左图）或同侧（右图）。无论光栅为透射式还是反射式，两相邻光束的光程差为$\Delta =dsin\alpha \pm dsin\beta$，因此光栅方程的一般形式为$$d(sin\alpha \pm sin\beta )=m\lambda m=0,\pm 1,\pm 2,\cdots$$若衍射光与入射光在光栅法线同侧上式取正号，异侧取负号。

事实上，斜入射时多缝夫琅禾费衍射的强度分布公式[3]也随之为：$$I_{(}{\theta }_i,{\theta }_o)=I_0(\frac{sin{\alpha }^{\prime }}{{\alpha }^{\prime }}{)}^2(\frac{sinN{\delta }^{\prime }/2}{sin{\delta }^{\prime }/2}{)}^2$$其中 ${\alpha }^{\prime }=\frac{\pi a}{\lambda }(sin{\theta }_i\pm sin{\theta }_o)$, ${\delta }^{\prime }=\frac{2\pi d}{\lambda }(sin{\theta }_i\pm sin{\theta }_o)$

闪耀光栅[2,3]

前面讲到的光栅有一个很大缺点，作为分光元件，无色散的0级主极大占了总光能很大一部分，而其余每级光谱的强度很弱。这是由于单缝衍射中央极大与缝间干涉0级重叠造成的。实际上光栅使用时可能需要利用某一级干涉光谱。下面介绍的闪耀光栅能把单缝衍射0级与缝间干涉0级错开，将光能从干涉0级光谱转移并集中到所需的某一级光谱上，实现该级光谱的闪耀。

目前闪耀光栅多是平面反射光栅。在平坦光栅面上刻划出一系列锯齿状槽面。刻槽面与光栅面的夹角称为闪耀角 ${\theta }_B$，如下图[4]所示

有两种典型照明方式：

• 平行光沿刻槽面法线入射（Littrow自准式入射）
  这时单缝衍射主极大方向为几何光学的反射方向，即沿着原方向返回。而对缝间干涉，以光栅面法线为基准，入射角和反射角都为 ${\theta }_B$。由斜入射情形的光栅方程可知，决定各级干涉主极大的式子（闪耀条件）为$$2dsin{\theta }_B=m\lambda$$

当 $m=1$时，满足该式子的${\lambda }_1$为1级闪耀波长，光栅的单缝衍射0级主极大正好落在${\lambda }_1$光波的干涉1级谱线上。此外，通过设计闪耀角，可以改变干涉级 $m$，使光栅适用于某一特定波段的某级光谱的闪耀。

• 平行光沿光栅面法线入射
  这时入射角为$0$，反射角为 $2{\theta }_B$，这时相邻两束光的光程差为$\Delta L=dsin2{\theta }_B$。类似的，可用斜入射情形的光栅方程分析该情况。

需要指出的是，
（1) 单缝衍射主极大方向就是入射光的反射方向，若该方向上相邻两束光的光程差满足光栅方程的 $m$ 级，就能实现单缝衍射0级与缝间干涉 $m$ 级的重合。
（2）在闪耀光栅不变且入射光为某一特定波段的情况下，由于上述照明方式规定了入射角，只可能激发某级光谱的闪耀；而通过改变入射角 $\alpha$ 可改变相邻光束相位差，实现不同干涉级光谱的闪耀。

Python仿真[5,6]

这里主要参考的是Sebastien在他的Wavefront Shaping网站给出的Setting up a DMD: Diffraction effects[5]，事实上也是这篇教程激励我去重新回顾光栅衍射的知识并总结成博客。感谢大佬为Wavefront Shaping社区开源贡献，也希望光学社区多学习一下CS的开源精神（狗头）

不同入射角下的闪耀判据

仿真时需要注意以下几点：

1. 观察最开始的DMD微镜翻转状态示意图，由于微镜朝45度轴向排列，因此正常在与微镜垂直的横断面（水平面）上定义的光栅闪耀角、入射角和反射角需要变换到45度轴面上去。如下图所示，可按照角度正切值除以$\sqrt{2}$的方式进行校正（投影）。注意：如果将DMD沿45度摆放，这时微镜翻转就是水平向左或向右12度，无须做角度校正。
   

2. 观察闪耀光栅的示意图，考虑衍射光与入射光在光栅面法线同/异侧光程差取正负号的问题，一般情形下的入射角 $\alpha$ 和反射角 $\beta$ 有如下关系：$$\alpha +\beta =2{\theta }_B$$ 这样两相邻光束的光程差可直接写为$\Delta =dsin\alpha +dsin\beta$

3. 考虑一般情形的闪耀条件：$$m=\frac{d(sin\alpha +sin(2{\theta }_B-\alpha ))}{\lambda }$$在闪耀光栅不变且入射光为某一特定波段时，通过改变入射角 $\alpha$ 可改变衍射方向上相邻光束的相位差，若 $m$ 为整数，可实现$m$级光谱的闪耀。因此可定义闪耀判据为：$$\mu =\left|m\%1-0.5\right|$$ $\mu$ 在$0\sim 0.5$之间，当$\mu$接近 $0.5$ 时， $m$ 接近整数达到闪耀条件，该级光谱闪耀；当$\mu$接近 $0$ 时， $m$ 接近整数$+0.5$，光能分散在多级光谱上。

4. 由于对称性，仿真时只考虑X或Y轴的入射角，这样转变为一个1维问题。

以 $532nm$ 激光，DMD芯片DLP7000(微镜间距$d=13.68\mu m$，闪耀角${\theta }_B=12^0$)为例，不同入射角下的闪耀判据绘成曲线如下：
这里我们给出了几种入射情况：
a. $\alpha =0$，沿光栅面法线入射，此时衍射方向为从光栅面法线到刻槽面法线方向再偏离1个闪耀角；
b. $\alpha ={\theta }_B$，沿刻槽面法线入射，此时衍射方向为原路反向；
c. $\alpha =2{\theta }_B$，从光栅面法线到刻槽面法线方向再偏离1个闪耀角入射，此时衍射方向为沿光栅面法线出射
d. $\alpha =41.92^0$，做校准后${\alpha }_{1D}=32.41^0$，闪耀角${\theta }_{1D}=8.55^0$，该入射方向能实现$m=\frac{d(sin{\alpha }_{1D}+sin(2{\theta }_{1D}-{\alpha }_{1D}))}{\lambda }\approx 7$级干涉光谱的闪耀。这是该DMD对特定波段$532nm$激光在入射角由$0\sim 90^0$增加过程中的实现的第一次光谱闪耀，干涉级$m=7$。

光栅光谱的二维模拟

有了上述分析，这里直接给出入射角为 $\alpha =24^0$和 $\alpha =41.92^0$ 两种情况的DMD衍射光场分布和对应的光栅衍射图样。

（1）入射角为 $\alpha =24^0$，光能分散在各级光谱

（2）入射角为$\alpha =41.92^0$ ，第一次光谱闪耀，干涉级$m=7$

Source code[6]

# -*- coding: utf-8 -*-
# Written by Sebastien Popoff, adapted by Shengfu Cheng
# 08/06/2020
import matplotlib.pyplot as plt
import numpy as np

import matplotlib as mpl
from matplotlib import rc

## Parameters of DMD
d = 13.68 #micromirror pitch
gamma = 12./180*np.pi  #micromirror tilt angle relative to flat state

#Use sqrt(2) to calibrate the incident angle and tilt angle 
#the angle of rotation of the mirrors is at 45 degrees compared to the axis of the pixels
gamma_1D = np.arctan(np.tan(gamma)/np.sqrt(2))
_lambda = 532e-3 #wavelength

## A simple criterion matching the diffraction angle and the reflection angle
_beta = lambda x: 2*gamma_1D-x
m = lambda x: 1.*d/_lambda * (np.sin(x)+np.sin(_beta(x)))

## Test different incident angles
# when m is an integer, we are at a blazing angle with a maximum of energy at the order along the optical axis
# whem m is n+1/2, the enregy is spread over many diffraction orders not aligned with the optical axis
alpha_vec = np.linspace(-np.pi/2,np.pi/2,1000)
alpha_1D_vec = np.arctan(np.tan(alpha_vec)/np.sqrt(2))
criterion = np.abs(np.mod([m(a) for a in alpha_1D_vec],1)-0.5)
#find blaze angle
alpha_loc = np.linspace(0,5*np.pi/18,100)
blaze_angle = 0
for a in alpha_loc:
    if np.abs(np.mod(m(np.arctan(np.tan(a)/np.sqrt(2))),1)-0.5) > np.abs(np.mod(m(np.arctan(np.tan(blaze_angle)/np.sqrt(2))),1)-0.5):
        blaze_angle = a

plt.figure()
plt.plot(alpha_vec*180/np.pi,criterion,linewidth = 1)

plt.plot([0,0],[0.,0.5],color = 'C3',label = r'$\alpha = 0, \mu = %.4f$' \
    % np.abs(np.mod(m(0),1)-0.5), linewidth = 1) #At Normal line of grating surface 

plt.plot([gamma*180/np.pi,gamma*180/np.pi],[0.,0.5],color = 'C1',label = r'$\alpha = \theta_B, \mu = %.4f$' \
     % np.abs(np.mod(m(gamma_1D),1)-0.5),linewidth = 1) #At Normal line of groove surface 

plt.plot([2*gamma*180/np.pi,2*gamma*180/np.pi],[0.,0.5],color = 'C2',label = r'$\alpha = 2\theta_B, \mu = %.4f$' \
    % np.abs(np.mod(m(2*gamma_1D),1)-0.5), linewidth = 1) #2-time angles away from the normal of grating surface

plt.plot([blaze_angle*180/np.pi,blaze_angle*180/np.pi],[0.,0.5],color = 'C4',label = r'$\alpha = %.2f^{0}, \mu = %.4f$' \
    % (blaze_angle*180/np.pi, np.abs(np.mod(m(np.arctan(np.tan(blaze_angle)/np.sqrt(2))),1)-0.5)), linewidth = 1) #blaze angle

plt.xticks(np.arange(-90,90,10))
plt.yticks(np.arange(0,0.55,0.05))
plt.title(r'Blazing condition when $d=%.1f\mu m, \lambda =%g nm$' % (d,1e3*_lambda))
plt.legend(loc='center right')
plt.ylabel(r'Blazing criterion $\mu$')
plt.xlabel(r'Incident angle $\alpha$')

## A full numerical simulation of the Fourier plane for an all-on configuration
# alpha = 24/180*np.pi # Incident angle
alpha = blaze_angle
alpha_1D = np.arctan(np.tan(alpha)/np.sqrt(2))
beta = -alpha +2*gamma # reflection angle
alpha_1D = np.arctan(np.tan(alpha)/np.sqrt(2))
beta_1D = 2*gamma_1D-alpha_1D#np.arctan(np.tan(beta)/np.sqrt(2))
N = 20 # number of mirrors in each direction
g = 2 # gap between mirrors in micron
res = 10 # pixels per mirror 
Nx = N*res

## Pixelate image function
f = np.ones([N,N]) # all-on configuration
## Phase slope due to incident and reflection angle
X,Y = np.meshgrid(np.arange(N),np.arange(N))
phi = np.exp((X-Y)*1j*2*np.pi/_lambda*d*(np.sin(alpha_1D)+np.sin(beta_1D))) 

## cell unit
Cell = np.zeros([res,res])
gpix = int(np.round(g/(2.*d)*res))
Cell[gpix:res-gpix,gpix:res-gpix] = 1.
## Mirror image
MI = np.zeros([Nx,Nx],dtype='complex')
for i in range(N):
    for j in range(N):
        MI[i*res:(i+1)*res,j*res:(j+1)*res]= f[i,j]*phi[i,j]*Cell
plt.figure()
plt.imshow(np.real(MI),interpolation = 'None')

## In the Fourier plane
coeff = 5
FP = np.fft.fftshift(np.fft.fft2(MI,s=[coeff*Nx,coeff*Nx]))
ROIsize = 500
ROI = 1
plt.figure()
plt.imshow((np.abs(FP[coeff*Nx//2-ROIsize//2:coeff*Nx//2+ROIsize//2,coeff*Nx//2-ROIsize//2:coeff*Nx//2+ROIsize//2])),interpolation = 'None')
plt.clim([0,(np.max(np.abs(FP)))/1.5])
plt.scatter(ROIsize//2, ROIsize//2, s=200, edgecolors='red',alpha = 0.5,c='yellow',linewidths= 1.,  marker='x')

plt.show()

Reference

[1] DLP7000芯片手册
[2] 郁道银等.《工程光学》第4版
[3] 赵凯华 《新概念物理教程-光学》
[4] Wikipedia: Blazed_grating
[5] [tutorial] Setting up a DMD: Diffraction effects
[6] blazing_angle_DMD.ipynb