Open Access
Add to BibSonomy
Abstract
A stepwise angular spectrum method (SASM) for curved interfaces is presented to calculate the wave propagation in planar lens-like integrated optical structures based on photonic slab waveguides. The method is derived and illustrated for an effective 2D setup first and then for 3D slab waveguide lenses. We employ slab waveguides of different thicknesses connected by curved surfaces to realize a lens-like structure. To simulate the wave propagation in 3D including reflection and scattering losses, the stepwise angular spectrum method is combined with full vectorial finite element computations for subproblems with lower complexity. Our SASM results show excellent agreement with rigorous numerical simulations of the full structures with a substantially lower computational effort and can be utilized for the simulation-based design and optimization of complex and large scale setups.
© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Calculating the diffraction of waves in optical systems has been intensively investigated by means of different theoretical approaches and numerical methods [1–16]. A central challenge is the calculation of light propagation through apertures of various forms, i.e. the propagation of the optical field from an input to an output plane or the propagation of light through a lens as an essential part of any imaging system. Almost all of these concepts have their origin from the scalar diffraction theory, which was mainly influenced by the Huygens-Fresnel principle and Kirchhoff’s and Rayleigh-Sommerfeld’s diffraction theory [17–19]. These concepts are some of the oldest and most general formulations to describe the propagation of light from a plane interface to a point of interest behind that surface. The methods can easily be extended to more complex structures, e.g. diffraction at shifted or tilted planes [6–11] or light propagation from/to curved surfaces [12–16].
In this work we want to investigate the diffraction of a slab waveguide lens, as illustrated in Fig. 1. The structure consists of a slab waveguide with an additional lens-shaped part of different, but constant thickness on top of the structure. The incoming wave is semi-guided by the slab waveguide on the input side and limited in lateral direction. The shape of the lens determines the focal point of the outgoing beam. When viewed on the wavelength scale relevant for numerical wave simulations, these configurations must be considered huge, potentially leading to a (prohibitively) large effort of rigorous numerical simulations. Hence here, we present a numerical method to calculate the wave propagation in the structure of Fig. 1, using a variation of the angular spectrum method.
Fig. 1. Convex (a) and concave (b) integrated-optical lens consisting of slab waveguides with different thicknesses connected by curved surfaces. The incoming/outgoing beam is semi-guided by the slab waveguide.
The angular spectrum method (ASM) is a rigorous technique for modeling the diffraction of a wave in homogeneous media by expanding the field into a series of plane and evanescent waves. The mathematical background lies in the field of Fourier Optics [17], but has been further investigated in the area of integrated optics. In [16], a stepwise discretized variant of the ASM was used to determine the diffraction from/to a curved surface, where the curved interface is discretized into a finite number of small, flat elements. In this paper, we propose an extension of this method, which makes calculations of full vectorial 3D light propagation in slab waveguide lenses feasible.
Existing techniques to simulate wave propagation in optical waveguides include, for example, the effective index method, projecting the problem to 2D, or variants of the beam propagation method, relying on the slowly varying envelope assumption. For small radiation losses and low refractive index contrasts, this offers a good way to approximately describe the field behavior. However, if the refractive index contrast is comparably large, pronounced radiation losses and back reflections occur, leading to inconsistent results. Here, we focus on high refractive index contrast Si/SiO$_2$ lens-configurations with non-negligible radiation losses and back reflections.
Earlier concepts [20–25] already discussed integrated optical lenses, of e.g. mode-index, geodesic or grating types, which were mostly investigated experimentally. The ASM method introduced in this work should, in principle, be able to handle a wide range of step-index structures, like those presented in [24,25]. However, those results are not suitable for a quantitative assessment of our method. Instead, we gain confidence by comparisons with rigorous direct numerical simulation methods.
Mainly for illustration of our approach, but in principle applicable as well within an effective-index approach in low-index contrast setups, we start our analysis by considering simple 2D configurations consisting of two different materials separated by a curved surface. To calculate the wave propagation for such a system, we present the theory of the adapted stepwise angular spectrum method (SASM) in Section 2. Additionally, in Section 3 the method is extended to two curved interfaces in 2D, resembling a simple lens, to simulate plane wave propagation through these surfaces. To validate our approach, some examples are given and compared to corresponding 2D finite element method (FEM) solutions using COMSOL Multiphysics [26]. Furthermore, to analyze more realistic structures, the whole method is transferred to 3D, meaning configurations of slab waveguides of two distinct thicknesses, connected by curved surfaces as illustrated in Fig. 1. The theory is explained in Section 4. A series of examples show the validity of the approach by again comparing the results with full vectorial 3D FEM solutions.
2. Single curved interface in 2D
Initially, a 2D interface is considered as illustrated in Fig. 2. Equivalent to a 3D setup with translational invariance in $y$, it consists of two different materials, $n_1$ and $n_2$, separated by an almost arbitrary curved interface (it must be ensured that at this interface no multiple reflections occur) given by the function $g(x)$. The incoming plane wave, laterally restricted in the $x$-direction, propagates in positive $x$-$z$-direction and reaches the surface at an angle of incidence $\varphi _0$ with respect to the $z$-axis. We assume time harmonic fields $\sim \mathrm {exp}(\mathrm {i}\omega t)$, with frequency $\omega =kc=2\pi c/\lambda$ for vacuum wavelength $\lambda$, wavenumber $k$ and speed of light $c$. Furthermore, the whole structure is infinitely extended in $y$-direction and the fields are also constant along $y$ such that we have to consider electromagnetic fields
in the frequency domain.To calculate the outgoing transmitted and reflected fields for a given incoming wave, the surface is discretized in $n$ segments equidistant along the transverse ($x$-axis) direction of length $\Delta x$ with center $x_n$ [16]. Each segment is approximated by a piecewise flat surface $g_n(x)$, tangential to the interface at $x_n$. The red dashed line, shown in Fig. 2, illustrates the approximation of the curved interface. The SASM is then used to calculate the reflected and transmitted outgoing fields from each element. Roughly the method can be outlined as follows: in a first step the field at the approximated interface is Fourier transformed to obtain its angular spectrum. In a second step a modified back-transformation leads to the desired outgoing fields. Adding up the intermediate results for each element leads to the total field solution. Note that the piecewise composite function to approximate the curved interface is, most likely, discontinuous, but this does not have a strong negative influence on the field behavior, since on the one hand the smooth scattered field from all individual elements are superimposed for the overall solution, and on the other hand we chose the number of elements $n$ large enough to observe convergence. Furthermore, the curved elements are less discontinuous than the staircase approximation used in the original method [16], where already good results are achieved.
Fig. 2. 2D sketch of a curved interface $g(x)$ with two different materials $n_1,n_2$. The incoming plane wave reaches the surface at an incidence angle $\varphi _0$. The dashed red line represents the discretization of the curved interface.
By using oblique approximation elements, each solution is calculated in a new local coordinate system $(\xi _n,\eta _n)$ rotated by the inclination angle $\alpha$ of the approximated cell. For each element this angle is given by $\mathrm {tan}\alpha =g'(x_n)$, where $g'(x)=\frac {\mathrm {d}}{\mathrm {dx}}g(x)$ is the spatial derivative of $g(x)$. The local, rotated and shifted coordinate system is then specified by
2.1 Theory
The incoming wave bundle reaches the interface at an oblique primary angle of incidence $\varphi _0$ with respect to the $z$-axis. We assume that it is analytically given as a superposition of plane waves propagating in various directions $(k_x,k_z)$ with amplitude $w_g(k_x;\varphi _0)$ as
To calculate the angular spectrum, the incoming field has to be determined in the respective local coordinate system $(x',z'; k'_x,k'_z)$. Separately, for each element, the incoming field at the interface is then given by
For the known incoming field on each of the $n$ elements in local coordinates $\Psi _{\mathrm {in},n}(x')$, the next step is to calculate its angular spectrum by the inverse Fourier transformation
The radiated fields (in forward or backward direction), generated from each of the $n$ segments, are calculated by propagating the fields from each element using the angular spectrum. Hence, the outgoing field from the $n$-th element in global coordinates $(x,z)$ is given by the Fourier transformation
2.2 Examples for a single interface
We now discuss some examples using the 2D SASM as shown in Fig. 3. Here, we assume a simple circular curvature function $g(x)=R-\sqrt {R^{2}-x^{2}}$, with curvature radius $R>0$ and material parameters $n_1=1$ and $n_2=3$. The curved interface is divided into $n=50$ discretizing elements (we chose this value for all following simulation examples). Such a configuration focuses the incident light and thus represents a simple lens. The incoming wave is a TE-mode with principal field component $\Psi =E_y$, which reaches the curved surface at an incidence angle of $\varphi _0$. The Gaussian beam width is set to $W_c=10\mu m$ and the structure is operating at telecommunication wavelength $\lambda =1.55\mu m$. The integrals (here and in the following) required to calculate the solutions were evaluated using numerical quadrature.
Fig. 3. Field plots of the absolute electric field component $|E_y|$ for different curvature radii $R$ and incidence angles $\varphi _0$. Shown are the results using the SASMo, SASMs and the corresponding COMSOL (C) solution. Additionally, line plots for fixed $z$-coordinate at focal point and fixed $x$-coordinate are shown. The contour lines mark the levels of the absolute field at 2%, 5% and 10%.
For the 2D case, we compare the SASM with oblique elements with the approach from [16] where a staircase approximation for $g(x)$ was used. To indicate the appropriate method, we will use the abbreviations SASMo (SASM with oblique elements) and SASMs (SASM with staircase elements) in the following.
Furthermore, a full 2D COMSOL simulation (FEM) will serve as a reference solution. In COMSOL, we use port boundary conditions to excite a plane wave of Gaussian shape with specific width $W_c=10\,\mu \textrm {m}$ and incidence angle $\varphi _0$. The mesh is given by a maximum element size of $\lambda /10$. Additionally, we use perfectly matched layers to simulate an infinite structure.
Figure 3 shows the absolute field value of the electric field component $|E_y|$ for different curvature radii $R$ and incidence angles $\varphi _0$ and views along one fixed axis ($x$- or $z$-axis fixed). The SASM solutions with oblique (SASMo) and perpendicular (SASMs) elements and the corresponding COMSOL solution (C) are compared. Additional contour lines are added at 2%, 5% and 10% of the absolute field to make the results more comparable. If necessary, the individual contour lines can be identified on the basis of the declining field behavior. Generally, both angular spectrum methods show very good agreement with the COMSOL results. However, using straight elements the amplitude in the focal point is not exactly the same as in the reference FEM solution, especially for small curvature radii, which can be best observed in the line plots. In contrast, when comparing the line plots using the SASMo and the scale of these figures, the solutions match perfectly with the reference solutions. Only very minor differences can be observed on the left of the focal point and also in the contour lines of the transmitted field. It can be stated that both approximations lead to a good agreement with the COMSOL solution, but the method with oblique elements works slightly better for smaller curvature radii and approximates the curved surface more precisely.
Using geometrical optics, the focus distance $f$ of the outgoing beam for this structure can be predicted by the Lensmaker’s equation [18] for a single interface
The focus of the curved interface with material parameters $n_1=1, n_2=3$ and radius $R=20\,\mu \textrm {m}$ or $R=40\,\mu \textrm {m}$ is calculated as $f=30\mu m$ or $f=60\mu m$, respectively. The field plots in Fig. 3 confirm this, even if the focal point is slightly distorted due to spherical aberration.3. Multiple reflections between two curved surfaces
The more interesting case is a configuration with two curved boundaries, e.g. a concave or convex lens as illustrated in Fig. 4. Now, multiple reflections take place between the interfaces, resulting in a superposition of forward and backward propagating waves.
Fig. 4. Configuration with two curved interfaces given by $g^{(1)}(x)$ and $g^{(2)}(x)$. The incoming wave reaches the first surface under an incidence angle $\varphi _0$ and material parameters are $n_1$ in the outer region and $n_2$ between the two interfaces. The red dashed lines represent the approximation of the curved interfaces.
The solution of the electromagnetic field is again calculated by using the SASM with oblique elements from Section 2.1, but now in an adapted cascaded way. This means that the method has to be extended for multiple reflection: after calculating the scattered fields from the first boundary, the transmitted field serves as the new incoming field on the second boundary. Using the SASM for the second interface again, the calculated reflected field serves as the new incoming field on the first interface. This process has to be repeated several times in order to achieve a converged solution.
3.1 Theory
We discuss the 2D case again, meaning a structure that is infinitely expanded in the $y$-direction as considered before in Section 2.1. But now a configuration with materials $n_1$, $n_2$ and two curved interfaces is considered as illustrated in Fig. 4. The boundaries are defined by $g^{(1)}(x)$ and $g^{(2)}(x)$, respectively. Again, the incoming wave reaches the first interface at an oblique incidence angle $\varphi _0$. It should be noted that the method can analogously be applied to a configuration with three different materials.
Discretizing the two interfaces is done in a similar manner as before: the interfaces are divided into $n$ (first interface) or $m$ (second interface) equidistant, tangential, piecewise flat surfaces $g_n^{(1)}(x)$ and $g_m^{(2)}(x)$ with center $x_n$ or $x_m$ and width $\Delta x_n$ or $\Delta x_m$, respectively. Note, that it is not necessary to use the same discretization for both surfaces. $\alpha _1$ refers to the angle of inclination for the first interface and $\alpha _2$ refers to the second. Furthermore, we introduce two different local coordinate systems, $(x',z')$ for the first and $(\hat {x},\hat {z})$ for the second boundary, defined by Eq. (2) with the parameters accordingly adjusted for the considered interface. The relation between the two local coordinate systems is further given by
Additionally, the superimposed angular spectrum on the second interface consists only of the angular spectra of the scattered waves from the first surface propagating in positive $z$-direction
The individual angular spectra on the elements of the second interface are3.2 Examples for convex and concave lenses
To show the correctness of the cascaded approach, we compare the SASMo results with 2D COMSOL solutions. Again, we assume a circular curvature function for both interfaces, that either generate a convex or concave focusing lens. The interface functions are of the form
where $z_{p}$ moves the circle along the $z$-axis to position the interface. The equation has to be adjusted, if either a convex ($R_1>0$, $R_2<0$, $g^{(1)}$ with $-$, $g^{(2)}$ with $+$) or concave lens is considered ($R_1<0$, $R_2>0$, $g^{(1)}$ with $+$, $g^{(2)}$ with $-$). To focus the incoming beam, the refractive index between the interfaces has to be higher than in the outside region for a convex lens. Hence, we choose $n_1=1$ and $n_2=3$ for such a lens. Additionally, it is exactly the other way around for a concave lens, which is why we choose $n_1=3$ and $n_2=1$ for that case. Furthermore, the Gaussian beam width is again given by $W_c=10\mu m$ and the operating wavelength by $\lambda =1.55\,\mu \textrm {m}$. The incoming wave is the TE-mode with principal field component $\Psi =E_y$, which reaches the curved surface at an incidence angle of $\varphi _0$.To keep the simulation time as low as possible, we calculated only a finite number $\nu < \infty$ of back reflections between the two interfaces.The iterative process is repeated until the maximum value of the field amplitude on the second interface (Eq. (28)) is below 0.5% of the initial input field. In numbers, this means that about 6-10 iterations were performed, depending on the properties of the considered structure.
Figure 5 shows some results for a convex (a-c) and concave (d) lens for different curvature radii and normal incidence angle $\varphi _0=0^{\circ }$. The upper row shows the calculated results using the SASMo and the bottom row the corresponding 2D COMSOL solutions. For all configurations, the results agree well; even the contour lines at low level show the same qualitative behavior. Figure 6 compares results for several convex lenses with oblique incidence angles $\varphi _0 > 0$ for different curvature radii $R_1,R_2$ (upper row). Again, we observe good agreement with the corresponding COMSOL solutions (bottom row).
Fig. 5. Field plots of the absolute electric field component $|E_y|$ for convex (a-c) and concave (d) lens configurations with different curvature radii $R_1,R_2$ and fixed incidence angle $\varphi _0=0^{\circ }$. Shown are the results using the SASMo (upper row) and the corresponding COMSOL (C) solution (bottom row). The contour lines mark the levels of the absolute field at 2%, 5% and 10%.
Fig. 6. Field plots of the absolute electric field component $|E_y|$ for convex lens configurations with different curvature radii $R_1,R_2$ and incidence angle $\varphi _0>0^{\circ }$. Shown are the results using the SASMo (upper row) and the corresponding COMSOL (C) solution (bottom row). The contour lines mark the levels of the absolute field at 2%, 5% and 10%.
The focus $f$ of the outgoing beam can be calculated by the Lensmaker’s equation for a thick lens [18] with thickness $d=|z_2-z_1|$ by
4. Curved surfaces in 3D: the slab waveguide lens
We now extend our 2D stepwise angular spectrum method with oblique elements to a practically more relevant 3D setting. Full 3D simulations are often computationally expensive and simulation time as well as memory requirements are high. We show that by using the SASM in combination with vectorical cross-sectional 2D field solutions, the 3D solution can be computed with substantially less effort when compared to full 3D FEM simulations using COMSOL.
4.1 Theory
The curved interface from Section 2.1 is again considered (Fig. 7(a)), but now without the infinite extension in $y$-direction. Instead, slab waveguides of different thicknesses for the incoming and outgoing part are assumed (Fig. 1 and Fig. 7(b)). The fields are now semi-guided waves in the corresponding slabs and as such localized in the core. Thus, the fields are $y$-dependent and we need to consider fully vectorial electromagnetic fields
Fig. 7. a) Curved interface between two different slab regions; b) cross section of an interface between two slab waveguides of thicknesses $d_1$ and $d_2$ with refractive index $n_1$ in the core and $n_2,n_3$ in the claddings.
The incoming waveguide is assumed to be of thickness $d_1$ and the outgoing waveguide of thickness $d_2$. Both waveguides consist of media with refractive indices $n_1$ in the core and $n_2,n_3<n_1$ in the claddings. The curved interface is again given by an almost arbitrary function $g(x)$. The incoming wave, guided by the corresponding slab waveguide, reaches the interface at a primary in-plane incidence angle $\varphi _0$.
To calculate the 3D solution, the SASM is used in combination with full vectorial 2D field solutions of the cross section in the ($y,z$)-plane for different incidence angles $\varphi$. More precisely, the previously presented method described for the 2D structure in Section 2.1 is applied to the configuration of Fig. 7 with some modifications. Generally, the SASM is defined as a superposition of plane waves with specific amplitudes and different directions. These planar waves are now replaced by fully vectorial 2D solutions of the cross section (cf. the discussion of oblique propagation of semi-guided waves in [27–35]). The advantage of using the vectorial 2D solutions is that out-of-plane losses are fully taken into account. For calculating the constituting fields the structure is assumed to be infinite in $x$-direction, i.e. the fields are assumed to be constant along $x$, and the incoming semi-guided mode excites the configuration at oblique angles of incidence. In our earlier works we referred to these solutions as 2.5D. The FEM method implemented in COMSOL is used to calculate these fields. In principle, any other rigorous method that yields accurate vectorial fields for 2.5D wave propagation problems in slab waveguide structures, should be suitable as well, such as the RCWA [36] or eigenmode expansion techniques [28,37]. Since, the structure and fields are no longer $y$-independent, the subsequent SASM steps are performed with the effective refractive indices $n_{\mathrm {eff}_1},n_{\mathrm {eff}_2}$ of the corresponding slab waveguides (these replace the actual refractive indices of the media) instead of $n_1,n_2$, modeling the propagation of the (potentially several) guided modes supported by the slabs. Then, by applying the modified SASM to the 2D solutions the whole 3D field is determined, where the $y$-dependence is now given by the COMSOL solution.
Since only the plane waves are replaced with COMSOL field solutions, the ansatz for the incoming field on the curved surface $\Psi _{\mathrm {in},n}(x')$ and the angular spectrum $A_n(k'_x)$ can be adopted from Eq. (6) and Eq. (7), replacing $n_1$ by $n_{\mathrm {eff_1}}$ and $n_2$ by $n_{\mathrm {eff_2}}$. The outgoing field for each of the $n$ elements is then given by the SASM in combination with the 2D COMSOL solutions for different incidence angles $\varphi$:
Similarly, the method can also be applied to structures with two interfaces by applying the theory from Section 3.1. In contrast to the 2D case, some modifications have to be taken into account: First, due to vectorial 3D solutions, TE and TM modes are excited simultaneously at the interfaces, leading to diffracted fields for each polarization. Therefore, two angular spectra have to be considered for each interface. Second, the constituting 2D COMSOL solution is still calculated for each interface separately for either TE or TM incidence and then assembles to generate the full field. Calculating the constituting solutions including two interfaces, requires a significantly larger number of 2D fields since in addition to the angle of incidence $\varphi$, the distance $d$ between the interfaces has to be varied depending on the position $x$. Therefore, we preferred the way of appropriately assembling the single interface solutions. The drawback is that for small distances $d$ between the two curved surfaces, radiated or evanescent waves resulting from one interface do not lead to further reflections/transmissions on the other interface. Meaning, that only guided modes generate the diffracted fields from each surface. But we shall see that this does not have a huge influence on the field solution.
4.1.1 Expansion of the 2D COMSOL solution
For some particular structures of interest, the focal point may be far from the origin of the curvature. Depending on the interface and the material parameters, large computational windows are required to include the focal point. This would result in an inconveniently high computational effort for the many 2D COMSOL solutions used for the SASM. Hence, it is advantageous to simplify the approach further. This can be done by extracting only the guided 1D slab modes from the COMSOL solution and propagate them analytically in the respective slab waveguides with the corresponding amplitude and wavenumber. To achieve this, the fields of the full vectorial 2D solution are extracted at two cross-section lines along the $y$-direction, at arbitrary positions $z<0$, for the fields in the incoming slab and a position $z>0$, for the outgoing fields, when the discontinuity is positioned in $z=0$ (Fig. 7(b)).
Generally, in both slab waveguides (incoming and outgoing) the fields can be written as a superposition of all guided slab waveguide modes $\boldsymbol {\Psi }_m$ with specific amplitudes $a_m$, wavenumbers $k_{z,m}$ and scattered field $\boldsymbol {\Psi }_{\mathrm {sc}}$
4.1.2 Excitation by laterally guided rib waveguide modes
Instead of using a Gaussian wave packet $w_g(k_x)$, as presented in Eq. (5) for the 2D case to excite the structure with a laterally ($x$-) restricted mode, we now make use of a rib waveguide, which is placed in front of the structure. Such a waveguide, of comparably large width and low etch depth, guides a laterally wide mode that is limited in $x$-$y$-direction. When connecting the rib waveguide to a slab, it can be expected that a wave, which is guided by the rib waveguide and limited in lateral direction, will remain reasonably confined after entering the slab waveguide without spreading much in the adjacent slab. The mode profile of the rib waveguide is calculated with COMSOL’s mode analysis solver and projected onto the infinitely expanded slab waveguide structure. Exciting the 3D structure by a rib waveguide seems to be more realistic, since this type of excitation is potentially easier to implement experimentally than to generate an incident semi-guided Gaussian beam, and is therefore briefly discussed in this section. The theory of exciting a slab waveguide with a laterally restricted mode by using a rib waveguide was already introduced in [32]. Equation (4) determines the incoming field as a superposition of slab waves at different angles with specific amplitudes. The weight function $w_g(k_x)$ is now replaced by an unknown weight function $w_0(k_x)$ – determined by the TE$_0$ mode of the rib waveguide –, thus the incoming field is given by
4.2 Examples: curved slab waveguide discontinuity
To present some examples for a 3D curved interface, two slabs with thicknesses $d_1=0.05\mu m$ and $d_2=0.22\mu m$ and material parameters $n_1=3.45$, $n_2=n_3=1.45$ for silicon dioxide in the cladding and silicon in the core are considered. Both slab waveguides support the fundamental TE and TM modes with different effective refractive indices. The interface function $g(x)$ is again given by $g(x)=R-\sqrt {R^{2}-x^{2}}$ with curvature radius $R>0$. The incoming wave is the fundamental TE$_0$-mode of the incoming rib waveguide of width $15\mu m$, height $0.05\mu m$ (same height as the adjoining incoming slab waveguide) and etch depth $0.02\mu m$.
Figure 8 shows the calculated field solutions for the absolute electric field value $|E_x|$ for curvature radii $R=20\mu m$ (a) and $R=30\mu m$ (b) and normal incidence angle (top view (upper row) and side view (bottom row)). Here, the full vectorial 2D COMSOL solutions are calculated in the range $z\in [-5,20]\mu m$ (with discontinuity in the point $z=0$; hence the 2D solution has to be moved and rotated depending on the position of $g(x)$ and the inclination angle $\alpha$) and for incidence angles from $\varphi =-85^{\circ }$ to $\varphi =85^{\circ }$ with steps of $0.5^{\circ }$.
Fig. 8. Field plots of the absolute electric field component $|E_x|$ for curvature radius $R=20\,\mu \textrm {m}$ (a) and $R=30\mu m$ (b) and normal incidence angle $\varphi =0^{\circ }$; top view (upper row) and side view (bottom row). The corresponding full vectorial 3D COMSOL solution is shown on the right side (C). The contour lines mark the levels of the absolute field at 2%, 5% and 10%.
To show the validity of the 3D SASM approach, we compare the results with full vectorial 3D FEM solutions calculated with COMSOL. In COMSOL, we also use an etched rib waveguide, with the same parameters as mentioned before, to excite the structure with a laterally restricted mode that is guided by the slab waveguide. The incoming field – the mode of the rib waveguide – is calculated by using port boundary conditions and perfectly matched layers are again used to simulate open boundaries. For the single interface solutions, the volume to be calculated is large, since the focus point is quite far away from the boundary. Hence, we simulated the 3D solution only in the range from $z\in [-5,20]\mu m$ and extended the 3D solution with the method from Section 4.1.2. Scattering losses are almost fully decayed at this point, hence the extension is a good approximation of the full 3D solution. The SASM results show good agreement with the 3D FEM solutions (C) for both cases, side and top view. Examples for oblique incidence angles $\varphi _0>0$ are shown in Fig. 9 for different curvature radii $R$. Note that the incoming rib waveguide, rotated by the appropriate incidence angle, is not explicitly shown in the field plots. Only the field in the adjacent slab waveguides is presented after the structure got excited by the rib waveguide mode.
Fig. 9. Field plots of the absolute electric field component $|E_x|$ for different curvature radii $R$ and incidence angles $\varphi _0>0^{\circ }$; top view (upper row) and side view (bottom row). The contour lines mark the levels of the absolute field at 2%, 5% and 10%.
Already for a single interface in 3D, comparably high losses occur. To get an impression of the magnitude of loss, Fig. 10 shows the relative power transmittance T, the reflection R and the total outgoing power to guided modes P$_{\mathrm {out}}$ depending on the incidence angle $\varphi$ for the constituting 2D solutions. It can be stated from the plots that losses can not be neglected, especially for small incidence angles and incoming TM mode. However, for incoming TE mode, losses are fully suppressed beyond a specific higher angle due to the oblique propagation of semi-guided waves as discussed in detail in some of our previous works [27–34].
Fig. 10. Reflection R and transmittance T for fundamental guided TE-/TM-polarization versus the incidence angle $\varphi$ of the 2D constitutive solutions for either incoming TE (a) or TM (b) mode.
4.3 Examples: convex slab waveguide lenses
To complete the work, some examples of slab waveguide structures with two interfaces – representing a lens –, analogous to the Section 3.2, follow next. The waveguide parameters are the same as before, i.e. the slab waveguide thicknesses are $d_1=0.05\mu m$ for the incoming and outgoing slab and $d_2=0.22\mu m$ for the middle part. Material parameters are chosen for silicon $n_1=3.45$ and silicon dioxide $n_2=1.45$. The incoming wave is guided by an adjoining rib waveguide of width $15\mu m$, height $0.05\mu m$ and etch depth $0.02\mu m$. The constituting 2D COMSOL solutions are again calculated in the range $z\in [-5,20]\mu m$ and extended with the method from Section 4.1.1. For different radii $R_1,R_2$ and angle of incidence $\varphi _0$ the electromagnetic field solution is calculated.
To show the validity of the 3D SASM approach, we again compare the results with full vectorial 3D FEM solutions calculated with the FEM in COMSOL. The program was executed on a SMP node with 32 cores and 1TB RAM. The calculation time of a single 3D simulation took about 20 hours and the entire storage space was required. Nevertheless, the 3D problem does not fully converge, because the mesh could not be sufficiently refined. Contrarily, the calculation of the SASM solution on a computer with 16 cores and 32GB RAM took approximately 10 minutes. It should be noted that the 3D FEM computes the full 3D field, while our SASM evaluates the field only on a slice or on cut planes as required for the plots.
Figure 11 and Fig. 12 show the top view (upper row) and side view (bottom row) of the absolute electric field component $|E_x|$ for different incidence angles and radii compared to 3D FEM solutions (right panel). For oblique incidence angle (Fig. 12), it was easier to move the interfaces than to rotate the incident rib port in COMSOL. This shows the flexibility of the presented method in this work, as it is also possible to simulate the rotated structure with the SASM. For clarity, the contour lines are only shown at the levels of 5% and 10% in the following, as the radiated fields produce a rather confusing field pattern. It can be stated from the field plots that the results match well despite the only partially converged mesh of the FEM simulation and the piecewise linear approximations of the interface in the SASM. Further examples for increased incidence angles are presented in Fig. 13, where the side view in the lower panel is shown along the propagation direction of the beam. As can be clearly seen from the side views of the plots, losses can not be neglected for these high refractive index contrast waveguide lenses.
Fig. 11. Field plots of the absolute electric field component $|E_x|$ for $R_1=30\,\mu \textrm {m}, R_2=-30\mu m$ (a) and $R_1=40\mu m, R_2=-40\mu m$ (b) and normal incidence angle $\varphi =0^{\circ }$; top view (upper row) and side view (bottom row). The corresponding full vectorial 3D COMSOL solution is shown on the right side (C). The contour lines mark the levels of the absolute field at 5% and 10%.
Fig. 12. Field plots of the absolute electric field component $|E_x|$ for $R_1=30\,\mu \textrm {m}, R_2=-30\mu m$ and incidence angles $\varphi =15^{\circ }$ (a), $\varphi =30^{\circ }$ (b); top view (upper row) and side view (bottom row). The corresponding full vectorial 3D COMSOL solution is shown on the right side (C). The contour lines mark the levels of the absolute field at 5% and 10%.
Fig. 13. Field plots of the absolute electric field component $|E_x|$ for different radii $R_1,R_2$ and oblique incidence angles $\varphi _0$; top view (upper row) and side view (bottom row). The contour lines mark the levels of the absolute field at 5% and 10%.
5. Concluding remarks
We presented a method to calculate the wave propagation in 3D slab waveguide lenses using the stepwise angular spectrum method (SASM) in combination with full vector 2D solutions. The curved interface is divided into a finite number of plane elements tangential to the surface, from which the radiated field is calculated using the SASM. By superimposing the partial solutions, the total field is determined. The method can be applied to simple 2D settings, consisting of two different materials separated by one or two curved interfaces, as well as for 3D slab waveguide lenses, consisting of slab waveguides of different thicknesses, which are also separated by one or two curved interfaces. In all examples shown in this work we found that our method gives results of comparably accuracy when compared to full 3D FEM solutions with substantially lower computational effort.
Furthermore, our SASM approach might also be interesting for the simulation of wave propagation in other areas of photonics, e.g. surface plasmon polaritons or Bloch surface waves [39–41]. As far as we can assess, the application should be straightforward.
Due to its efficiency, the method could be used to optimize the 3D curved slab waveguide lens regarding any arbitrary figure of merit. Here, we may think of an optimization concerning a small focal length, narrow focal point, no aberrations, etc. The advantage of our method is that sectional areas can be easily evaluated without having to calculate the entire 3D solution, which is useful for an efficient and fast optimization.
Furthermore, the method is feasible enough to take oblique incidence on the curved boundary surface into account, although this is not discussed in detail. With the help of oblique incidence at large angles it is possible to suppress losses completely and to prevent coupling to the other modes, e.g. TM modes. In some of our earlier papers [27–34] the theory behind the oblique incidence has been explained in detail. This effect can potentially also be used as a tool for optimization.
Funding
Deutsche Forschungsgemeinschaft (HA7314/1, TRR142).
Acknowledgments
The authors of this paper acknowledge computing time at PC$_2$ Paderborn Center for Parallel Computing.
Disclosures
The authors declare no conflicts of interest.
References
1. S. Kou, C. Sheppard, and J. Lin, “Calculation of the volumetric diffracted field with a three-dimensional convolution: the three-dimensional angular spectrum method,” Opt. Lett. 38(24), 5296–5298 (2013). [CrossRef]
2. T. Shimobaba, N. Masuda, and T. Ito, “Arbitrary shape surface fresnel diffraction,” Opt. Express 20(8), 9335–9340 (2012). [CrossRef]
3. S. Odate, C. Koike, H. Toba, T. Koike, A. Sugaya, K. Sugisaki, K. Otaki, and K. Uchikawa, “Angular spectrum calculations for arbitrary focal length with a scaled convolution,” Opt. Express 19(15), 14268–14276 (2011). [CrossRef]
4. Y. Wang, S. Yan, A. T. Friberg, D. Kuebel, and T. D. Visser, “Electromagnetic diffraction theory of refractive axicon lenses,” J. Opt. Soc. Am. A 34(7), 1201–1211 (2017). [CrossRef]
5. P. Wu and T. Stepinski, “Extension of the angular spectrum approach to curved radiators,” J. Acoust. Soc. Am. 105(5), 2618–2627 (1999). [CrossRef]
6. K. Yamamoto, Y. Ichihashi, T. Senoh, R. Oi, and T. Kurita, “Calculating the fresnel diffraction of light from a shifted and tilted plane,” Opt. Express 20(12), 12949–12958 (2012). [CrossRef]
7. K. Matsushima, “Shifted angular spectrum method for off-axis numerical propagation,” Opt. Express 18(17), 18453–18463 (2010). [CrossRef]
8. S. De Nicola, A. Finizio, G. Pierattini, P. Ferraro, and D. Alfieri, “Angular spectrum method with correction of anamorphism for numerical reconstruction of digital holograms on tilted planes,” Opt. Express 13(24), 9935–9940 (2005). [CrossRef]
9. J. Stock, N. G. Worku, and H. Gross, “Coherent field propagation between tilted planes,” J. Opt. Soc. Am. A 34(10), 1849–1855 (2017). [CrossRef]
10. K. Matsushima, H. Schimmel, and F. Wyrowski, “Fast calculation method for optical diffraction on tilted planes by use of the angular spectrum of plane waves,” J. Opt. Soc. Am. A 20(9), 1755–1762 (2003). [CrossRef]
11. L. Onural, “Exact solution for scalar diffraction between tilted and translated planes using impulse functions over a surface,” J. Opt. Soc. Am. A 28(3), 290–295 (2011). [CrossRef]
12. S. H. Wiersma, P. Török, T. D. Visser, and P. Varga, “Comparison of different theories for focusing through a plane interface,” J. Opt. Soc. Am. A 14(7), 1482–1490 (1997). [CrossRef]
13. S. Guha and G. D. Gillen, “Vector diffraction theory of refraction of light by a spherical surface,” J. Opt. Soc. Am. B 24(1), 1–8 (2007). [CrossRef]
14. E. Şahin and L. Onural, “Scalar diffraction field calculation from curved surfaces via gaussian beam decomposition,” J. Opt. Soc. Am. A 29(7), 1459–1469 (2012). [CrossRef]
15. G. B. Esmer, L. Onural, and H. M. Ozaktas, “Exact diffraction calculation from fields specified over arbitrary curved surfaces,” Opt. Commun. 284(24), 5537–5548 (2011). [CrossRef]
16. C.-Y. Hwang, S. Oh, I.-K. Jeong, and H. Kim, “Stepwise angular spectrum method for curved surface diffraction,” Opt. Express 22(10), 12659–12667 (2014). [CrossRef]
17. J. Goodman, Introduction to Fourier Optics, 2nd ed (MaGraw-Hill, 1996).
18. M. Born and E. Wolf, Principles of Optics, 7th. ed. (Cambridge University, Cambridge, UK, 1999).
19. G. Kirchhoff, “Zur Theorie der Lichtstrahlen,” Ann. Phys. 254(4), 663–695 (1883). [CrossRef]
20. P. K. Tien, S. Riva-Sanseverino, R. J. Martin, and G. Smolinsky, “Two-layered construction of integrated optical circuits and formation of thin-film prisms, lenses, and reflectors,” Appl. Phys. Lett. 24(11), 547–549 (1974). [CrossRef]
21. G. Hatakoshi, H. Inoue, K. Naito, S. Umegaki, and S. Tanaka, “Optical waveguide lenses,” Opt. Acta 26(8), 961–968 (1979). [CrossRef]
22. W. Chang and P. Ashley, “Fresnel lenses in optical waveguides,” IEEE J. Quantum Electron. 16(7), 744–754 (1980). [CrossRef]
23. D. W. Hewak and J. W. Y. Lit, “Solution deposited optical waveguide lens,” Appl. Opt. 28(19), 4190–4198 (1989). [CrossRef]
24. K. E. Spaulding and G. M. Morris, “Achromatic waveguide lenses,” Appl. Opt. 30(18), 2558–2569 (1991). [CrossRef]
25. K. E. Spaulding, G. M. Morris, and J. C. Brazas, “Achromatic mode-index/diffractive waveguide lens,” J. Lightwave Technol. 10(12), 1855–1859 (1992). [CrossRef]
26. K. E. Spaulding, COMSOL Multiphysics GmbH, Göttingen, Germany. https://www.comsol.de.
27. M. Hammer, A. Hildebrandt, and J. Förstner, “Full resonant transmission of semi-guided planar waves through slab waveguide steps at oblique incidence,” J. Lightwave Technol. 34(3), 997–1005 (2016). [CrossRef]
28. M. Hammer, “Oblique incidence of semi-guided waves on rectangular slab waveguide discontinuities: A vectorial QUEP solver,” Opt. Commun. 338, 447–456 (2015). [CrossRef]
29. M. Hammer, A. Hildebrandt, and J. Förstner, “How planar optical waves can be made to climb dielectric steps,” Opt. Lett. 40(16), 3711–3714 (2015). [CrossRef]
30. A. Hildebrandt, A. Alhaddad, M. Hammer, and J. Förstner, “Oblique incidence of semi-guided waves on step-like folds in planar dielectric slabs: Lossless vertical interconnects in 3-D integrated photonic circuits,” Proc. SPIE 9750, 97501B (2016). [CrossRef]
31. M. Hammer, L. Ebers, A. Hildebrandt, A. Alhaddad, and J. Förstner, “Oblique semi-guided waves: 2-d integrated photonics with negative effective permittivity,” 2018 IEEE 17th International Conference on Mathematical Methods in Electromagnetic Theory, (MMET) pp. 5–9 (2018).
32. L. Ebers, M. Hammer, and J. Förstner, “Oblique incidence of semi-guided planar waves on slab waveguide steps: effects of rounded edges,” Opt. Express 26(14), 18621–18632 (2018). [CrossRef]
33. M. Hammer, L. Ebers, and J. Förstner, “Oblique evanescent excitation of a dielectric strip: A model resonator with an open optical cavity of unlimited Q,” Opt. Express 27(7), 9313–9320 (2019). [CrossRef]
34. M. Hammer, L. Ebers, and J. Förstner, “Oblique quasi-lossless excitation of a thin silicon slab waveguide: a guided-wave variant of an anti-reflection coating,” J. Opt. Soc. Am. B 36(9), 2395–2401 (2019). [CrossRef]
35. L. Ebers, M. Hammer, M. B. Berkemeier, A. Menzel, and J. Förstner, “Coupled microstrip-cavities under oblique incidence of semi-guided waves: a lossless integrated optical add-drop filter,” OSA Continuum 2(11), 3288–3298 (2019). [CrossRef]
36. E. Silberstein, P. Lalanne, J.-P. Hugonin, and Q. Cao, “Use of grating theories in integrated optics,” J. Opt. Soc. Am. A 18(11), 2865–2875 (2001). [CrossRef]
37. S. Thongrattanasiri, J. Elser, and V. A. Podolskiy, “Quasi-planar optics: computing light propagation and scattering in planar waveguide arrays,” J. Opt. Soc. Am. B 26(12), B102–B110 (2009). [CrossRef]
38. M. Lohmeyer and R. Stoffer, “Integrated optical cross strip polarizer concept,” Opt. Quantum Electron. 33(4/5), 413–431 (2001). [CrossRef]
39. L. Yu, E. Barakat, T. Sfez, L. Hvozdara, J. D. Francesco, and H. Herzig, “Manipulating bloch surface waves in 2D: a platform concept-based flat lens,” Light: Sci. Appl. 3(1), e124 (2014). [CrossRef]
40. E. A. Bezus, L. L. Doskolovich, N. L. Kazanskiy, V. A. Soifer, and S. I. Kharitonov, “Design of diffractive lenses for focusing surface plasmons,” J. Opt. 12(1), 015001 (2010). [CrossRef]
41. H. Kim, J. Hahn, and B. Lee, “Focusing properties of surface plasmon polariton floating dielectric lenses,” Opt. Express 16(5), 3049–3057 (2008). [CrossRef]
