Nonlinear Plasmonics: Four-photon Near-field Photolithography using Optical Antennas

Abstract

Recent experiments on four-photon nonlinear exposure of photoresist near nanoplasmonic structures raise a question: What nonlinear processes are responsible for the observed profile? Here, we study the nonlinear exposure of poly(methyl methacrylate) (PMMA) in the near-field of gold nanoantennas. We consider six possible nonlinear processes and study them in terms of the developed volumes and the exposure profiles in photoresist. We find that the direct fourth harmonic generation (4HG) in gold is the dominant nonlinear process. The developed volume from 4HG and the exposure profiles both match closely with the experiments. The next strongest process is direct four-photon absorption (4PA) in PMMA. The strength of 4PA process is about one order of magnitude weaker than 4HG. The developed volume and exposure profiles predicted from 4PA process clearly deviate from experiments.

Appendix: FDTD Simulations

FDTD is a time domain simulation method to solve the Maxwell’s equations [49]. We implemented FDTD to calculate the near-field profile for all the processes using Lumerical FDTD Solutions (version 8.5). Our target is to obtain the intensity profile of all frequencies that may participate in the studied nonlinear process, i.e., $I_{\mathrm FF}(\mathbf {r})$ for $\omega _{0}$, $I_{\mathrm 2HG}(\mathbf {r})$ for $2\omega _{0}$, $I_{\mathrm 3HG}(\mathbf {r})$ for $3\omega _{0}$, and $I_{\mathrm 4HG}(\mathbf {r})$ for $4\omega _{0}$. For conveniently calculating exposure dose for varying laser power, we carried out frequency domain normalization in each simulation to obtain the field profile excited by a continuous harmonic wave with constant intensity $I_{0}=1\times 10^{10}~{\mathrm W/cm^{2}}$ and $\lambda =860~{\mathrm nm}$. Then the intensity for different laser powers can be scaled and the exposure dose can be calculated according to Table 1. Our fundamental methodology to simulate the nonlinear processes is: we first simulate the linear response, obtain the first-order field and calculate the higher-order nonlinear source terms; then we launch the source terms as electric dipoles in a second simulation and normalize the simulated field intensity properly. In the following, we present the details of calculating each nonlinear process.

Simulations for Process (I): Direct 4PA

To calculate the exposure dose under 4PA process in PMMA, we only need to simulate the linear response of gold nanoantennas. Through Fourier transform and source normalizations, we obtained the complex electric field $\tilde {\textbf {E}}{^{(1)}}(\textbf {r})$ for fundamental frequency $\omega _{0}$, where the superscript (1) represents the first order, i.e., linear response. The corresponding harmonic field in time domain is expressed as

$$\mathbf{E}^{(1)}(\mathbf{r},t)=\tilde{\mathbf{E}}^{(1)}(\mathbf{r})e^{-i\omega_{0}t}+c.c. $$

(18)

Under this convention, the amplitude of the harmonic field is $|2\tilde {\textbf {E}}^{(1)}(\textbf {r})|$ and the intensity profile $I_{\mathrm FF}(\textbf {r})=2n\epsilon _{0}c|\tilde {\textbf {E}}^{(1)}(\textbf {r})|^{2}$, where n is the refractive index of the photoresist.

Simulations for Process (II): Direct 4HG

Due to the inversion symmetry inside the bulk gold, the fourth-order nonlinear susceptibility $\chi ^{(4)}$ is zero in dipole approximation. However, the symmetry breaks at the surface potential of the gold and multiple harmonics can be generated from the gold surface [50–52]. To simulate the 4HG from the gold surface, we first obtained the linear field $\tilde {\textbf {E}}^{(1)}(\textbf {r})$ and then launched electric field dipole sources (center frequency $4\omega _{0}$) on the computational nodes located at the metal–dielectric boundary.

The nonlinear complex electric field dipole source is calculated by

$$ \tilde{E}_{\mathrm s\bot}^{(4)}(4\omega_{0})=\frac{\chi^{(4)}\left(\tilde{E}_{\bot}^{(1)}\right)^{4}}{\epsilon(4\omega_{0})-1}, $$

(19)

where the subscript $\mathrm \bot$ denotes the component that is at the interface and perpendicular to the surface. Because all the surfaces of gold antenna were meshed into staircase in FDTD simulations, $\mathrm \bot$ could be x, y, or z, as shown in Fig. 7. The time domain electric dipoles were launched with amplitude and phase angle determined from $\tilde {E}_{\mathrm s\bot }^{(4)}(4\omega _{0})$.

Fig. 7

Configurations of electric dipole sources on the boundary nodes for calculating the direct 4HG from gold surface. Because all the curvature of the rod was meshed into staircase shape in FDTD, the field components normal to the surface is pointing to the x, y or z direction. The polarization of each dipole has to be consistent with the field component at that node in the Yee’s cell configuration [49]. For example, the x-polarized dipole has to be added on the node point from which FDTD calculates $E_{x}$

In order to estimate the value of $\chi ^{(4)}$, we used our approach to simulate the 4HG excited by laser pulses on a planar gold surface. According to the theoretical results by Georges et al. [52], a peak intensity of $1\times 10^{10}~{\mathrm W/cm}^{2}$ can produce 4HG intensity about $10^{-1}~{\mathrm W/cm}^{2}$. We simulated the intensity of the reflected 4HG wave, compared our results to the theoretical results, and determined $\chi ^{(4)}=9.44\times 10^{-25}~{\mathrm m^{3}/V^{3}}$ for a 2-nm mesh.

Simulations for Process (III) and (IV)

Process (III) and (IV) are both two-step nonlinear processes to produce a 4HG photons. In process (III), two steps of 2HG produce 215-nm photons. In process (IV), one step of 3HG followed by one step of SFG produce 215-nm photons. We calculated both processes using hydrodynamic Drude model, following the approaches given by Zeng et al. [44].

Maxwell’s equations for hydrodynamic Drude model:

$$\begin{array}{rll} \frac{\partial{\mathbf{B}}}{\partial{t}}&=&-\nabla\times\mathbf{E}\notag\\ \frac{\partial{\mathbf{E}}}{\partial{t}}&=&c^{2}\nabla\times\mathbf{B}-\frac{\mathbf{J}}{\epsilon_{0}}\notag\\ \rho&=&\epsilon_{0}\nabla\cdot\mathbf{E}\notag\\ \frac{\partial{\mathbf{J}}}{\partial{t}}&=&\frac{e^{2}n_{0}}{m}\mathbf{E}-\gamma\mathbf{J}+\sum\limits_{k}\frac{\partial}{\partial r_{k}}\left(\frac{\mathbf{J}J_{k}}{en_{0}-\rho}\right)\notag\\ &&-\frac{e}{m}\big[\rho\mathbf{E}+\mathbf{J}\times\mathbf{B}\big]\notag \end{array} $$

$\omega _{p}^{2}=\frac {e^{2}n_{0}}{m\epsilon _{0}}$ where $\omega _p=1.367\times 10^{16}~{\mathrm s}^{-1}$ is the bulk plasma frequency and the damping frequency $\gamma =6.478\times 10^{13}~{\mathrm s}^{-1}$ [44].

The hydrodynamic Drude model contains nonlinear equations which are linearized by perturbative expansions

$$ \mathbf{E}=\sum\limits_{j}\mathbf{E}^{(j)} $$

(20)

$$ \mathbf{B}=\sum\limits_{j}\mathbf{B}^{(j)}$$

(21)

$$ \mathbf{J}=\sum\limits_{j}\mathbf{J}^{(j)} $$

(22)

$$ \rho=\epsilon_{0}\sum\limits_{j}{\nabla\cdot\mathbf E}^{(j)} $$

(23)

where the superscript (j) indicates the field of the jth order.

The first-order equations are the local Drude model which can be directly simulated in FDTD simulations, given as

$$ \frac{\partial \mathbf{B}^{(1)}}{\partial t}=-\nabla\times\mathbf{E}^{(1)} $$

(24)

$$ \frac{\partial\mathbf{E}^{(1)}}{\partial t}=c^{2}\nabla\times\mathbf{B}^{(1)}-\frac{\mathbf{J}^{(1)}}{\epsilon_{0}} $$

(25)

$$ \frac{\partial{\mathbf{J}^{(1)}}}{\partial{t}}=\frac{e^{2}n_{0}}{m}\mathbf{E}^{(1)}-\gamma\mathbf{J}^{(1)} $$

(26)

The second-order equations are given by

$$ \frac{\partial \mathbf{B}^{(2)}}{\partial t}=-\nabla\times\mathbf{E}^{(2)} $$

(27)

$$ \frac{\partial\mathbf{E}^{(2)}}{\partial t}=c^{2}\nabla\times\mathbf{B}^{(2)}-\frac{\mathbf{J}^{(2)}}{\epsilon_{0}} $$

(28)

$$ \begin{array}{rll} \frac{\partial{\mathbf{J}^{(2)}}}{\partial{t}}&=&\frac{e^{2}n_{0}}{m}\mathbf{E}^{(2)}-\gamma\mathbf{J}^{(2)}+\mathbf{S}^{(2)}\\ \mathbf{S}^{(2)}&=&\sum\limits_{k}\frac{\partial}{\partial r_{k}}\left(\frac{\mathbf{J}^{(1)}J_{k}^{(1)}}{en_{0}}\right)-\frac{e}{m} \end{array}$$

(29)

$$ \times \big[\epsilon_{0}(\nabla\cdot\mathbf{E}^{(1)})\mathbf{E}^{(1)}+\mathbf{J}^{(1)}\times\mathbf{B}^{(1)}\big] $$

(30)

where subscript k represents the x, y, and z coordinates. $\mathbf {S}^{(2)}$ is the nonlinear source term driving the second-order response. From Eq. 30, we calculated the complex electric field dipole source for the 2HG to be

$$\begin{array}{rll} \tilde{\mathbf{E}}_{\mathrm{s}}^{(2)}(2 \omega_{0})&=& {} \frac{\tilde{\mathbf{S}}^{(2)}(2\omega_{0})}{{\epsilon_{0}}{\omega_{p}^{2}}} \notag\\ &=& {} -\frac{\epsilon_{0}}{en_{0}} {} \left\{\left[1+\frac{\omega_{p}^{2}}{(\omega_{0}+i\gamma)^{2}}\right] \left[\left(\nabla\cdot\tilde{\mathbf{E}}^{(1)}\right)\tilde{\mathbf{E}}^{(1)}\right]\right.\notag\\ && {\kern2.5pc} \left.+{} \frac{\omega_{p}^{2}}{2\omega_{0}(\omega_{0}+i\gamma)}\nabla\left(\tilde{\mathbf{E}}^{(1)}\cdot\tilde{\mathbf{E}}^{(1)}\right) {} \right\} \end{array} $$

(31)

In process (III), the first step of 2HG is calculated by launching the time domain electric dipole sources according to $\tilde {\textbf {E}}_{\mathrm s}^{(2)}(2\omega _{0})$ at all computational nodes inside gold.

The fourth-order equations are given by

$$ \frac{\partial \textbf{B}^{(4)}}{\partial t}=-\nabla\times\textbf{E}^{(4)} $$

(32)

$$ \frac{\partial\textbf{E}^{(4)}}{\partial t}=c^{2}\nabla\times\textbf{B}^{(4)}-\frac{\textbf{J}^{(4)}}{\epsilon_{0}} $$

(33)

$$ \frac{\partial{\textbf{J}^{(4)}}}{\partial{t}}=\frac{e^{2}n_{0}}{m}\textbf{E}^{(4)}-\gamma\textbf{J}^{(4)}+\textbf{S}_{1}^{(4)}+\textbf{S}_{2}^{(4)} $$

(34)

$$\begin{array}{rll} \textbf{S}_{1}^{(4)}&=&\sum\limits_{k}\frac{\partial}{\partial r_{k}}\left(\frac{\textbf{J}^{(2)}J_{k}^{(2)}}{en_{0}}\right)\\ &&-\frac{e}{m}\left[\epsilon_{0}\left(\nabla\cdot\textbf{E}^{(2)}\right)\textbf{E}^{(2)} +\textbf{J}^{(2)}\times\textbf{B}^{(2)}\right] \end{array} $$

(35)

$$\begin{array}{rll} \textbf{S}_{2}^{(4)}&=& \sum\limits_{k}\frac{\partial}{\partial r_{k}}\left(\frac{\textbf{J}^{(1)}J_{k}^{(3)}+\textbf{J}^{(3)}J_{k}^{(1)}}{en_{0}}\right)\\ &&-\frac{e}{m}\left[\epsilon_{0}\left(\nabla\cdot\textbf{E}^{(1)}\right)\right.\textbf{E}^{(3)}+\epsilon_{0}\left(\nabla\cdot\textbf{E}^{(3)}\right)\textbf{E}^{(1)}\\ &&\left. \qquad \quad +\textbf{J}^{(1)}\times\textbf{B}^{(3)}+\textbf{J}^{(3)}\times\textbf{B}^{(1)}\right] \end{array} $$

(36)

There are two nonlinear source terms: $\textbf {S}_{1}^{(4)}$ describes 2HG of two 2HG photons in process (III); $\textbf {S}_{2}^{(4)}$ describes the sum frequency generation of one FF photon and one 3HG photon in process (IV). From Eqs. 35 and 36, we calculated the complex electric field dipole source for the two processes, respectively.

$$\begin{array}{rll} \tilde{\textbf{E}}_{\mathrm s1}^{(4)}(4\omega_{0})&=&\frac{\tilde{\textbf{S}}_{1}^{(4)}(4\omega_{0})}{\epsilon_{0}\omega_{p}^{2}}= -\frac{\epsilon_{0}}{en_{0}}\\ &&\left\{\left[1+\frac{\omega_{p}^{2}}{(2\omega_{0}+i\gamma)^{2}}\right]\left[\left(\nabla \cdot\tilde{\textbf{E}}^{(2)}\right)\tilde{\textbf{E}}^{(2)}\right]\right. \\ &&\left. {\kern5pt} +\frac{\omega_{p}^{2}}{4\omega_{0}(2\omega_{0}+i\gamma)}\nabla\left(\tilde{\textbf{E}}^{(2)} {} \cdot\tilde{\textbf{E}}^{(2)}\right)\right\} \end{array} $$

(37)

$$\begin{array}{rll} \tilde{\textbf{E}}_{\textrm s2}^{(4)}(4\omega_{0}) {} &=& {} \frac{\tilde{\textbf{S}}_{2}^{(4)}(4\omega_{0})}{\epsilon_{0}\omega_{p}^{2}} \\ &=& {} - {} \frac{\epsilon_{0}}{en_{0}} {} \left\{ {} \left[1\,+\,\frac{\omega_{p}^{2}}{(\omega_{0}+i\gamma)(3\omega_{0}+i\gamma)}\right] \right. \\ &&\qquad \quad\left[\left(\nabla\cdot\tilde{\textbf{E}}^{(1)}\right)\tilde{\textbf{E}}^{(3)}\,+\,\left(\nabla\cdot\tilde{\textbf{E}}^{(3)}\right)\tilde{\textbf{E}}^{(1)}\right] \\ &&\left. \qquad \quad+\frac{\omega_{p}^{2}}{(\omega_{0}\,+\,i\gamma)(3\omega_{0}\,+\,i\gamma)}\nabla {} \left(\tilde{\textbf{E}}^{(1)} {} \cdot {} \tilde{\textbf{E}}^{(3)}\right) {} \right\}\\ \end{array} $$

(38)

To calculate the second step of 2HG in process (III), time domain electric dipole sources (center frequency $4\omega _{0}$) are launched according to $\tilde {\textbf {E}}_{\mathrm s1}^{(4)}(4\omega _{0})$ given in Eq. 37.

The source terms in Eq. 38 also require third-order field $\tilde {\textbf {E}}^{(3)}$ and we implemented the third-order nonlinear susceptibility $\chi ^{(3)}$ directly into the Lumerical FDTD Solutions “$\chi ^{(2)}/\chi ^{(3)}$” model by adding the nonlinear third-order polarization $\textbf {P}^{(3)}$ to the linear material model

$$ P_{k}^{(3)}=\epsilon_{0}\chi^{(3)}\left(E_{k}^{(1)}\right)^{3}, $$

(39)

where subscript k represents the x, y, and z coordinates, $\chi ^{(3)}=7.56\times 10^{-19}~{\mathrm m}^{2}/{\mathrm V}^{2}$, according to Ref. [44].

To calculate the second step nonlinear process in (IV), time domain electric dipole sources (center frequency $4\omega _{0}$) are launched according to $\tilde {\textbf {E}}_{\mathrm s2}^{(4)}(4\omega _{0})$ given in Eq. 38.

Simulations for Process (V)

For process (V), we need to calculate the intensity profile of the 2HG, $I_{\mathrm 2HG}(\mathbf {r})$. We have applied two models to calculate 2HG: (1) hydrodynamic Drude model, according to Eq. 31; and (2) nonlinear surface $\chi ^{(2)}$ model. The surface $\chi ^{(2)}$ model was implemented in the similar way as the surface $\chi ^{(4)}$ model used in Section “Simulations for Process (II): Direct 4HG.” $\chi ^{(2)}=7.11\times 10^{-9}~{\mathrm m/V}$ was determined by simulating a planar gold surface and comparing the simulated 2HG intensity with the experiments and theories [52, 53]. Figure 8 compares the developed volume calculated from both models. Since the volume is slightly larger in the hydrodynamic Drude model, we followed the hydrodynamic Drude model when comparing the processes in Fig. 3.

Fig. 8

The developed volume around gap antenna calculated from hydrodynamic Drude (HD) model and the surface $\chi ^{(2)}$ model. The curve calculated from the two models are very close to each other. Hydrodynamic Drude model was used in the article

Simulations for Process (VI)

For process (VI), we need to calculate the intensity profile of the fundamental frequency following the method in Section “Simulations for Process (II): Direct 4HG” and the third-order field as described in Section “Simulations for Process (III) and (IV).”