We present a novel analytical approach for efficient sensitivity analysis of surface plasmon polaritons (SPPs) waveguide-based structures using the beam propagation method (BPM). Our approach exploits the adjoint variable technique to extract the response sensitivities with respect to all the design parameters regardless of their number. No extra BPM simulations are required. The accuracy of the results are comparable to those obtained using the expensive central finite difference approximations applied at the response level. Our approach is successfully applied to different SPPs structures for different applications.
© 2008 Optical Society of America
Some metals such as gold and silver have negative dielectric constant at optical frequncies. This is similar to the behavior of electron gas or cold plasma -. At the interface of these metals with any material with positive dielectric constant the plasma oscillators (plasmons) interact with light to produce plaritons. The light is propagating and guided to the interface as a result of this strong coupling in the form of transverse magnetic (TM) surface wave .
Surface plasmon polaritons (SPPs) based structures attracted much attention in the last decade [1-10]. This is mainly due to their ability to guide light in subwavelength dimensions [1-6]. These structures can thus combine the small physical dimensions of electronic circuits with the high operational speed of optical circuits. These features allow the SPPs circuits to be utilized in subwavelength optical interconnects and photonic circuits. SPPs guide the light at the interface between the dielectric and the metal with exponential decaying in both media. The guided optical field is very sensitive to the properties of the media. Consequently, SPPs have been also successfully exploited in biosensing application [7-9].
Investigating the behavior and the sensitivity of the SPPs structures is essential to determine its potential in the different proposed applications. The sensitivity information of these structures with respect to its design parameters is also useful for yield and tolerance analyses. Various modeling techniques have been utilized for studying the behavior and the sensitivity of the optical field in different SPPs structures. They include method of lines , the finite difference frequency domain (FDFD) , and the finite element (FEM) method . Most of these techniques are utilized for solving the modes of the SPPs structures to investigate the behavior and the sensitivity of the propagation constants and the losses of the different modes. Very little research is carried out to study the behavior of the propagating optical field for these structures . For both applications, all of the utilized techniques are based on performing many simulations to study the sensitivity of the SPPs structures with respect to the design parameters. This approach, however, is inefficient especially for a large number of design parameters. For example, for N design parameters N extra simulations are need to extract the sensitivity information if the forward finite difference approximation is used.
In order to efficiently model the SPPs waveguide structures, the beam propagation method (BPM) can be exploited. It is capable of both calculating the modal field and simulating the propagating optical field assuming that the utilized BPM technique suits the structure under investigation. This technique has been recently applied for both applications for SPPs structures [10,11]. Accordingly, efficient sensitivity extraction of SPPs structures using BPM is of prime importance.
In order to reduce the computational cost of obtaining sensitivity information, the adjoint variable method (AVM) is exploited. This method has been recently utilized to extract the sensitivity of various photonic structures using different numerical methods such as the FDFD  and the finite difference time domain FDTD . It has been also utilized for efficient sensitivity analysis using BPM [14,15]. However, the application of the AVM to the BPM was limited to the lossless case. In order to calculate the sensitivity of any objective function using AVM based on BPM simulation, the derivatives of the system matrices need to be calculated before starting the simulation. The design parameters may include dimensions of discontinuities or material properties. In addition, the utilized approaches in [14,15] are based on exploiting the perturbation theory to estimate a finite difference approximation of these derivatives. The accuracy of this approach is, however, limited by the accuracy of the utilized finite difference approach. In addition, for large structures and large number of design parameters, the computational cost of calculating these derivatives can be high.
In this paper, we propose an analytical approach to extract the sensitivity of any objective function with respect to all the design parameters of SPPs structures. This approach is based on calculating the analytical derivatives of the system matrices with respect to the refractive index. This analytical derivative is then utilized to extract the sensitivity of any type of design parameters. In addition, we utilize the Cauchy-Riemann relationship to reduce the computational cost of obtaining this analytical derivative of the system matrices for the case of SPPs structures.
We start by briefly reviewing the 3D BPM using the alternative direction implicit method (ADI BPM) governing equations in Section 2. Our sensitivity approach and its practical implementation are presented in Section 3. In Section 4, our approach is illustrated through different examples. Finally, the conclusions are given in Section 5.
2. 3D ADI BPM
The governing equations of the semivectorial 3D BPM can be given by 
where Ψ y is the slowly varying envelope of the transverse electric field component E y. Similarly, the BPM equation for the envelope of the field component in the x direction Ψ x can be obtained . In Eq.(1), k is the free space wave number, n is the refractive index, and no is the reference refractive index.
For the TM case we have:
Similarly the operator P xx can be separated to x-dependent and y-dependent components for the TE case. The solution of the system of equations given in Eq. (2) is performed by solving two tridiagonal system of equations at each propagation step .
3. Sensitivity analysis using adjoint variable method (AVM)
In this section, the procedure to extract the sensitivity information using our analytical approach is proposed. This procedure is described mainly for the TM case. It can, however, be extended to the TE in a straight forward way.
The sensitivity of any objective function f with respect to the design parameter pi is given by 
where 1, i=2…,N is the index of the design parameters pi. These design parameters may be dimensions of discontinuities or constitutive parameters. In Eq. (6), the objective function is assumed to depend on the field states at different propagation steps zk i.e. , lk∈I∀k, where I is the index set of the associated field states. By differentiating (2) with respect to the ith design parameter pi, we get
The derivatives of the system matrices in Eq. (8) are given by
The sensitivity of the desired objective function can thus be calculated using the following expression in :
The derivatives of the system matrices required in Eqs. (7)-(9) are estimated using perturbation theory in [14,15]. However, this technique may be inefficient for large number of design parameters and large structures. This is mainly due to the need to construct the system matrices of the perturbed system for each design parameters. By utilizing the explicit dependence of these system matrices on the refractive index, we provide an analytical derivative. This can be done by assuming a change in the refractive index value of Q grid points due to the change in the design parameter pi. Following a similar approach to that in , the derivatives of the system matrices can be rewritten using the following formula
The derivative of the system matrix with respect to the refractive index value of the qth grid point ∂R/∂nq is analytically calculated. This derivative is sparse with non zero elements at the elements corresponding to the grid point q. This derivative ∂R/∂nq can be calculated in terms of ∂Bx/∂nq and ∂By/∂nq which can be calculated analytically. For example, in the case of finite difference discritization of the operators given in Eq.(5), the matrix Bx can be written as
Here Δx is the spatial grid size in the x direction.
By differentiating Bx with respect to the refractive index of grid point q we obtain
where I is the identity matrix. The term ∂By/∂nq can be calculated using a similar approach. It follows that the derivative ∂R/∂nq can be calculated analytically in an efficient manner. It is calculated only once for all the design parameters. This approach reduces the computational cost of obtaining the derivative of the system matrices especially for structures with large number of design parameters.
For SPPs structures, the metal is represented by a complex refractive index. Hereafter, we propose an efficient approach to calculate the derivatives of the system matrix with respect to design parameters which are related to the change of the dimensions or the material parameters of the metal. This complex refractive index can be written as
where nr represents the real refractive index of the material and nim represents the losses of this material. Using the theory of complex analysis, the derivative of the system matrix with respect to the refractive index for the complex case is given by:
where Rkr and Rkm are the real part and the imaginary part of the system matrix Rk. Here, the matrix Rk is assumed to be analytical. In addition, the derivatives of the system matrices with respect to the imaginary part of the refractive index can be efficiently calculated using the Cauchy-Riemann formula as follow :
It follows that the sensitivity with respect to any design parameter can be obtained efficiently by utilizing only one analytical derivative of the system matrix.
4. Numerical examples
In order to verify the accuracy of the implemented 3D ADI BPM, the obtained results are compared with those obtained in  using a similar method. This comparison is based on the imaginary distance version of the 3D ADI BPM for a structure containing a metal loaded on a dielectric substrate as shown in Fig.1. The same parameters in  are utilized to extract the effective index and the attenuation coefficient. An excellent agreement is obtained between the implemented method and the results in . As an example, the effective index of the given structure is calculated using our implemented 3D ADI BPM and compared with the one in  as shown in Fig. 2. The results obtained using our approach are also compared with those obtained using FEMLAB . This commercial software package utilizes the finite element method (FEM) to solve optical structures. The results are also shown in Fig. 2. It is clear from these results that a good agreement is also obtained between our implemented 3D ADI BPM and the FEM.
In the following examples, the implemented 3D ADI BPM is utilized in our proposed AVM approach given in Section 3. In these examples, two different applications of the SPPS waveguide structures are studied using 3D ADI BPM. The sensitivities of these structures are calculated using our AVM approach and compared with the central finite difference CFD applied directly at the level of response.
4.1 A metal loaded on a channel dielectric waveguide
In this example, a SPP waveguide design is proposed for subwavelength applications as shown in Fig. 3. In this design, the silicon on insulator (SOI) material is utilized due to its wide application in electronic circuits. It also allows for strong guiding and hence subwavelength light confinement. The utilized metal is gold (Au) with refractive index nm=0.18-j10.2 at a wavelength of 1.55 µm. The refractive indexes of silicon (Si) and the insulator (SiO2) are ns=3.46, and ni=1.46, respectively. The thickness of the metal layer and the silicon layer are tm=0.10 µm, and ts=0.50 µm, respectively. The width of the metal and silicon layers are taken as Wm=Ws=0.5 µm. The propagation length of this waveguide structure is given as Lp=1/(2Im(β)), where β is the complex propagation constant of the fundamental TM mode of this structure. This length is calculated using the imaginary distance 3D ADI BPM to be 32.0 µm. The vector of design parameters of this structure is p=[ts tm Ws Wm]T. The sensitivity of the propagation length and the effective index of the fundamental mode are calculated using the proposed AVM. The results are compared with the sensitivity information obtained using the CFD approach applied directly at the response level. A very good agreement is obtained between our approach and the CFD as shown in Figs. 4-6. The CFD requires 8 additional simulations. Our approach, however, requires no additional simulations. In general, to obtain an optimal response, the sensitivity should be ideally zero. Thus our target is to minimize the sensitivity to enhance the response. As clear from Fig. 4, the sensitivity of the propagation length is decreasing with the increase of the thickness of the silicon layer ts. Thus increasing the silicon layer thickness reduces the loss of the fundamental mode. On the other hand, Fig. 5 shows that the sensitivity of propagation length approaching zero when the width of the metal layer Wm approaches 0.5 µm which is the same value for the width of the silicon layer Ws. The same behavior can be seen in Fig. 6, where the sensitivity of the propagation length decreases when Ws approaches the nominal value of Wm which is 0.5 µm. Figs. 5 and 6 indicate that in order to reduce the losses and hence increase the propagation length, the width of both the silicon layer and the metal layer should be similar or ideally the same. The width mismatch between the two layers forces the field lines at the edges of the metals to pass through the air before reaching the silicon region. This results in an additional discontinuity. Thus we may conclude that the presence of the additional discontinuity due to width mismatch increases the losses of the whole structure.
4.2 A compact 1×3 power splitter
In this example, a compact 1×3 power splitter is designed using SPPs waveguide structure as shown in Fig. 7. This simple structure consists of a thin layer of metal (Au) with refractive index nm=0.18-j10.2 deposited on an insulator (SiO2) with refractive index ns=1.459. The various widths of the metal sections define the regions of single mode and multimode operation. The design of this structure is based on the multimode interference in the multimode waveguide by utilizing the self imaging principle . The width of the metal in the multimode section Wm is taken to be 15.0 µm. The width of the metal in the single waveguides Ws is taken to be 3.0 µm to ensure single mode operation at 1.55 µm . The thickness of the metal layer of the multimode region tm is taken to be 0.2 µm. The thickness of the metal layer of the single mode region ts is also taken to be 0.2 µm The sensitivity of the coupling coefficient of the power splitter with respect to the design parameters of the structure p=[Wm tm Lm ts]T is studied using our AVM technique. The coupling coefficient is defined as
where Φ N and Ψ are the normalized modal field of the output waveguides and the normalized field distribution at the output waveguides, respectively. The length of the multimode section Lm is calculated to be 78.2 µm. This length, which is calculated using the 3D ADI BPM, maximizes the power coupling coefficient at the output waveguides. A good agreement is obtained between the sensitivity information obtained using our AVM approach and the sensitivity obtained using the time consuming CFD approach. In order to illustrate the accuracy of our approach the normalized sensitivity of the power coupling with respect to width of the multimode region Wm is shown in Fig. 8. The normalized sensitivity of the power coupling sensitivity is also calculated with respect to the thickness of the metal in the multimode region tm as shown in Fig 9. As shown from this figure, the sensitivity of the coupling power approaches zero when the thickness reaches the value of 0.2 µm. Hence, the power coupling is maximum and the loss is minimized once the thickness exceeds this value. This result agrees with the result obtained in .
By utilizing the CFD approach, 8 additional simulations are needed to obtain the sensitivity of all the design parameters. On the other hand, our approach is utilized to obtain the sensitivity information with no additional simulation.
A novel analytical approach is proposed for efficient sensitivity analysis of SPPs structures. The proposed approach is based on exploiting AVM with analytical derivative of the system matrices. This approach is applied for mode calculations as well as propagation problems. The approach is simple and easy to implement. The obtained sensitivities have comparable accuracy to those obtained using the CFD. The universality of the propose approach is successfully illustrated through different examples.
References and links
1. S. A. Maier, Plasmonics: Fundamentals and Applications, (Springe, 2007).
2. P. Berini, “Plasmon-polariton waves guided by thin lossy metal films of finite width: Bound modes of symmetric structures,” Phys. Rev. B, Condens. Matter. 61, 10484–10503 (2000). [CrossRef]
3. G. Veronis and S. Fan, “Modes of Subwavelength Plasmonic Slot Waveguides,” J. Lightwave Technol. 25, 2511–2521 (2007). [CrossRef]
5. I. Breukelaar, R. Charbonneau, and P. Berini, “Long-range surface plasmon-polariton mode cutoff and radiation in embedded strip waveguides,” J. Appl. Phys. 100, 043104-1-043104-9 (2006). [CrossRef]
6. N. N. Feng, M. L. Brongersma, and L. D Negro, “Metal-Dielectric slot-waveguide structures for the propagation of surface plasmon polaritons at 1.55 µm,” IEEE J.Quantum Electron. 43, 479–485 (2007). [CrossRef]
7. J. Homola, S. S. Yee, and G. Gauglitz, “Surface plasmon resonance sensors: review,” Sens. Actuators B 54, 3–15 (1999). [CrossRef]
8. R. D. Harris and J. S. Wilkinson, “Waveguide surface plasmon resonance sensors,” Sens. Actuators B 29, 261–267 (1995). [CrossRef]
10. J. Shibayama, S. Takagi, T. Yamazaki, J. Yamauchi, and H. Nakano, “Numerical analysis of waveguide-based surface Plasmon resonance sensor with adsorbed layer using two- and three-dimensional beam-propagation methods,” IEICE Trans. Electron.E 90-C, 95–100 (2007). [CrossRef]
11. J. Shibayama, S. Takagi, T. Yamazaki, J. Yamauchi, and H. Nakano, “Eigenmode analysis of a light-guiding metal line loaded on a dielectric substrate using the imaginary-distance beam-propagation method,” J. Lightwave Technol. 23, 1533–1539 (2005). [CrossRef]
13. M. A. Swillam, M. H. Bakr, and X. Li, “Accurate sensitivity analysis of photonic devices exploiting the finite-difference time-domain central adjoint variable method,” Appl. Opt. 46, 1492–1499.(2007). [CrossRef] [PubMed]
14. M. A. Swillam, M. H. Bakr, and X. Li, “Efficient adjoint sensitivity analysis exploiting the FD-BPM,” J. Lightwave Technol. 25, 1861–1869 (2007). [CrossRef]
15. M. A. Swillam, M. H. Bakr, and X. Li, “Full vectorial 3D sensitivity analysis and design optimization using BPM,” J. Lightwave Technol. 26, 528–536 (2008). [CrossRef]
16. Y. Hsueh, M. Yang, and H. Chang,“Three-dimensional noniterative full-vectorial beam propagation method based on the alternating direction implicit method,” J. Lightwave Technol. 19, 2389–2397 (1999). [CrossRef]
17. P. A. W. Basl, M. H. Bakr, and N. K. Nikolova, “Efficient estimation of sensitivities in TLM with dielectric discontinuities,” IEEE Microwave Wirel. Compon. Lett. 15, 89–91 (2005). [CrossRef]
18. J. W. Brown and R. V. Churchil, Complex Variables and Applications (McGraw-Hill, 2003).
19. FEMLAB, 2.3 ed. COMSOL AB, Sweden, 2002. http://www.comsol.com
20. L. B. Soldano and E. C. Pennings, “Optical multi-mode interference devices based on self-imaging: principles and applications,” J. Lightwave Technol. 13, 615–627 (1995). [CrossRef]