Identification of Plasmonic Modes in Parabolic Cylinder Geometry by Quasi-Separation of Variables

This paper describes the plasmonic modes in the parabolic cylinder geometry as a theoretical complement to the previous paper (J Phys A 42:185401) that considered the modes in the circular paraboloidal geometry. In order to identify the plasmonic modes in the parabolic cylinder geometry, analytic solutions for surface plasmon polaritons are examined by solving the wave equation for the magnetic field in parabolic cylindrical coordinates using quasi-separation of variables in combination with perturbation methods. The examination of the zeroth-order perturbation equations showed that solutions cannot exist for the parabolic metal wedge but can be obtained for the parabolic metal groove as standing wave solutions indicated by the even and odd symmetries.


Introduction
Metallic tapered waveguides have attracted considerable attention as a possible experimental structure for the waveguides in achieving deep subwavelength confinement not only in optical region [1,2] but also in terahertz (THz) region [3][4][5]. Theoretical investigation of the metallic tapered waveguides has been conducted so far in the context of superfocusing beyond the diffraction limit of focused electromagnetic waves [6][7][8]. Superfocusing in metallic tapered waveguides originates from two remarkable theoretical papers [9,10] in 1997 regarding the peculiarities of surface plasmon polaritons (SPPs) [11,12] known as electromagnetic waves propagating along a metal-dielectric interface. The first [9] predicted subwavelength waveguides of SPPs by J. Takahara (one of the present authors) and his colleagues, whereas the second [10] introduced the concept of superfocusing for SPPs at the apexes of metallic tapered waveguides. These papers together provide a complete conceptual description of SPPs in metallic tapered waveguides at both extreme sides. A substantive issue in theoretical physics is to identify the plasmonic modes in metallic tapered waveguides that clarify the intermediate behavior of SPPs propagating from an infinite distance to the apex of tapered waveguides.
Our approach to finding plasmonic modes in metallic tapered waveguides is to solve a wave equation by the quasiseparation of variables (QSOV) [13][14][15] in combination with the use of perturbation methods, which is based on the idea of smoothly connecting analytic solutions at both extreme sides of metallic tapered waveguides. The QSOV technique was first applied to conical geometry [13] and then to wedgeshaped geometry [14], with unprecedented results of plasmonic modes that enable us to consider superfocusing properties with respect to the taper angle in conical-and wedge-shaped geometries. A remaining important topic in metallic tapered waveguides is to investigate the effect of round tips of tapered waveguides, which is very difficult to consider using conventional analytical methods based on classical separation of variables (CSOV) such as the geometrical-optics (GO) approximation (also called adiabatic approximation) [16][17][18][19][20][21][22][23][24][25][26][27] and local separation of variables (LSOV) approximation [10,[28][29][30][31]. In a previous paper [15], we successfully applied the QSOV technique to a circular paraboloidal geometry, which is one of the simplest tapered geometries with a finite tip radius, allowing us to discover that despite the metallic tapered waveguide, superfocusing does not take place in this geometry. Some of the knowledge is practically useful for analysis of THz experimental results in the circular paraboloidal geometry [32]. To properly understand the effect of round tips of metallic tapered waveguides, more theoretical investigation using the QSOV technique is needed into the tapered geometries with finite tip radii.
In this paper, we describe a theoretical investigation of plasmonic modes in a parabolic cylinder geometry by solving the wave equation for the magnetic field by means of the QSOV in combination with perturbation methods. We observed that plasmonic modes cannot exist for the parabolic metal wedge but can be obtained for the parabolic metal groove as standing waves classified by even and odd symmetries. This result is remarkably similar to that in the circular paraboloidal geometry [15], where we previously observed that the plasmonic modes cannot exist for the metallic solid paraboloid but can be obtained for the metallic hollow paraboloid as standing waves. Importantly, the superfocusing does not take place in the parabolic cylinder geometry as well as in the circular paraboloidal geometry [15]; hence, from a parabolical standpoint, the present paper is interrelated and complementary to the previous paper [15] treating the circular paraboloidal geometry. The new knowledge obtained in this paper is a fundamental building block in plasmonics and is practically useful for explaining numerical simulation through metallic subwavelength structures consisting partly of parabolic cylinder geometry, as described in [33][34][35].

Quasi-Separation of Variables Applied to the Wave Equation for a Magnetic Field
We consider two parabolic cylinder structures: a metallic parabolic cylinder surrounded by a dielectric and a dielectric parabolic cylinder surrounded by metal. In this paper, these are simply called a parabolic metal wedge and a parabolic metal groove, respectively. As shown in Fig. 1, the parabolic cylindrical coordinates (ξ,η,z) can be used to describe an infinite parabolic cylinder with permittivity ε 1 , and surrounded by matter with permittivity ε 2 . The parabolic cylindrical coordinates are related to the Cartesian coordinates (x,y,z) according to the following transformations [36][37][38]: where λ 0 is the wavelength of interest in vacuum used for a scale factor. We here choose for the y-coordinate of (1) to take ffiffi ffi ξ p as positive; we must then take ffiffi ffi η p with the positive sign on the semi-parabola above the x-axis and with the negative sign on the lower semi-parabola (see Fig. 1). Curves of constant ffiffi ffi ξ p and ffiffi ffi η p in the cross-section of the x-y plane for (1) are shown in Fig. 2, where all of the parabolas that open in the direction of the negative x-axis are full parabolas to which a single value of ffiffi ffi ξ p within the interval 0 ≤ ffiffi ffi ξ p < ∞ , i.e., 0≤ ξ < ∞, is assigned, while the parabolas that open in the direction of the positive x-axis are regarded as consisting of two semi-parabolas: the upper one is described by þ ffiffi ffi η p and the lower one by − ffiffi ffi η p so that ffiffi ffi η p has the interval −∞ < ffiffi ffi η p < þ ∞ [38]. Therefore, all the points in Fig. 2 are uniquely assigned to both ξ and ffiffi ffi η p , where 0≤ξ<∞ and −∞ < ffiffi ffi η p < þ∞ . For the parabola of Fig. 1, the upper and lower semi-parabolas are Fig. 1 Geometry of a parabolic cylinder structure for the parabolic metal wedge and the parabolic metal groove. ρ 0 is the radius of curvature for the parabolic cylinder structure at the apex. ε 1 and ε 2 are the permittivities inside and outside the parabolic cylinder, respectively. H z is the magnetic field of SPPs with translational symmetry about the z-axis. E ξ and E η are the electric fields of the radial and angular components, respectively, induced by the magnetic filed H z . þ ffiffi ffi η p ¼ ffiffiffiffi ffi η 0 p and − ffiffi ffi η p ¼ ffiffiffiffi ffi η 0 p are the equations for the upper and lower surfaces, respectively, of the parabolic cylinder structure. Geometric dimensions of the xand y-axes are normalized by the wavelength in vacuum, λ 0 In the case where the imaginary parts of the permittivities are ignored, we have the conditions ε 1 <0, ε 2 >0 for the parabolic metal wedge and the conditions ε 1 >0, ε 2 <0 for the parabolic metal groove. Here, we examine simple plasmonic modes with translational symmetry about the z-axis, whose magnetic field at time t at the point located by the coordinate vector x is H(x,t) directed along the z-axis and depend on the intervals ξ∈[0,∞) and ffiffi ffi η p ∈ −∞; ∞ ð Þ ; therefore, it is written as in the parabolic cylindrical coordinate system. Assuming the harmonic time dependence e −iωt , we can write where ω is the angular frequency of interest. Combining the two Maxwell curl equations without any sources, we find the wave equation for the magnetic field in (4): where c is the velocity of light in vacuum. After using the relation λ 0 ω/c=2π, Eq. (5) can be rewritten as [37] ∂ ffiffi with the definition of the magnetic field H z ξ; ffiffi ffi η p À Á in the regions of permittivities ε 1 and ε 2 as H z1 ξ; ffiffi ffi η p À Á and H z2 ξ; ffiffi ffi η p À Á , respectively. In (6), the magnetic field may be classified into even (+) or odd (−) symmetries with respect to ffiffi ffi η p , because the geometrical structure in Fig. 1

is invariant
ffiffi ffi η p ; z→zÞ and the differential operator of the wave Eq. (6) is even under its reflection transformation. Accordingly, we express the magnetic field for even (+) and odd (−) symmetries such as H þ z ξ; ffiffi ffi η p À Á and H − z ξ; ffiffi ffi η p À Á , respectively. As in previous papers [13][14][15], here we look for solutions in the following form by using quasi-separation of variables: with the boundary conditions to determine Η j s (η,ξ) uniquely, where s indicates either the even (+) or odd (−) symmetry with respect to ffiffi ffi η p , i.e., (double-sign corresponds). Substituting (7) into (6), we get for j=1,2 that Since the left side of (10) depends on ξ alone, while the right side depends on both ξ and η, both sides must depend on ξ alone. Setting both sides equal to ζ j (ξ) for j=1,2 and rearranging terms in each equation, we obtain the radial equations and the angular equations Here, ζ j s (ξ) for j=1,2 are called quasi-separation invariants, or simply separation quantities, which are analogous to the separation constants in classical separation of variables (CSOV). Now, we can replace (6) with (11) and (12), which satisfy the boundary conditions of (8).

Unification of the Radial Equations in the Two Regions
The radial equations (11) are separately expressed in two different regions j=1,2. Remarkably, those equations can be transformed into a unified form independent of the regions by considering the limiting cases ξ→0+ and ξ→∞.
in which the material terms ε j are lost and therefore the different expressions in the regions j=1,2 are a superfluity of equations. Then, by setting ζ u in the unified notation for the two different regions. Equation (15) is used to investigate asymptotic solutions for Ξ u s (ξ) as ξ→0+. Two nontrivial linearly independent particular solutions to (15) are given by (4), we see that exp iν ffiffi ffi ξ p ð Þ and exp −iν ffiffi ffi ξ p ð Þ in (16) correspond to the outgoing and incoming waves, respectively (the use of this condition may seem strange for a discussion of the condition ξ→0+), but it is very useful for characterization as seen in "Boundary Conditions for the Radial and Extended Angular Functions" section for ξ→∞ when ν>0.
For ξ→∞, since the physical situation is considered the same as for SPPs in planar geometry, (11) can be described with a unified notation for the two regions j=1,2 in the form Ξ u s (ξ)=Ξ j s (ξ) and this unified radial function satisfies the Sommerfeld radiation conditions [40]: , and k p is the wave number of SPPs in the planar geometry, given [41,42] with the wave number in vacuum, k 0 , defined by k 0 =ω/c. Using (1), we can transform (17) into the form This equation allows us to choose simple candidates for the unified radial function in the form, Ξ u s (ξ)=exp(±iλ 0 k p ξ/2) for ξ→∞, which, unfortunately, do not belong to a class of solutions for (11). Then, carefully choosing for (18) that Ξ u s (ξ)=ξ −1/4 exp(±iλ 0 k p ξ/2) for ξ→∞, we obtain their differential equation under the acceptable condition, ξ→∞. Since (20) can be derived from (11) by using the as the limiting equation of (11) for ξ→∞. Now, we are ready to look for a unified form of (11) by following an inversion process of finding the limiting Eqs. (15) and (20) as ξ→0+ and ξ→∞, respectively. A simple candidate for the unified radial equation is given by which approaches (15) and (20) as ξ→0+ and ξ→∞, respectively. Then, we can find without difficulty that more general candidate is given by where ζ u s (ξ) satisfies the condition ζ u s (ξ)/ξ→0 as ξ→∞. Comparing (23) with (11), we find that the unification conditions are given by which allows for the unified notation of the radial functions as follows: Further discussion on the unified radial Eq. (23) would require more detailed information on the unified separation quantity ζ u s (ξ), which can be determined from the boundary conditions.

Boundary Conditions for the Radial and Extended Angular Functions
In the preceding section, we showed that original radial Eq. (11) can be simplified into the unified radial Eq. (23), in which the unified radial function for ξ→ 0+ can be expressed by a linear combination of the outgoing wave, exp iν ffiffi ffi ξ p ð Þ , and the incoming wave, exp −iν ffiffi ffi ξ p ð Þ , as seen in (16). In this case, the unified radial function Ξ u s (ξ) is divided into an incoming part, Ξ in s (ξ), and an outgoing part, Ξ out s (ξ); that is which satisfies the boundary condition By choosing and substituting (28) into (26), we obtain the boundary condition of the unified radial function as follows: where H 0 is the complex amplitude of the magnetic field to be determined by the initial conditions. In order to obtain the boundary condition of the angular functions, we can use the continuity of the radial electric field at the metal-dielectric surface (þ ffiffi ffi η p ¼ ffiffiffiffi ffi η 0 p ). The Ampére-Maxwell equation in the absence of current density allows us to describe the electric field as a function of the magnetic filed. Since the magnetic field is given by (3), the electric field E(x,t) can be described as in which we see from (4) that Each component of the electric field is described by using the magnetic field as where the superscript j=1,2 indicates the region of permittivities ε 1 and ε 2 , respectively. Substituting (7) into (32), and using (25), we can rewrite the radial electric field for j=1,2 as which is continuous at the metal-dielectric surface (þ ffiffi ffi η p ¼ ffiffiffiffi ffi η 0 p ). Hence, we find the following expression for the boundary conditions of the angular functions: Note that the continuity of the radial electric field on the other metal-dielectric surface (− ffiffi ffi η p ¼ ffiffiffiffi ffi η 0 p ) is automatically satisfied when the magnetic field has either even (+) or odd (−) symmetries with respect to ffiffi ffi η p .

Application of Perturbation Methods for Solving the Extended Angular Equations
According to (24) and (25), the extended angular Eqs. (12) and (13) for j=1,2 can be simplified to the forms In spite of this simplification, we still face difficulty with the fact that (36) and (37) include the as-yet-to-be-determined unified radial function Ξ u s (ξ). Moreover, we still face a basic difficulty in solving such a complicated partial differential equation (PDE) containing first-and second-order partial derivatives with respect to ffiffi ffi ξ p and a second-order partial derivatives with respect to ffiffi ffi η p . To overcome these difficulties, we apply perturbation methods to the extended angular Eq. (36) by treating F s j ffiffi ffi η p ; ξ À Á on the right side as a perturbing term because an exact solution can be found for the left side. According to the perturbation theory [43], we introduce the perturbation parameter 0≤δ≤1 into (36) by replacing F s j ffiffi ffi η p ; ξ À Á with δ F s j ffiffi ffi η p ; ξ À Á on the right side and look for solutions of the perturbed equation of the form Accordingly, ζ u s (ξ) and Ξ u s (ξ) should be described as respectively. For Ξ u (ξ), the incoming Ξ in (ξ) and the outgoing Ξ out (ξ) should also be described as respectively. Substituting (38)-(40) into the perturbed equation, and setting the coefficients of the powers of δ equal to each other, we have a system of equations for Η ,…, in the power series of (38); for coefficient of (38) and using (8), we obtain for j=1,2 that which leads to the following Dirichlet boundary conditions by setting the coefficients of the powers of δ equal to each other in (44). In a similar manner, by substituting (38) into (35), we obtain the following Neumann boundary conditions Note that the imposition of both the Dirichlet (45) and Neumann boundary conditions (46) is referred to as the Cauchy boundary conditions [44]; that is, we specify the values and normal derivatives of the functions Η Perturbation methods should also be applied to the unified radial Eq. (23). Substituting (39) and (40) Substituting (40) into (29), and setting the coefficients of the powers of δ equal to each other, we have the following Dirichlet boundary conditions: In a similar manner, by substituting (40)- (42) into (26), we have the following relations: Similarly, by substituting (41) and (42) into (27) Summing up the points we have discussed in this section, we arrive at a sequence of problems, P 0 , P 1 , P 2 ,…, from which we can find the function 1, 2, … in sequence; the problem P 0 is given by The higher-order problems, P 1 , P 2 ,…, are not shown here because they are not treated in this paper due to some difficulties involved in numerical calculations.

Extended Angular Functions of the Zeroth-Order for the Parabolic Metal Groove
The zeroth-order extended angular Eq. (43) are considered for the parabolic metal groove, in which we can write ε 1 = ε d and ε 2 = ε m by introducing notations of metallic permittivity ε m < 0 and dielectric permittivity ε d >0. The solutions to (43) must be divided into even (+) and odd (−) symmetries with respect to ffiffi ffi η p , de- Substituting ε 1 =ε d (>0) into the zeroth-order extended angular Eq. (43) for j=1, we have . Dividing (52) by 4πε d 1/2 and rearranging terms, we obtain which is a special case of the modified parabolic cylinder Eq. (A.7) corresponding to the parameter values Since the solutions to (53) must be divided into even (+) and odd (−) symmetries with respect to ffiffi ffi η p , written (51), we can describe the solutions to (53) as where w 1 and w 2 are even and odd functions, respectively (see Appendix A).
Characteristic Equations for Determining the Unified Separation Quantity of the Zeroth-Order for the Parabolic Metal Groove Á in (60) into the boundary condition for s=+, ε 1 =ε d , and ε 2 =ε m in (35), we obtain where we use notations of w 1 0 a; z ð Þ ¼ ∂ z w 1 a; z ð Þ and U ' (a,z)=∂ z U(a,z). Equation (61) is the characteristic equation for determining the zeroth-order unified separation quantity for the even (+) symmetry, ζ u where we use notation of w 2 ' (a,z)=∂w 2 (a,z)/∂z. Equation (62)  (ξ), respectively. In this case, we choose to solve them numerically. In numerical calculations, we use the dielectric permittivity ε d =1 and the metallic permittivity ε m = −20 by assuming that the dielectric matter is air with a permittivity of ε d =1 and the metallic matter is gold with a permittivity of ε m =−20.6+1.57i at a wavelength of 750 nm [45]; the imaginary part of ε m is significantly smaller than the real part and thus can be ignored for the sake of simplicity. For the numerical calculations, we used a set of Fortran subroutines provided by Gil et al. [46] for computing U(a,x) to numerically solve characteristic Eqs. (61) and (62) in Mathematica (Wolfram Research, Inc.). Without this set, we were not able to obtain the whole profiles of Figs. 3 and 4 from (61) and (62) by employing the usual methods of defining U(a,x) as in Mathematica. The elementary outline of how to call a Fortran subroutine from Mathematica can be found in [47].  Fig. 4 are fitting curves for the respective solid lines, though they are unclear due to their overlapping with the solid lines. In Fig. 4a, the red and black lines are used for 0.001 ≤ η 0 ≤ 0.1 and 0.2 ≤ η 0 ≤ 1.0, respectively, because the behavior of profiles is clearly expressed. For the curve fitting, taking into account the approximate behavior of ζ u s(0) (ξ), s=+,− for large values of ξ (described in Appendix B), we select a figure-of-merit function with four parameters, p 0 , p 1 , p 2 , and p 3 , expressed as follows: In Tables 1 and 2, we show several values of p 0 , p 1 , p 2 , and p 3 for specific values of η 0 obtained by fitting (63) to the solid curves in Fig. 4a, b, respectively (to maintain clarity in the diagram, some of the curves are not shown).

Unified Radial Functions of the Zeroth-Order for the Parabolic Metal Groove
In the parabolic metal groove, the zeroth-order unified separation quantities ζ u (61) and (62), respectively, are quite accurately estimated by the figureof-merit function in (63). Even with it, we cannot solve the zeroth-order unified radial Eq. (47) as an already-known differential equation. In this subsection, we demonstrate that (47) roughly approximates the modified parabolic cylinder Eq. (A.7) through transformations.
By defining the modified wave numbers of SPPs in the planar geometry, k mp s (ξ) for s=+,−, as from (65). The rough approximation in (67) can be numerically estimated by the following ratio:  Table 1 Values of the parameters of the even (+) symmetry in (63) for the specific value of η 0 (=ρ 0 /λ 0 ) obtained by fitting the solid curves in Fig. 4a η 0 ξ u (0) are numerically calculated in the same manner as that employed for drawing the thick line in Fig. 3

Table 2
Values of the parameters of the odd (−) symmetry in (63) for the specific value of η 0 (=ρ 0 /λ 0 ) obtained by fitting the curves in Fig. 4 . By using (31)- (33) for (80), we obtain a differential equation expressed in terms of the total derivatives with respect to ξ and η as follows: which becomes the exact differential equation dψ(ξ,η)= 0 with the solution ψ(ξ,η)=constant, if we set ψ ξ; η ð Þ ¼ H z ξ; ffiffi ffi η p ; t−π=2ω À Á . Then, for the time-varying scalar field, f(ξ,η,t), of the electric field-line representation, we can set The field-line pattern at time t=t 0 is described by the scalar field with the contour f(ξ,η,t 0 )=C, where C indicates the contour level. The contour interval can be controlled by the interval value of the contour level. The evolution of field-line patterns can be investigated for different moments: t=t 0 +nτ with n=1,2,3,⋅⋅⋅, where τ denotes a suitable chosen duration between two neighboring snapshots.
For actual calculations, using (4), (7), and (25), we rewrite (82) in a form that contains the unified radial function and the extended angular function with even (+) and odd (−) symmetries. The scalar field of the zeroth-order approximation for (82) is given by Figures 6 and 7 show the electric field-line patterns of the zeroth-order approximated plasmonic modes of (a) even and (b) odd symmetries in the parabolic metal groove for the surfaces η 0 =0.3 and η 0 =0.1, respectively, where the surface η 0 is characterized in (2). The electric field-line patterns were obtained by using (83) with t=π/2ω. In the calculations, for the even symmetry, we employed Η  Figure 8 shows the electric field-line patterns of the zerothorder approximated plasmonic modes of odd symmetry in the parabolic metal groove for the surface η 0 =0.204785; this is the specific value of η 0 for ζ u Fig. 6 Electric field lines of the zeroth-order plasmonic mode of the a even and b odd symmetries in the parabolic metal groove for the surface η 0 (=ρ 0 /λ 0 )=0.3, obtained using (83) at t=π/2ω under the conditions ε d =1 and ε m =−20. The blue and red lines indicate the clockwise and counterclockwise loops, respectively, of the line of electric force. The geometric dimensions of the horizontal and vertical axes are normalized by the wavelength in vacuum λ 0 Fig. 7 Electric field lines of the zeroth-order plasmonic mode of the a even and b odd symmetries in the parabolic metal groove for the surface η 0 (=ρ 0 /λ 0 )=0.1, obtained using (83) at t=π/2ω under the conditions ε d =1 and ε m =−20. The blue and red lines indicate the clockwise and counterclockwise loops, respectively, of the line of electric force. The geometric dimensions of the horizontal and vertical axes are normalized by the wavelength in vacuum λ 0 to Dr. F. Peper of NICT for his useful discussions and comments on calling a Fortran subroutine from Mathematica. This work was supported by KAKENHI (21560073).