圆柱形原子团簇中等离激元的本征模式和本征电荷分布

The Eigen-Mode and Eigen-Charge of Plasmons in Cylindrical Atomic Clusters

摘要: 基于自由电子气体模型和无规相近似下的线性响应理论，我们给出了寻找圆柱形原子团簇中可维持的等离激元模式的方法。本征函数被用来寻找等离激元，电荷分布被用来展现等离激元。此外，我们也给出了寻找被外场激发的等离激元模式的方法。

Abstract: The computational method of plasmons in cylinder is studied, based on the linear response theory in the random-phase approximation and the free-electron gas model. The eigen-function is offered to find the plasmons maintained. The energy absorption spectrum is offered to find the plasmons excited by the external fields.

1. 引言

21世纪是信息时代，快速、高密度和大容量是其显著特征。光子学器件和电子器件相比，其有无法比拟的低能耗和高速等特点，将成为信息时代的基石。随着器件集成化的发展趋势，器件将被做的越来越小。这必然要求器件的尺寸和空间距离均要突破光学衍射极限的限制。因受衍射极限的限制，传统光学的基本原理和技术难以为纳米尺度上的光子操控提供有效指引。表面等离激元是光与金属表面的电子相互作用，形成的一种表面电磁模式，其横向可以实现对电磁场亚波长尺度的束缚，从根本上突破了衍射极限的限制 [1] [2] [3] 。表面等离激元一般可以通过光学耦合或电子束轰击的方法激发 [4] [5] [6] 。目前，基于表面等离激元的纳米激光器 [7] 、生物探测和传感技术 [8] 以及突破衍射极限的光刻技术 [9] 等都已被开发出来。

最近，基于自由电子气体模型和无规相近似下的线性响应理论，我们提供了一种研究等离激元的新方法 [10] [11] 。利用这种新方法，我们发现在耦合原子链中，纵向等离激元随原子链长的增大而发生红移，横向等离激元随原子链长的增加而发生蓝移；纵向等离激元的电荷分布成体分布特征，而横向等离激元则分化成了中心模式(电荷主要分布在原子链的中心)和两端模式(电荷主要分布在原子链的两端) [11] 。用新方法得到的这些研究结果和以前人们的研究结果一致 [12] [13] ，这充分说明了新方法的有效性和正确性。

本文中，我们把研究等离激元的新方法推广到圆柱系统。在圆柱中，电子的本征函数是个复变函数，这导致圆柱中等离激元模式的查找要比耦合原子链中等离激元的查找要复杂的多。本文中，我们推导了圆柱系统中能量损耗谱和电荷密度的具体表达式。

2. 研究方法

由含时密度泛函理论可知，一个系统因外界扰动所产生的诱导电荷密度可以写成

${\rho }^{ind}\left(r,\omega \right)={\displaystyle \int \prod \left(r,r^{\prime },\omega \right)V\left(r^{\prime },\omega \right)\text{d}{r}^{\prime }}$ (1)

其中 $\prod \left(r,r^{\prime },\omega \right)$ 是密度-密度响应函数，在无规则相近似下，其可写为

$\prod \left(r,r^{\prime },\omega \right)=2e^2{\displaystyle \underset{nm}{\sum }\frac{f\left(E_n\right)-f\left(E_m\right)}{E_n-E_m-\omega -i\gamma }{\psi }_n^{\ast }\left(r\right){\psi }_m\left(r\right){\psi }_m^{\ast }\left(r^{\prime }\right){\psi }_n\left(r^{\prime }\right)}$ (2)

其中， $f\left(E_n\right)$ 是费米分布函数。本文我们只考虑零温情况，在费米能级以下费米分布函数为1，在费米能级以上费米分布函数为0。 ${\psi }_n$ 和 $E_n$ 分别代表系统未受外界扰动时的本征函数和本征能量。V为作用于系统的总扰动势，其可看做由外加电势 $V^{ext}$ 和诱导电势组成，即 $V=V^{ext}+V^{ind}$ 。

由库伦定律知，诱导电荷在空间产生的电势为

$V^{ind}\left(r,\omega \right)=\frac{1}{4\text{π}{\epsilon }_0}{\displaystyle \int \frac{{\rho }^{ind}\left(r^{\prime },\omega \right)}{\left|r-r^{\prime }\right|}\text{d}{r}^{\prime }}$ (3)

将电荷密度(1)式代入上式，可得

$\begin{array}{c}{V}^{ind}\left(r,\omega \right)=2e^2{\displaystyle \underset{nm}{\sum }\frac{f\left(E_n\right)-f\left(E_m\right)}{E_n-E_m-\omega -i\gamma }}\\ \times \frac{1}{4\text{π}{\epsilon }_0}{\displaystyle \int \frac{{\psi }_n^{\ast }\left(r^{\prime }\right){\psi }_m\left(r^{\prime }\right)}{\left|r-r^{\prime }\right|}\text{d}{r}^{\prime }}\left[V_{nm}^{ext}\left(\omega \right)+V_{nm}^{ind}\left(\omega \right)\right]\end{array}$ (4)

其中，

$V_{nm}\left(\omega \right)={\displaystyle \int {\psi }_n^{\ast }\left(r\right)V\left(r,\omega \right){\psi }_m\left(r\right)}$ (5)

在(4)式两边，同乘以 ${\psi }_{n^{\prime }}^{\ast }\left(r\right){\psi }_{m^{\prime }}\left(r\right)$ 并对整个系统积分，可得：

$\begin{array}{c}{V}_{n^{\prime }{m}^{\prime }}^{ind}\left(\omega \right)=2e^2{\displaystyle \underset{nm}{\sum }\frac{f\left(E_n\right)-f\left(E_m\right)}{E_n-E_m-\omega -i\gamma }}\\ \times \frac{1}{4\text{π}{\epsilon }_0}{\displaystyle \int \frac{{\psi }_n^{\ast }\left(r^{\prime }\right){\psi }_m\left(r^{\prime }\right){\psi }_{n^{\prime }}^{\ast }\left(r\right){\psi }_{m^{\prime }}\left(r\right)}{\left|r-r^{\prime }\right|}\left[V_{nm}^{ext}\left(\omega \right)+V_{nm}^{ind}\left(\omega \right)\right]\text{d}{r}^{\prime }\text{d}r}\end{array}$ (6)

令

$M_{n^{\prime }{m}^{\prime },nm}\left(\omega \right)=\frac{2e^2}{4\text{π}{\epsilon }_0}{\displaystyle \underset{nm}{\sum }\frac{f\left(E_n\right)-f\left(E_m\right)}{E_n-E_m-\omega -i\gamma }}{\displaystyle \int \frac{{\psi }_{n^{\prime }}^{\ast }\left(r\right){\psi }_{m^{\prime }}\left(r\right){\psi }_n^{\ast }\left(r^{\prime }\right){\psi }_m\left(r^{\prime }\right)}{\left|r-r^{\prime }\right|}\text{d}{r}^{\prime }\text{d}r}$ (7)

则(6)式可写成

${\displaystyle \underset{nm}{\sum }\left[{\delta }_{n^{\prime }{m}^{\prime },nm}-M_{n^{\prime }{m}^{\prime },nm}\left(\omega \right)\right]V_{nm}^{ind}\left(\omega \right)}={\displaystyle \underset{nm}{\sum }{M}_{n^{\prime }{m}^{\prime },nm}\left(\omega \right)V_{nm}^{ext}\left(\omega \right)}$ (8)

当外加电场知道时，根据电场和电势的关系 $V^{ext}\left(r,\omega \right)=-{\displaystyle {\int }_r^{r_0}E\left(r,\omega \right)\cdot \text{d}r}$ 可以获得 $V^{ext}\left(\omega \right)$ ( $r_0$ 是电势参考点)，进而利用(5)式可获得 $V_{nm}^{ext}\left(\omega \right)$ 。当 $V_{nm}^{ext}\left(\omega \right)$ 知道时，利用(8)式可获得 $V_{nm}^{ind}\left(\omega \right)$ ，进而利用(1)式可获得诱导电荷密度。

等离激元是电荷的集体振荡，当其激发时会使系统的能量损耗增大，从而导致能量损耗谱中出现损耗峰。因而，我们可以认为能量损耗(Energy absorption)峰所对应的频率就是等离激元的频率。通常情况下，一个系统的能量损耗可以通过下式来计算

$G\left(\omega \right)=-\frac{1}{2}\omega \mathrm{Im}\left\{{\displaystyle \int {\rho }^{ind}\left(r,\omega \right){\left[V^{ext}\left(r,\omega \right)\right]}^{\ast }}\right\}$ (9)

下面，我们来具体推导圆柱系统的能量损耗谱和电荷分布。在自由电子气体模型下，圆柱内电子的本征函数和本征能量分别为

${\psi }_{nkl}\left(r,\theta ,z\right)=\sqrt{\frac{2}{\text{π}{L}_0{R}_{\text{0}}^{\text{2}}}}\frac{J_n\left(\frac{x_{nk}}{R_0}r\right)}{J_{n-1}\left(x_{nk}\right)}\sin \left(\frac{l\text{π}}{L_0}z\right){\text{e}}^{\pm in\theta }$ (10)

$E_{nkl}=\frac{{\text{π}}^2{\hslash }^2}{2m{}_{e}a{}^2}\left[\frac{x_{nk}^2}{{\text{π}}^2{R}^2}+\frac{l^2}{L_0^2}\right]$ (11)

这里，圆柱长为 $L_0=L\times a$ ，圆柱半径为 $R_0=R\times a$ ，a代表晶格常数， $J_n$ 代表n阶贝塞尔函数， $x_{nk}$ 代表n阶贝赛尔函数的第k个零点， $m_e$ 是电子的质量，l是圆柱轴向的量子数。在后面推导中，为了表示方便，电子的波函数被简写成

${\psi }_{nkl}\left(r,\theta \right)=C{R}_n\left(N_{k}r\right)\sin \left(G_{l}z\right){\text{e}}^{\pm in\theta }$ (12)

这里 $C=\sqrt{2/\text{π}{L}_0{R}_{\text{0}}^{\text{2}}}$ ， $N_k={x_{nk}/R}_{\text{0}}$ ， $R_n\left(N_{k}r\right)=J_n\left(\frac{x_{nk}}{R_{\text{0}}}r\right)/J_n\left(x_{nk}\right)$ ， $G_l=l\text{π}/L_0$ 。将(12)式代入(1)式，

并利用电荷密度和电势的展开式

${\rho }^{ind}\left(r,\theta ,\omega \right)={\displaystyle \underset{m}{\sum }{\rho }_m^{ind}\left(r,\omega \right){\text{e}}^{im\theta }}$ (13)

$V^{ind}\left(r,\theta ,\omega \right)={\displaystyle \underset{m}{\sum }{V}_m^{ind}\left(r,\omega \right){\text{e}}^{im\theta }}$ (14)

可得

$\begin{array}{c}{\rho }_m^{ind}\left(r,\omega \right)=4\text{π}{e}^2\left[{\displaystyle \underset{nklm{k}^{\prime }{l}^{\prime }}{\sum }\left(\frac{1}{E_{nkl}-E_{\left(n+m\right)k^{\prime }{l}^{\prime }}-\omega -i\gamma }+\frac{1}{E_{nkl}-E_{\left(n+m\right)k^{\prime }{l}^{\prime }}+\omega +i\gamma }\right)}\right]\\ \times C^2{R}_n\left(N_{k}r\right)R_{n+m}\left(N_{k^{\prime }}r\right)\sin \left(G_{l}z\right)\sin \left(G_{l^{\prime }}z\right)V_{nklm{k}^{\prime }{l}^{\prime }}\left(\omega \right)\end{array}$ (15)

$V_m^{ind}\left(r,z,\omega \right)={\displaystyle \int \frac{r^{\prime }{\rho }_m^{ind}\left(r^{\prime },z^{\prime },\omega \right){\text{e}}^{im\left({\theta }^{\prime }-\theta \right)}}{4\text{π}{\epsilon }_0\sqrt{r^2+{r^{\prime }}^2-2r{r}^{\prime }\cos \left(\theta -{\theta }^{\prime }\right)+{\left(z-z^{\prime }\right)}^2}}\text{d}{r}^{\prime }\text{d}{\theta }^{\prime }\text{d}{z}^{\prime }}$ (16)

这里

$V_{nklm{k}^{\prime }{l}^{\prime }}\left(\omega \right)=C^2{\displaystyle \int R_n\left(N_{k}r\right)\sin \left(G_{l}z\right)V_m\left(r,z,\omega \right)R_{n+m}\left(N_{k^{\prime }}r\right)\sin \left(G_{l^{\prime }}z\right)r\text{d}r\text{d}z}$ (17)

将电荷密度(15)式代入电势方程(16)式后，两边同乘以 $C^{2}r{R}_{n_1+m}\left(N_{{k^{\prime }}_1}r\right)R_{n_1}\left(N_{k_1}r\right)\sin \left(G_{{l^{\prime }}_1}z\right)\sin \left(G_{l_1}z\right)$ ，并对r，z积分，则有

$\begin{array}{c}{V}_{n_1{k}_1{l}_{1}m{k^{\prime }}_1{l^{\prime }}_1}^{ind}\left(\omega \right)=4\text{π}{e}^2{C}^4{\displaystyle \underset{nklm{k}^{\prime }{l}^{\prime }}{\sum }\left(\frac{1}{E_{nkl}-E_{\left(n+m\right)k^{\prime }{l}^{\prime }}-\omega -i\gamma }+\frac{1}{E_{nkl}-E_{\left(n+m\right)k^{\prime }{l}^{\prime }}+\omega +i\gamma }\right)}\\ \times {\displaystyle \int [\frac{r{r}^{\prime }{\text{e}}^{im\left({\theta }^{\prime }-\theta \right)}{R}_{n_1+m}\left(N_{{k^{\prime }}_1}r\right)R_{n_1}\left(N_{k_1}r\right)R_n\left(N_k{r}^{\prime }\right)R_{n+m}\left(N{}_{k^{\prime }}{r}^{\prime }\right)}{\text{4π}{\epsilon }_{\text{0}}\sqrt{r^2+{r^{\prime }}^2-2r{r}^{\prime }\cos \left(\theta -{\theta }^{\prime }\right)+{\left(z-z^{\prime }\right)}^2}}}\\ \underset{\stackrel{}{{\displaystyle }}}{\overset{\stackrel{}{{\displaystyle }}}{}}\times \sin \left(G_{{l^{\prime }}_1}z\right)\sin \left(G_{l_1}z\right)\sin \left(G_l{z}^{\prime }\right)\sin \left(G_{l^{\prime }}{z}^{\prime }\right)\text{d}r\text{d}z\text{d}{r}^{\prime }\text{d}{z}^{\prime }\text{d}{\theta }^{\prime }]\times \left[V_{nkm{k}^{\prime }{l}^{\prime }}^{ind}+V_{nkm{k}^{\prime }{l}^{\prime }}^{ext}\right]\end{array}$ (18)

令

$\begin{array}{c}{M}_{n_1{k}_1{l}_{1}m{k^{\prime }}_1{l^{\prime }}_1,nklm{k}^{\prime }{l}^{\prime }}\left(\omega \right)=4\text{π}{e}^2{C}^4\left(\frac{1}{E_{nkl}-E_{\left(n+m\right)k^{\prime }{l}^{\prime }}-\omega -i\gamma }+\frac{1}{E_{nkl}-E_{\left(n+m\right)k^{\prime }{l}^{\prime }}+\omega +i\gamma }\right)\\ \times {\displaystyle \int [\frac{r{r}^{\prime }{\text{e}}^{im\left({\theta }^{\prime }-\theta \right)}{R}_{n_1+m}\left(N_{{k^{\prime }}_1}r\right)R_{n_1}\left(N_{k_1}r\right)R_n\left(N_k{r}^{\prime }\right)R_{n+m}\left(N{}_{k^{\prime }}{r}^{\prime }\right)}{\text{4π}{\epsilon }_{\text{0}}\sqrt{r^2+{r^{\prime }}^2-2r{r}^{\prime }\cos \left(\theta -{\theta }^{\prime }\right)+{\left(z-z^{\prime }\right)}^2}}}\\ \underset{\stackrel{}{{\displaystyle }}}{\overset{\stackrel{}{{\displaystyle }}}{}}\times \sin \left(G_{{l^{\prime }}_1}z\right)\sin \left(G_{l_1}z\right)\sin \left(G_l{z}^{\prime }\right)\sin \left(G_{l^{\prime }}{z}^{\prime }\right)\text{d}r\text{d}z\text{d}{r}^{\prime }\text{d}{z}^{\prime }\text{d}{\theta }^{\prime }]\end{array}$ (19)

则

${\displaystyle \underset{nklm{k}^{\prime }{l}^{\prime }}{\sum }\left[{\delta }_{n_1{k}_1{l}_{1}m{k^{\prime }}_1{l}^{\prime },nklm{k}^{\prime }{l}^{\prime }}-M_{n_1{k}_1{l}_{1}m{k^{\prime }}_1{l^{\prime }}_1,nklm{k}^{\prime }{l}^{\prime }}\left(\omega \right)\right]V_{nklm{k}^{\prime }{l}^{\prime }}^{ind}\left(\omega \right)}={\displaystyle \underset{nklm{k}^{\prime }{l}^{\prime }}{\sum }{M}_{n_1{k}_1{l}_{1}m{k^{\prime }}_1{l^{\prime }}_1,nklm{k}^{\prime }{l}^{\prime }}\left(\omega \right)V_{nklm{k}^{\prime }{l}^{\prime }}^{ext}\left(\omega \right)}$ (20)

由(14)式可知

$V_m^{ext}\left(r,\omega \right)=\frac{1}{\text{2π}}{\displaystyle \int V^{ext}\left(r,\theta ,\omega \right)\cos \left(m\theta \right)\text{d}\theta }$ (21)

当外加电势 $V^{ext}\left(r,\theta ,\omega \right)$ 已知时，通过(21)式可求得 $V_m^{ext}\left(r,\omega \right)$ ；将 $V_m^{ext}\left(r,\omega \right)$ 代入(17)式，进而可求得 $V_{nklm{k}^{\prime }{l}^{\prime }}^{ext}\left(\omega \right)$ ；将 $V_{nklm{k}^{\prime }{l}^{\prime }}^{ext}\left(\omega \right)$ 代入(20)式，可求得 $V_{nklm{k}^{\prime }{l}^{\prime }}^{ind}\left(\omega \right)$ 。当获得 $V_{nklm{k}^{\prime }{l}^{\prime }}^{ext}\left(\omega \right)$ 和 $V_{nklm{k}^{\prime }{l}^{\prime }}^{ind}\left(\omega \right)$ 后，利用(15)式可获得诱导电荷密度 ${\rho }_m^{ind}\left(r,\omega \right)$ ，进而通过(13)式获得等离激元激发时的电荷分布。当 $V_m^{ext}\left(r,\omega \right)$ 和 ${\rho }_m^{ind}\left(r,\omega \right)$ 知道后，利用(21)式可获得能量损耗谱。

当外加电势 $V^{ext}\left(r,\theta ,\omega \right)=0$ 时，利用(20)式可以得到 $e^{im\theta }$ 分量上的等离激元的本征方程，即：

${\displaystyle \underset{nklm{k}^{\prime }{l}^{\prime }}{\sum }\left[{\delta }_{n_1{k}_1{l}_{1}m{k^{\prime }}_1{l}^{\prime },nklm{k}^{\prime }{l}^{\prime }}-M_{n_1{k}_1{l}_{1}m{k^{\prime }}_1{l^{\prime }}_1,nklm{k}^{\prime }{l}^{\prime }}\left(\omega \right)\right]V_{nklm{k}^{\prime }{l}^{\prime }}^{ind}\left(\omega \right)}=0$ (22)

等离激元的本征频率可通过令上式中的系数行列式为零而确定，即：

$A\left(\omega \right)=\det \left[{\delta }_{n_1{k}_1{l}_{1}m{k^{\prime }}_1{l}^{\prime },nklm{k}^{\prime }{l}^{\prime }}-M_{n_1{k}_1{l}_{1}m{k^{\prime }}_1{l^{\prime }}_1,nklm{k}^{\prime }{l}^{\prime }}\left(\omega \right)\right]=0$ (23)

3. 计算结果

在实际计算中，(19)式中的虚部 $i\gamma$ 通常是不为零的。这是因为 $i\gamma =0$ 暗示示电子在振荡过程中不受阻尼，进而导致等离激元的寿命为无限大。在实际中，等离激元的寿命都是有限的，所以要取一个不为零的 $i\gamma$ ，以表示等离激元寿命的有限性和电子在振荡过程中有阻尼存。这样，本征函数 $A\left(\omega \right)$ 通常情况下都会是个复函数。利用(23)式求等离激元的频率时，不但保证本征函数 $A\left(\omega \right)$ 的实部 $\mathrm{Re}\left(A\left(\omega \right)\right)=0$ ，还要保证本征函数 $A\left(\omega \right)$ 的虚部 $\text{Im}\left(A\left(\omega \right)\right)\approx 0$ 。这样， $\left|\text{Im}\left(1/A\left(\omega \right)\right)\right|$ 在等离激元频率处会显示一个有限峰，可以利用这个峰来寻找等离激元。

图1中展示了 $m=0$ 时，本征函数对等离激元本征模式的预言，这里圆柱长为 $L_0=\text{6}\times a$ ，圆柱半径

Figure 1. The spectrum functions $\left|\text{Im}\left(1/A\left(\omega \right)\right)\right|$ and $\mathrm{Re}\left(A\left(\omega \right)\right)$ change with frequency

图1. 普函数 $\left|\text{Im}\left(1/A\left(\omega \right)\right)\right|$ 和 $\mathrm{Re}\left(A\left(\omega \right)\right)$ 随频率的变化

Figure 2. The eigen-charge distributions

图2. 本征电荷分布

为 $R_0=\text{1}\times a$ ，频率 $\omega$ 是以 ${\text{π}}^2{\hslash }^2/2m_e{a}^2$ 位单位， $m_e$ 代表电子的质量， $a=0.2857\text{nm}$ 代表金的晶格常数。这里普函数 $\left|\text{Im}\left(1/A\left(\omega \right)\right)\right|$ 和普函数 $\mathrm{Re}\left(A\left(\omega \right)\right)$ 被结合起来以寻找等离激元的本征模式。从图中可以看出在谱函数 $\left|\text{Im}\left(1/A\left(\omega \right)\right)\right|$ 的峰频处，普函数 $\mathrm{Re}\left(A\left(\omega \right)\right)$ 也为零。这说明普函数 $\left|\text{Im}\left(1/A\left(\omega \right)\right)\right|$ 的峰对应的频率就是等离激元的本征频率。

图2展示了在本征频率处，圆柱轴向积分后的本征电荷分布。从图中可以看到，本征电荷的分布显示出了明显的角对称特征，这种特征在其它的研究中也被发现 [14] 。其次，还可发现本征电荷的实部和虚部特性相同，但是它们存在一定的相位差。

4. 小结

本文提供了圆柱中等离激元的计算方法。有外场时，能量损耗谱和电荷分布被分别用来寻找激发的等离激元模式和展现相应模式的电荷分布。无外场时，本征函数谱和本征电荷分别被用来寻找系统可维持的等离激元模式和响应的电荷分布。本文中，我们主要讨论了系统可维持的等离激元模式。结果发现，本征电荷分布展现出了明显的角对称特征。

基金项目

本文为陕西省科技厅科研项目(2018JQ1091)，西安航空学院校级项目(2017KY0209)和西安航空学院横向项目(规则纳米器件中等离激元的模式及其激发方式研究)的资助课题。

NOTES

^*通讯作者。