The theoretical basis for Schottky junction-based SPPs and their behavior when exposed to a high-intensity dressing field is described in this section. The optical response of a dressed metal is analyzed, and expressions for its susceptibility and dielectric functions are derived. The Floquet–Fermi golden rule is used to examine the damping effects in dressed metal electrons, demonstrating the possibility of manipulating the metal dielectric function. Ultimately, expressions for the dispersion relations of potential SPP modes at a lossy Schottky junction are derived.
Dielectric function of dressed metals
This part of the analysis presents the derivation of an expression for the conductivity of a dressed metal film. The equilibrium linear response theory is extended to a driven system out of equilibrium, which is probed by a weak external potential. Moreover, an expression for the current response is derived, which depends not only on the frequency of the probe but also on the frequency of the drive. Finally, an analytical expression for the susceptibility and dielectric functions of a dressed metal is presented.
Floquet theory
We consider a metal film subjected to a high-intensity periodical dressing field. The wave function of a single electron in the dressed metal can be identified as Floquet states22,25,26
$$\begin{aligned} {|{\psi _{\alpha }(t)}\rangle } = \exp \left( -i\frac{\epsilon _{\alpha }}{\hbar } t\right) {|{u_{\alpha }(t)}\rangle }. \end{aligned}$$
(1)
In this case, \(\alpha\) represents a discrete set of quantum numbers, \({|{u_{\alpha }(t)}\rangle }\) are corresponding Floquet modes, and \(\epsilon _{\alpha }\) are the corresponding quasienergies. The periodicity of the Floquet modes, allow to rewrite them as a Fourier expansion as follows
$$\begin{aligned} {|{u_{\alpha }(t)}\rangle } = \sum _{n} \exp (-in\Omega t) {|{u_{\alpha }^n}\rangle }, \end{aligned}$$
(2)
where \(\Omega\) is the angular frequency of the dressing field.
The Floquet picture of the Kubo formula
In linear response theory, the average current response of an equilibrium system for a probe bias with a vector potential \({\textbf {A}}({\textbf {k}},\omega )\) can be expressed by the general Kubo formula in \({\textbf {k}}\)-momentum space27
$$\begin{aligned} \begin{aligned} {\langle {{\hat{J}}^a({\textbf {k}},\omega )}\rangle }{} & {} = \lim _{\mu \rightarrow 0^+} -\frac{1}{2\pi \hbar V} \sum _{b} \int _{-\infty }^{\infty } \int _{-\infty }^{\infty } \int _{-\infty }^{\infty } \bigg \langle { [{\hat{j}}^a({\textbf {k}},t),{\hat{j}}^b(-{\textbf {k}},t’)]\bigg \rangle }_0 A^b({\textbf {k}},t’) \exp (i\omega t) \frac{\exp [-i\omega ‘ (t-t’)]}{\omega ‘ + i\mu } ~ {\textrm{d}}\omega ‘ {\textrm{d}}t'{\textrm{d}}t \\{} & {} – \frac{e^2n_m}{m} A^a({\textbf {k}},\omega ). \end{aligned} \end{aligned}$$
(3)
Here, \({\hat{j}}\) is the single-particle current operator of the system, \(\omega\) is the frequency of the response current, V is the system volume, \(a,b \in \{x,y,z\}\) present the directional components in three-dimensional coordinate space, and \({\langle {\cdot }\rangle }_0\) denotes the statistical thermal average with respect to the system without the probe bias. Due to the presence of high-intensity external dressing in our case, the system will not be in equilibrium. Since the occupation number of Floquet states are independent, we can assume our system to be in a stationary state under the dressing field28. From second quantization formalism, we can expand the current operators \({\hat{j}}^{a,b}({\textbf {k}},t)\) using the Floquet states as the basis for our system as follows29
$$\begin{aligned} {\hat{j}}^a({\textbf {k}},t) = \sum _{\alpha \beta } {\hat{j}}^a_{\alpha \beta }({\textbf {k}},t) {\hat{a}}_{\alpha }^{\dagger }(t_0 = 0){\hat{a}}_{\beta }(t_0 = 0), \end{aligned}$$
(4)
where \({\hat{a}}_{\alpha }^{\dagger }(t)\) and \({\hat{a}}_{\alpha }(t)\) are creation and annihilation operators for \(\alpha\)-th Floquet state. These operators should satisfy the following relationships
$$\begin{aligned} {\hat{a}}_{\alpha }^{\dagger } (t){|{0}\rangle } = {|{\psi _{\alpha }(t)}\rangle }, \quad {\hat{a}}_{\alpha }(t){|{0}\rangle } = 0,\quad [{\hat{a}}_{\alpha }(t),{\hat{a}}_{\beta }^{\dagger }(t)]_{\pm } = \delta _{\alpha \beta }, \quad \text {and} \quad [{\hat{a}}_{\alpha }(t),{\hat{a}}_{\beta }(t)]_{\pm } = [{\hat{a}}_{\alpha }^{\dagger }(t),{\hat{a}}_{\beta }^{\dagger }(t)]_{\pm } = 0. \end{aligned}$$
(5)
Here, \({|{0}\rangle }\) is the vacuum state containing no particle, and positive (negative) subscripts refer to the fermionic anticommutators (commutators) for the Floquet state particles. Furthermore, \({\hat{j}}^a_{\alpha \beta }({\textbf {k}},t)\) is the single particle matrix element given in Schrödinger’s picture by
$$\begin{aligned} {\hat{j}}^a_{\alpha \beta }({\textbf {k}},t) = \langle {\psi _{\alpha }(t)}|{{\hat{j}}^a({\textbf {k}})}|{\psi _{\beta }(t)}\rangle = \langle {u_{\alpha }(t)} |{\exp \left( \frac{i\epsilon _{\alpha }t}{\hbar }\right) {\hat{j}}^a({\textbf {k}}) \exp \left( -\frac{i\epsilon _{\beta }t}{\hbar }\right) }|{u_{\beta }(t)}\rangle , \end{aligned}$$
(6)
where \({\hat{j}}^a({\textbf {k}}) = {\hat{j}}^a({\textbf {k}},t=0)\). Considering the Fourier series expand of the time-periodic Floquet modes, we can identify that
$$\begin{aligned} {\hat{j}}^a_{\alpha \beta }({\textbf {k}},t) = \sum _{n_1,n_2 = -\infty }^{\infty } \exp \left\{ \frac{i}{\hbar }[(\epsilon _{\alpha } – \epsilon _{\beta }) + (n_1 – n_2)\hbar \Omega ]t\right\} \Big \langle {u_{\alpha }^{n_1}}\Big |{{\hat{j}}^a({\textbf {k}})}\Big |{u_{\beta }^{n_2}}\Big \rangle , \end{aligned}$$
(7)
where \(n_1,n_2\) are integers. Next, we can evaluate the commutator operator in equation (3) using the above derived second quantization expansion
$$\begin{aligned}{}[{\hat{j}}^a({\textbf {k}},t),{\hat{j}}^b(-{\textbf {k}},t’)] = \sum _{\alpha \beta } {\hat{j}}^a_{\alpha \beta }({\textbf {k}},t) {\hat{j}}^a_{\beta \alpha }(-{\textbf {k}},t’) ({\hat{a}}_{\alpha }^{\dagger }{\hat{a}}_{\alpha } – {\hat{a}}_{\beta }^{\dagger }{\hat{a}}_{\beta }). \end{aligned}$$
(8)
Here, \({\hat{a}}_{\alpha }^{\dagger }{\hat{a}}_{\alpha }\) represent the number operator for the \(\alpha\)-th Floquet state. Thus, we introduce the distribution functions for the Floquet state particles in the system as follows
$$\begin{aligned} {\mathscr {F}}_{\alpha } {:}{=}{\langle {{\hat{a}}_{\alpha }^{\dagger }{\hat{a}}_{\alpha }}\rangle }, \quad \text {and} \quad {\mathscr {F}}_{\beta } {:}{=}\big \langle {{\hat{a}}_{\beta }^{\dagger }{\hat{a}}_{\beta }\big \rangle }. \end{aligned}$$
(9)
It is important to note that these distribution functions not necessarily to be equilibrium distribution functions, however it is assumed that these are time-independent. Thus, the statistical thermal expectation value of the above commutator becomes
$$\begin{aligned} \begin{aligned} \big \langle { [{\hat{j}}^a({\textbf {k}},t),{\hat{j}}^b(-{\textbf {k}},t’)]\big \rangle }_0&= \sum _{\alpha \beta } \sum _{n_1, \ldots ,n_4 = -\infty }^{\infty } \exp \left\{ \frac{i}{\hbar }[(\epsilon _{\alpha } – \epsilon _{\beta }) + (n_1 – n_2)\hbar \Omega ]t\right\} \exp \left\{ \frac{i}{\hbar }[(\epsilon _{\beta } – \epsilon _{\alpha }) + (n_3 – n_4)\hbar \Omega ]t’\right\} \\&\quad \times \big \langle {u_{\alpha }^{n_1}}\big |{{\hat{j}}^a({\textbf {k}})}\big |{u_{\beta }^{n_2}}\big \rangle \,\, \big \langle {u_{\beta }^{n_3}}\big |{{\hat{j}}^b(-{\textbf {k}})}\big |{u_{\alpha }^{n_4}}\big \rangle ({\mathscr {F}}_{\alpha } – {\mathscr {F}}_{\beta }). \end{aligned} \end{aligned}$$
(10)
Then, we can represent the vector potential corresponding to the probe bias in the frequency domain as follows
$$\begin{aligned} A^b({\textbf {k}},t’) = \frac{1}{2\pi } \int _{-\infty }^{\infty } A^b({\textbf {k}},\omega ”)\exp (-i\omega ”t’) ~{\textrm{d}}\omega ”, \end{aligned}$$
(11)
and submitting these back into the Kubo formula and evaluating the time integrals and \(\omega ‘\) integral, we can obtain
$$\begin{aligned} \begin{aligned} \langle {{\hat{J}}^a({\textbf {k}},\omega )}\rangle&= \lim _{\mu \rightarrow 0^+} -\frac{1}{\hbar V} \sum _{b} \int _{-\infty }^{\infty } \sum _{\alpha \beta } \sum _{n_1, \ldots ,n_4 = -\infty }^{\infty } \frac{\delta \left( \omega – \omega ” + (n_1 – n_2 + n_3 -n_4)\Omega \right) }{\left[ \omega + \frac{1}{\hbar }(\epsilon _{\alpha } – \epsilon _{\beta }) + (n_1 – n_2)\Omega + i\mu \right] }\\&\quad \times \big \langle {u_{\alpha }^{n_1}}\big |{{\hat{j}}^a({\textbf {k}})}\big |{u_{\beta }^{n_2}}\big \rangle \,\,\big \langle {u_{\beta }^{n_3}}\big |{{\hat{j}}^b(-{\textbf {k}})}\big |{u_{\alpha }^{n_4}}\big \rangle A^b({\textbf {k}},\omega ”) \left( {\mathscr {F}}_{\alpha } – {\mathscr {F}}_{\beta }\right) ~{\textrm{d}}\omega ” – \frac{e^2n_m}{m} A^a({\textbf {k}},\omega ). \end{aligned} \end{aligned}$$
(12)
We can
specify the probe electric field relation with vector potential using the electromagnetic theory
$$\begin{aligned} A^b({\textbf {k}},\omega ) = -\frac{i}{(\omega + i\gamma )} E^b({\textbf {k}},\omega ). \end{aligned}$$
(13)
Here, the factor \(\gamma\) is a phenomenological way to include electron scattering-caused damping effects in the quantum calculations. Then, we can rewrite the Kubo formula in a compact form
$$\begin{aligned} \langle {{\hat{J}}^a({\textbf {k}},\omega )}\rangle = \sum _{b} \int _{-\infty }^{\infty } \sigma ^{ab}({\textbf {k}},\omega ,\omega ”) E^b({\textbf {k}},\omega ”) ~{\textrm{d}}\omega ”, \end{aligned}$$
(14)
by initiating the conductivity tensor for the dressed metallic system
$$\begin{aligned} \begin{aligned} \sigma ^{ab}({\textbf {k}},\omega ,\omega ”)&= \lim _{\mu \rightarrow 0^+} \frac{i}{\hbar V} \sum _{\alpha \beta } \sum _{n_1, \ldots ,n_4 = -\infty }^{\infty } \frac{\delta \left( \omega – \omega ” + (n_1 – n_2 + n_3 -n_4)\Omega \right) }{\left[ \omega + \frac{1}{\hbar }(\epsilon _{\alpha } – \epsilon _{\beta }) + (n_1 – n_2)\Omega + i\mu \right] }\\&\quad \times \big \langle {u_{\alpha }^{n_1}}\big |{{\hat{j}}^a({\textbf {k}})}\big |{u_{\beta }^{n_2}}\big \rangle \,\, \big \langle {u_{\beta }^{n_3}}\big |{{\hat{j}}^b(-{\textbf {k}})}\big |{u_{\alpha }^{n_4}}\big \rangle \frac{\left( {\mathscr {F}}_{\alpha } – {\mathscr {F}}_{\beta }\right) }{\omega + i\gamma }+ \frac{ie^2n_m}{m(\omega + i\gamma )} \delta (\omega – \omega ”)\delta _{ab}, \end{aligned} \end{aligned}$$
(15)
Based on this expression, the current is no longer a simple product of conductivity and perturbation electric field as it is convoluted over the bias frequency \(\omega ”\), as opposed to the un-driven case. There is significance in the fact that the conductivity tensor depends on both the response frequency \(\omega\) and the bias frequency \(\omega ”\).
With the above-derived general expression for conductivity, we can determine a new expression for our Floquet system by specifying the conditions which affect the parameters. First, we assume that the response and bias frequency \(\omega\) and \(\omega ”\) are in the central Floquet zone
$$\begin{aligned} |\omega |,|\omega ”|< \left| {\Omega }/{2}\right| \quad \Rightarrow \quad |\omega -\omega ”| < \Omega . \end{aligned}$$
(16)
Under this condition, we can identify that the delta distribution in Eq. (15) only can be non-zero with the condition
$$\begin{aligned} n_1 – n_2 + n_3 – n_4 = 0. \end{aligned}$$
(17)
Focusing on more special case of conductivity by analyzing only the longitudinal conductivity where \(a = b = x\), and assuming that the current response is spatially homogeneous (\({\textbf {k}} \rightarrow 0\)), we can obtain
$$\begin{aligned} \begin{aligned} \sigma ^{xx}(\omega )&= \lim _{\mu \rightarrow 0^+} \frac{i}{\hbar V(\omega + i\gamma )} \frac{e^2}{m^2} \sum _{\alpha \beta } \sum _{n_1, \ldots ,n_4 = -\infty }^{\infty } \frac{\left( {\mathscr {F}}_{\alpha } – {\mathscr {F}}_{\beta }\right) }{\left[ \omega + \frac{1}{\hbar }(\epsilon _{\alpha } – \epsilon _{\beta }) + (n_1 – n_2)\Omega + i0^+\right] } \big \langle {u_{\alpha }^{n_1}}\big |{{\hat{p}}^x}\big |{u_{\beta }^{n_2}}\big \rangle \,\, \big \langle {u_{\beta }^{n_3}}\big |{{\hat{p}}^x}\big |{u_{\alpha }^{n_4}}\big \rangle \\&\quad + \frac{ie^2n_m}{m(\omega + i\gamma )}, \end{aligned} \end{aligned}$$
(18)
where \({\hat{p}}^x\) is the x-directional momentum operator. Polarization and current responses are physically indistinguishable effects of the conductivity tensor. Accordingly, susceptibility function of the system can be expressed in general as a function of linear conductivity30
$$\begin{aligned} \chi ^{xx} (\omega ) = \frac{\sigma ^{xx}(\omega )}{-i\omega }. \end{aligned}$$
(19)
Subsequently, we can identify the susceptibility function of the dressed quantum system as
$$\begin{aligned} \begin{aligned} \chi ^{xx} (\omega )&= \lim _{\mu \rightarrow 0^+} \frac{1}{\hbar \omega (\omega + i\gamma )} \frac{e^2}{m^2V} \sum _{\alpha \beta } \sum _{n_1, \ldots ,n_4 = -\infty }^{\infty } \frac{\left( {\mathscr {F}}_{\alpha } – {\mathscr {F}}_{\beta }\right) }{\left[ \omega + \frac{1}{\hbar }(\epsilon _{\alpha } – \epsilon _{\beta }) + (n_1 – n_2)\Omega + i\mu \right] } \big \langle {u_{\alpha }^{n_1}}\big |{{\hat{p}}^x}\big |{u_{\beta }^{n_2}}\big \rangle \,\, \big \langle {u_{\beta }^{n_3}}\big |{{\hat{p}}^x}\big |{u_{\alpha }^{n_4}}\big \rangle \\&\quad + \frac{e^2n_m}{m\omega (\omega + i\gamma )}. \end{aligned} \end{aligned}$$
(20)
The study involves a dressed metallic system that utilizes the free electron model to explain its transport properties. The particle distribution functions can be described using the Fermi–Dirac distribution. Additionally, very low-temperature conditions (\(T \rightarrow 0\)) are considered, leading to a re-writing of the Fermi–Dirac distribution
$$\begin{aligned} {\mathscr {F}}(\epsilon ) = \lim _{T \rightarrow 0} \frac{1}{\exp (\frac{\epsilon – \epsilon _F}{k_B T})+1} \approx \Theta (\epsilon _F – \epsilon ), \end{aligned}$$
(21)
where, \(\epsilon _F\) is the Fermi energy, \(k_B\) is the Boltzmann constant, T is the absolute temperature, and \({\Theta}(\cdot )\) is the Heaviside function. Now, we can restructure the general susceptibility function according to the dressed metallic system
$$\begin{aligned} \begin{aligned} \chi ^{xx} (\omega )&= \lim _{\mu \rightarrow 0^+} \frac{1}{\hbar \omega (\omega + i\gamma )} \frac{e^2}{m^2V} \sum _{\alpha \beta } \sum _{n_1, \ldots ,n_4 = -\infty }^{\infty } \frac{\left( \Theta (\epsilon _F – \epsilon _{\alpha })- \Theta (\epsilon _F – \epsilon _{\beta })\right) }{\left[ \omega + \frac{1}{\hbar }(\epsilon _{\alpha } – \epsilon _{\beta }) + (n_1 – n_2)\Omega + i\mu \right] } \big \langle {u_{\alpha }^{n_1}}\big |{{\hat{p}}^x}\big |{u_{\beta }^{n_2}}\big \rangle \,\, \big \langle {u_{\beta }^{n_3}}\big |{{\hat{p}}^x}\big |{u_{\alpha }^{n_4}}\big \rangle \\&\quad + \frac{e^2n_m}{m\omega (\omega + i\gamma )}. \end{aligned} \end{aligned}$$
(22)
Considering the transport properties of metallic electrons, it’s only necessary to analyze the behavior of the conduction electrons. This limits the analysis to electrons with energy similar to the Fermi energy. This leads to assume that \(\epsilon _{\alpha } = \epsilon ‘\), and \(\epsilon _{\beta } = \epsilon ”\), where \(\epsilon ‘,\epsilon ” \rightarrow \epsilon _F\). Next, we can apply these changes back into the susceptibility function and obtain
$$\begin{aligned} \begin{aligned} \chi ^{xx} (\omega )&= \lim _{\mu \rightarrow 0^+} \lim _{\epsilon ‘,\epsilon ” \rightarrow \epsilon _F} \frac{1}{\hbar \omega (\omega + i\gamma )} \frac{e^2}{m^2V} \sum _{\alpha \beta } \sum _{n_1, \ldots ,n_4 = -\infty }^{\infty } \frac{\left( \Theta (\epsilon _F – \epsilon ‘)- \Theta (\epsilon _F – \epsilon ”)\right) }{\left[ \omega + \frac{1}{\hbar }(\epsilon ‘ – \epsilon ”) + (n_1 – n_2)\Omega + i\mu \right] } \big \langle {u_{\alpha }^{n_1}}\big |{{\hat{p}}^x}\big |{u_{\beta }^{n_2}}\big \rangle \,\, \big \langle {u_{\beta }^{n_3}}\big |{{\hat{p}}^x}\big |{u_{\alpha }^{n_4}}\big \rangle \\&\quad + \frac{e^2n_m}{m\omega (\omega + i\gamma )}. \end{aligned} \end{aligned}$$
(23)
In our study, we assume that \(\omega \ne 0\) and \(\omega < |{\Omega }/{2}|\). It follows that the denominator of the above expression cannot be very small. Therefore, we can see that the first term of the above expression goes to zero. Finally, now we can assume that there is no contribution from the first term in the susceptibility function. Moreover, we can derive final expression for the susceptibility function for a dressed metallic quantum system as
$$\begin{aligned} { \chi ^{xx} (\omega ) = \frac{\epsilon _0\omega _{pm}^2}{\omega (\omega + i\gamma )}, \quad \text {where} \quad \omega _{pm} = \sqrt{\frac{e^2n_m}{\epsilon _0 m}}. } \end{aligned}$$
(24)
Here, \(\gamma _0\) is the un-driven damping factor. Finally, we can identify the dressed metal dielectric function
$$\begin{aligned} { \varepsilon _m (\omega ) = \varepsilon _{hm} – \frac{\chi ^{xx} (\omega )}{\epsilon _0} = \varepsilon _{hm} – \frac{\omega _{pm}^2}{\omega (\omega + i\gamma )}, } \end{aligned}$$
(25)
where \(\varepsilon _{hm}\) is the high-frequency permittivity of the metal. Although the derived expression is the same as the general Drude-Sommerfeld model description, we need to consider the effects on the damping factor \(\gamma\) induced by the dressing field. As the previous literature22,31 describes, the dressing field modifies the electron wave functions. As long as the electron scattering rate and the electron transport damping factor depend on the wave function, we can control the damping
factor as well as the optical properties of the system by applying an external dressing field. These effects can be analytically evaluated by applying the Floquet–Fermi golden rule for a dressed metallic system22. This enables us to achieve high-performance SPP modes in Schottky junction-based waveguides. Next, we can rewrite the x-directional metal dielectric function by introducing the normalized damping factor \({\tilde{\gamma }}\) as follows
$$\begin{aligned} \varepsilon _m (\omega ) = \varepsilon _{hm} – \frac{\omega _{pm}^2}{\omega (\omega + i{\tilde{\gamma }}\gamma _0)}, \quad \text {where} \quad {\tilde{\gamma }} = \frac{\gamma }{\gamma _0}. \end{aligned}$$
(26)
Here, \(\gamma _0\) is the un-driven damping factor. Moreover, \({\tilde{\gamma }}\) depends on the intensity and the frequency of the applied dressing field22. For details on the derivation of these results, see Section 1 and 2 of the Supplementary Information. Under the results section, we investigate numerically how metal dielectric function can be manipulated under various driving fields.
SPP modes at the Schottky junction
Our primary interest is to analyze the optical properties of Schottky junctions and achieve improved SPP modes at the interface. The interface between a semi-infinite metal and a semi-infinite n-type semiconductor film in the xy plane of the Cartesian coordinate space is considered as shown in Fig. 1a top section. The contact between metals and semiconductors leads to the flow of free charge carriers until their Fermi levels align according to thermodynamic principles, causing changes in charge density and creating a space charge region near the interface of the two materials. Notably, the width of the space charge region d can be manipulated by applying an external voltage across the interface. Due to this, the Schottky junction plays a crucial role in modern electronics.
In general, a complex charge density profile can be approximated by a piecewise-linear fluctuation. Therefore, to model the charge profile of the semiconductor region, a piecewise-linear variation is used as shown in Fig. 1a bottom section. This model has been used in earlier literature17,19, and is also employed in this study for comparison purposes. Additionally, the frequencies of the electromagnetic fields are carefully selected to be in an off-resonant regime to minimize photon absorption by the system.
The complex dielectric functions for n-type semiconductor \(\varepsilon _s\) can be described using the following form32:
$$\begin{aligned} \varepsilon _{s} (\omega ) = \varepsilon _{\text {intra}} (\omega ) + \varepsilon _{\text {inter}} (\omega ) \end{aligned}$$
(27)
The proposed approach explicitly distinguishes between the intraband effects \(\varepsilon _{\text {intra}}\) and the interband effects \(\varepsilon _{\text {inter}}\). In this context, the intraband effects are generated by the free electrons, whereas the interband effects are generated by the bound electrons. The intraband component of the dielectric function is can be characterized by the Drude model33
$$\begin{aligned} \varepsilon _{\text {intra}} (\omega ) = \varepsilon _{\infty } – \frac{\omega _{p}^2}{\omega ^2+i\omega \varsigma _s}. \end{aligned}$$
(28)
Here, \(\varepsilon _{\infty }\) is the high-frequency permittivity of the semiconductor, \(\omega _{p}\) is the plasma frequency of the semiconductor, and \(\varsigma _{s}\) is the damping factor of the semiconductor induced by the electron scattering in the material. The interband component of the dielectric function can be represented by both a real part and an imaginary part:\(\varepsilon _{\text {inter}} (\omega ) = \varepsilon ‘_{\text {inter}} (\omega ) + i\varepsilon ”_{\text {inter}} (\omega )\). Under the absence of interband absorptions we can neglect the imaginary part of the interband component. This leads to
$$\begin{aligned} \varepsilon _{s} (\omega ) = \varepsilon _{hs} – \frac{\omega _{p}^2}{\omega ^2+i\omega \varsigma _s}, \end{aligned}$$
(29)
where \(\varepsilon _{hs} = \varepsilon _{\infty } + \varepsilon ‘_{\text {inter}} (\omega )\). Generally, \(\varepsilon _{hs}\) exhibits dependence on the angular frequency. However, for certain spectral ranges, it can be approximately treated as constant34. Since the semiconductor plasma frequency is contingent upon the electron density, the dielectric function is also influenced by variations in electron density. In our system, the electron density of the semiconductor varies with the coordinates along the z-direction. Hence, the dielectric function of the semiconductor relies on both the angular frequency of the SPP excitation radiation and the spatial coordinate in the z-direction. Therefore, we can determine the dielectric function of the semiconductor in our Schottky junction as follows35
$$\begin{aligned} \varepsilon _{s} (z,\omega ) = \varepsilon _{hs} \left[ 1 – \frac{\omega _{ps}^2(z)}{\omega ^2+i\omega \varsigma _s} \right] , \quad \text {where} \quad \omega _{ps}(z) = \sqrt{\frac{n(z)e^2}{\varepsilon _{hs} \varepsilon _{0} m}}. \end{aligned}$$
(30)
Here, n(z) is position-dependent free charge carrier density, \(\varepsilon _0\) is the vacuum permittivity, and m is the effective mass of the free charge carriers. In addition, we assume that the damping factor is not dependent on the free charge carrier density. Moreover, we define the plasma frequency of the bulk semiconductor region as
$$\begin{aligned} \omega _{ps0} = \sqrt{ \frac{n_se^2}{\varepsilon _{hs}\varepsilon _0 m}}, \end{aligned}$$
(31)
and it is independent of the position. Since the charge density inside the semiconductor of the Schottky junction varies, we can identify four main regions inside a Schottky junction, as depicted in Fig. 1b. For the metal region, we use the general dielectric function expression derived in Eq. (26). With the dressing field removed, one is left with the well-known Drude model-based metal dielectric function expression.Then, we summarize these characteristics of each region under the Table 1.
Next, we try to identify electromagnetic wave modes \({\textbf {E}}(x,y,z,t)\) and \({\textbf {H}}(x,y,z,t)\) that can exist on the Schottky junction. Here, we represent the electric field of the SPP mode with \({\textbf {E}}\) and the magnetic field with \({\textbf {H}}\). Then, we hope to find solutions with two properties: (i) wave modes propagate through the surface towards the x-direction, and (ii) wave modes decay through both mediums in a perpendicular direction to the surface. Without loss of generality, we can assume that \(|{\textbf {E}}|\) and \(|{\textbf {H}}|\) are independent of y-coordinates. Therefore, we can find the solutions for \({\textbf {E}}\) and \({\textbf {H}}\) by defining them as vector components of the Cartesian coordinate system in the following form
$$\begin{aligned} {\textbf {E}} (x,y,z,t) = \left( {\mathscr {E}}_x(z)e^{ik_x x – i\omega t },0,{\mathscr {E}}_z(z)e^{ik_x x – i\omega t }\right) ^{\textsf{T}}, \quad \text {and} \quad {\textbf {H}} (x,y,z,t) = \left( 0,{\mathscr {H}}_y(z)e^{ik_x x – i\omega t },0\right) ^{\textsf{T}}. \end{aligned}$$
(32)
Here, we restricted our study to the TM-polarization mode. Furthermore, \(\omega\) is the frequency of the SPP mode, and \(k_x\) represents the x-directional wavenumber. Using the well-known Maxwell’s equations, we can derive two differential equations for the possible electromagnetic modes in a Schottky junction
$$\begin{aligned} {\mathscr {E}}_x(z) = \frac{i}{k_x} \left[ \frac{\partial {\mathscr {E}}_z(z)}{\partial z} + {\mathscr {E}}_z(z) \frac{\partial \ln {\varepsilon (z,\omega )}}{\partial z} \right] , \end{aligned}$$
(33)
and
$$\begin{aligned} \begin{aligned} \frac{\partial ^2{\mathscr {E}}_z(z)}{{\partial z}^2} + \frac{\partial \ln {\varepsilon (z,\omega )}}{\partial z} \frac{\partial {\mathscr {E}}_z(z)}{\partial z} – \left[ k_x^2 – \frac{\omega ^2 \varepsilon (z,\omega )}{c^2} – \frac{1}{\varepsilon (z,\omega )} \frac{\partial ^2\varepsilon (z,\omega )}{\partial z} + \left( \frac{1}{\varepsilon (z,\omega )} \frac{\partial \varepsilon (z,\omega )}{\partial z} \right) ^2 \right] {\mathscr {E}}_z(z) = 0, \end{aligned} \end{aligned}$$
(34)
where c is the speed of light in vacuum. Applying these two differential equations, we can find the wave mode solutions for each region of the Schottky junction.
Region A
In this region \(\varepsilon (z,\omega ) = \varepsilon _m(\omega )\) is not depending on the z-directional space coordinates, and we can identify that
$$\begin{aligned} \frac{\partial ^2{\mathscr {E}}_z(z)}{{\partial z}^2} – \kappa _A^2 {\mathscr {E}}_z(z) = 0, \quad \text {where} \quad \kappa _A^2 = \left[ k_x^2 – \frac{\omega ^2 \varepsilon _m(\omega )}{c^2} \right] . \end{aligned}$$
(35)
Then we can solve this equation and obtain the solutions
$$\begin{aligned} {\mathscr {E}}_z(z) = {\mathscr {E}}_A e^{\kappa _A z}, \quad \text {and} \quad {\mathscr {E}}_x(z) = \frac{i\kappa _A}{k_x} {\mathscr {E}}_A e^{\kappa _A z} \end{aligned}$$
(36)
Here, we only considered the \({\text{Re}} [\kappa _A] > 0\) solutions (with \(z<0\) values) as we need only the decaying solution on the outside of the interface, and \({\mathscr {E}}_A\) is an unknown real constant.
Region B
In this region, \(\varepsilon (z,\omega ) = \varepsilon _{hs}\) is a real-valued constant. Thus,
$$\begin{aligned} \frac{\partial ^2{\mathscr {E}}_z(z)}{{\partial z}^2} – \kappa _B^2 {\mathscr {E}}_z(z) = 0, \quad \text {where} \quad \kappa _B^2 = \left[ k_x^2 – \frac{\omega ^2 \varepsilon _{hs}}{c^2} \right] . \end{aligned}$$
(37)
For \(\kappa _B^2 > 0\), we can express the solutions as
$$\begin{aligned} {\mathscr {E}}_z(z) = {\mathscr {E}}_{B1} e^{\kappa _B z} + {\mathscr {E}}_{B2} e^{-\kappa _B z} \quad \text {and} \quad {\mathscr {E}}_x(z) = \frac{i \kappa _B}{k_x} \left[ {\mathscr {E}}_{B1} e^{\kappa _B z} – {\mathscr {E}}_{B2} e^{-\kappa _B z} \right] , \end{aligned}$$
(38)
where \({\mathscr {E}}_{B1}\) and \({\mathscr {E}}_{B2}\) are unknown constants. However, under the scenario \(\kappa _B^2 < 0\), we need to introduce different type of solutions such as
$$\begin{aligned} {\mathscr {E}}_z(z) = {\mathscr {E}}_{B1} \cos \left( \sqrt{-\kappa _B^2} z\right) + {\mathscr {E}}_{B2} \sin \left( \sqrt{-\kappa _B^2} z\right) , \quad \text {and} \quad {\mathscr {E}}_x(z) = \frac{i \sqrt{-\kappa _B^2}}{k_x} \left[ -{\mathscr {E}}_{B1} \sin \left( \sqrt{-\kappa _B^2} z\right) + {\mathscr {E}}_{B2} \cos \left( \sqrt{-\kappa _B^2} z\right) \right] . \end{aligned}$$
(39)
Region C
In this region, charge carrier density varies linearly. Thus, the permittivity function of this region is also dependent on the z-directional coordinates. Before further analysis, we should consider the practical aspect of these parameters. Our analysis uses an excitation field frequency \(\omega\) in \(\sim 10^{15} {{\hbox {s}}^{-1}}\) order. In addition, we model the semiconductor using the parameters of n-type doped gallium arsenide material, and we can identify that its damping factor is in \(\sim 10^{12} {{\hbox {s}}^{-1}}\) order. Under these practical conditions, we can assume that \(\varsigma _{s} \ll \omega\), and neglect the effects of damping in the semiconductor region. It is important to notice that for \(d_1 < z \le d_2\), we can assume that
$$\begin{aligned} k_x^2 \ll \left[ \frac{1}{(z-\nu )^2} – \frac{\omega ^2 \varrho }{c^2} (z – \nu )\right] , \quad \text {where} \quad \varrho = – \varepsilon _{hs} \left[ \frac{r\omega _{ps0}^2}{\omega ^2} \right] , \quad \text {and} \quad \nu = – \frac{\varepsilon _{hs}}{\varrho } \left[ 1 + \frac{h\omega _{ps0}^2}{\omega ^2} \right] . \end{aligned}$$
(40)
This leads to
$$\begin{aligned} \frac{\partial ^2{\mathscr {E}}_z(z)}{{\partial z}^2} + \frac{1}{(z-\nu )} \frac{\partial {\mathscr {E}}_z(z)}{\partial z} – \left[ – \frac{\omega ^2 \varrho }{c^2} (z – \nu ) + \frac{1}{(z-\nu )^2} \right] {\mathscr {E}}_z(z) = 0, \end{aligned}$$
(41)
and we can present the solutions as
$$\begin{aligned} {\mathscr {E}}_z(\zeta ) = {\mathscr {E}}_{C1} \frac{{\text {Ai}}'(\zeta )}{\zeta } + {\mathscr {E}}_{C2} \frac{{\text {Bi}}'(\zeta )}{\zeta }, \quad \text {and} \quad {\mathscr {E}}_x(z) = \frac{i}{k_x} \left[ \frac{{\mathscr {E}}_{C1}}{\vartheta } \text {Ai}(\zeta ) + \frac{{\mathscr {E}}_{C2}}{\vartheta } \text {Bi}(\zeta ) \right] . \end{aligned}$$
(42)
Here, \(\zeta = {(z -\nu )}/{\vartheta }\), \(\vartheta = \root 3 \of {{c^2}\big /{\varepsilon _{hs} r\omega _{ps0}^2}}\), \({\text {Ai}}'(\zeta )\) is first derivate with respect to \(\zeta\) of the first kind of Airy function \({\text {Ai}}(\zeta )\), \({\text {Bi}}'(\zeta )\) is first derivate with respect to \(\zeta\) of the second kind of Airy function \({\text {Bi}}(\zeta )\), and \({\mathscr {E}}_{C1}\) and \({\mathscr {E}}_{C2}\) are unknown constants.
Region D
Since the dielectric function in this region is a contact respect to the z-directional space coordinates, we can observe that
$$\begin{aligned} \frac{\partial ^2{\mathscr {E}}_z(z)}{{\partial z}^2} – \kappa _D^2 {\mathscr {E}}_z(z) = 0, \quad \text {and} \quad \kappa _D^2 = \left[ k_x^2 – \frac{\omega ^2 \varepsilon _s(\omega )}{c^2} \right] . \end{aligned}$$
(43)
Same as the region A, we can present the solution as
$$\begin{aligned} {\mathscr {E}}_z(z) = {\mathscr {E}}_D e^{-\kappa _D z} \quad \text {and} \quad {\mathscr {E}}_x(z) = -\frac{i\kappa _D}{k_x} {\mathscr {E}}_D e^{-\kappa _D z} \end{aligned}$$
(44)
Here, we only considered the \({\text{Re}} [\kappa _D] > 0\) solutions (with \(z > d_2\) values) as we need only the decaying solution on the outside of the interface, and \({\mathscr {E}}_D\) is an unknown real constant.
Dispersion relation of SPP modes at the Schottky junction
The dispersion relation for possible SPP modes relates the excitation light’s angular frequency to the SPP mode’s in-plane wavenumber magnitude. The dispersion relation can be found by applying the self-consistent boundary conditions required by Maxwell’s equations36
$$\begin{aligned} \varepsilon _{\text {R}1} {\mathscr {E}}_{z,\text {R}1}(z) = \varepsilon _{\text {R}2} {\mathscr {E}}_{z,\text {R}2}(z), \quad \text {and} \quad {\mathscr {E}}_{x,\text {R}1}(z) = {\mathscr {E}}_{x,\text {R}2}(z), \end{aligned}$$
(45)
where \(\text {R}1,\text {R}2 \in \{\text {A, B, C, D}\}\) at \(z=0\), \(z=d_1\), and \(z=d_2\) to ensure the continuity of the SPP electromagnetic field. Using the previously identified wave modes for each region in the Schottky junction, we can derive six linear simultaneous equations for these boundary conditions, and represent them in a matrix equation
$$\begin{aligned} {M}_{6\times 6}(\omega , k_x) {{\mathscr {E}}}_{6\times 1} = 0. \end{aligned}$$
(46)
Here, \({M}_{6\times 6}(\omega , k_x)\) is the matrix of coefficients, and it can be represented as
$$\begin{aligned} {M}_{6\times 6}(\omega , k_x) = \begin{bmatrix} \varepsilon _m(\omega ) &{} -\varepsilon _{hs} &{} -\varepsilon _{hs} &{} 0 &{} 0 &{} 0 \\ \kappa _A &{} -\kappa _B &{} \kappa _B &{} 0 &{} 0 &{} 0 \\ 0 &{} e^{\kappa _B d_1} &{} e^{-\kappa _B d_1} &{} – \frac{{\text {Ai}}'(\zeta _{d_1})}{\zeta _{d_1}} &{} – \frac{{\text {Bi}}'(\zeta _{d_1})}{\zeta _{d_1}} &{} 0 \\ 0 &{} \kappa _B e^{\kappa _B d_1} &{} – \kappa _Be^{-\kappa _B d_1} &{} -\frac{{\text {Ai}}(\zeta _{d_1})}{\vartheta } &{} -\frac{{\text {Bi}}(\zeta _{d_1})}{\vartheta } &{} 0 \\ 0 &{} 0 &{} 0 &{} \frac{{\text {Ai}}'(\zeta _{d_2})}{\zeta _{d_2}} &{} \frac{{\text {Bi}}'(\zeta _{d_2})}{\zeta _{d_2}} &{} -e^{-\kappa _D d_2} \\ 0 &{} 0 &{} 0 &{} \frac{{\text {Ai}}(\zeta _{d_2})}{\vartheta } &{} \frac{{\text {Bi}}(\zeta _{d_2})}{\vartheta } &{} \kappa _D e^{-\kappa _D d_2} \\ \end{bmatrix}, \end{aligned}$$
(47)
for \(\kappa _B^2 > 0\). However, under the \(\kappa _B^2 < 0\) condition, we can find that
$$\begin{aligned} {M}_{6\times 6}(\omega , k_x) = \begin{bmatrix} \varepsilon _m(\omega ) &{} -\varepsilon _{hs} &{} 0 &{} 0 &{} 0 &{} 0 \\ \kappa _A &{} 0 &{} -\kappa _{B’} &{} 0 &{} 0 &{} 0 \\ 0 &{} \cos (\kappa _{B’} d_1) &{} \sin (\kappa _{B’} d_1) &{} – \frac{{\text {Ai}}'(\zeta _{d_1})}{\zeta _{d_1}} &{} – \frac{{\text {Bi}}'(\zeta _{d_1})}{\zeta _{d_1}} &{} 0 \\ 0 &{} -\kappa _{B’}\sin (\kappa _{B’} d_1) &{} \kappa _{B’}\cos (\kappa _{B’} d_1) &{} -\frac{{\text {Ai}}(\zeta _{d_1})}{\vartheta } &{} -\frac{{\text {Bi}}(\zeta _{d_1})}{\vartheta } &{} 0 \\ 0 &{} 0 &{} 0 &{} \frac{{\text {Ai}}'(\zeta _{d_2})}{\zeta _{d_2}} &{} \frac{{\text {Bi}}'(\zeta _{d_2})}{\zeta _{d_2}} &{} -e^{-\kappa _D d_2} \\ 0 &{} 0 &{} 0 &{} \frac{{\text {Ai}}(\zeta _{d_2})}{\vartheta } &{} \frac{{\text {Bi}}(\zeta _{d_2})}{\vartheta } &{} \kappa _D e^{-\kappa _D d_2} \\ \end{bmatrix}, \end{aligned}$$
(48)
where \(\kappa _{B’} = \sqrt{-\kappa _B^2}\), and \({{\mathscr {E}}}_{6\times 1} = [{\mathscr {E}}_{A} , {\mathscr {E}}_{B1} , {\mathscr {E}}_{B2} , {\mathscr {E}}_{C1} , {\mathscr {E}}_{C2} , {\mathscr {E}}_{D} ]^{{\textsf{T}}}\). A nontrivial solution to \({{\mathscr {E}}}_{6\times 1}\) can be obtained by ensuring that the determinant of the coefficient matrix \({M}_{6\times 6}(\omega , k_x)\) is equal to zero
$$\begin{aligned} \det ({M}_{6\times 6}(\omega , k_x) ) = 0, \end{aligned}$$
(49)
and this leads to a secular equation relating \(\omega\) and \(k_x\). This is the dispersion relation for the possible SPP modes in a Schottky junction.