Light diffraction in slab waveguide lenses simulated with the stepwise angular spectrum method

: 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 ﬁrst 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 reﬂection and scattering losses, the stepwise angular spectrum method is combined with full vectorial ﬁnite 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.


Introduction
Calculating the diffraction of waves in optical systems has been intensively investigated by means of different theoretical approaches and numerical methods [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][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][18][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][7][8][9][10][11] or light propagation from/to curved surfaces [12][13][14][15][16].
In this work we want to investigate the diffraction of a slab waveguide lens, as illustrated in Figure 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 Figure 1, using a variation of the angular spectrum method.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][21][22][23][24][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 Figure 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.

Single curved interface in 2D
Initially, a 2D interface is considered as illustrated in Figure 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 ϕ 0 with respect to the z-axis.We assume time harmonic fields ∼ exp(iωt), with frequency ω = kc = 2πc/λ for vacuum wavelength λ, 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 ∆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 Figure 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.
By using oblique approximation elements, each solution is calculated in a new local coordinate system (ξ n , η n ) rotated by the inclination angle α of the approximated cell.For each element this angle is given by tanα = g (x n ), where g (x) = d dx g(x) is the spatial derivative of g(x).The local, rotated and shifted coordinate system is then specified by where (x, z) is the global coordinate system.This leads to a rotation of the respective wavenumbers Here, we modified the approach from [16], where straight elements, perpendicular to the longitudinal direction (z-axis), were used to approximate the curved surface.We compare the two approximation techniques in Section 2.2.

Theory
The incoming wave bundle reaches the interface at an oblique primary angle of incidence ϕ 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 ; ϕ 0 ) as For the 2D setting Ψ indicates the principal electric or magnetic field component E y or H y respectively, for TE or TM modes.This integral represents a wave packet calculated over a range of incidence angles ϕ (wavenumbers k x ).The last two exponential terms depict the propagation in xand z-direction and (x 0 , z 0 ) positions the focus of the beam.The first weighting factor defines the shape of the incoming field in lateral direction.Here, we assume a plane wave of Gaussian shape with primary wavenumber k x 0 = k 0 n 1 sinϕ 0 (or incidence angle ϕ 0 ) and Gaussian weight of half width w k = 4cosϕ 0 /W c , where W c is the full lateral beam width perpendicular to the propagation direction [27].
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 Note that in the local coordinate system the incoming wave reaches the surface at an incidence angle ϕ 0 + α.Therefore, the primary incident wavenumber is given by k x 0 = k 0 n 1 sin(ϕ 0 + α) and the Gaussian weight by w k = 4cos(ϕ 0 + α)/W c .The focus (x 0 , z 0 ) of the incident wave is adapted to the new coordinate system and since the origin of the coordinate system lies in (x n , g(x n )), the field on the interface is given at the point z = 0. Furthermore, the Jacobian matrix of the coordinate transformation is unity and thus omitted.
For the known incoming field on each of the n elements in local coordinates Ψ in,n (x ), the next step is to calculate its angular spectrum by the inverse Fourier transformation with ∆x = ∆x/cosα and the rectangular window function is defined as As already mentioned in [16], the rectangular window function could also be replaced by functions that achieve similar behavior, e.g. a function of Gaussian shape.
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 where the negative sign refers to forward (transmitted) propagating waves and the positive sign to backward (reflected) propagating waves with amplitudes a i as introduced below.Finally, the overall solution is then determined by summing up the intermediate results Here, we always indicated with i the reflected i = r or transmitted i = t field, but note that this approach is also valid to additionally approximate the incoming i = in field from Eq. ( 4) via detours.Depending on which field is considered, the wavenumbers are assumed to be given by k x = k 0 n i sinϕ i and k z = k 0 n i cosϕ i .This implies that for the incoming and reflected field n in = n r = n 1 , ϕ in = ϕ r = ϕ and for the transmitted field Analogously to a simple plane wave interface the wave amplitude functions a i (k x ) for this 2D setting are defined by the Fresnel coefficients.Since the shape of the incoming wave is given by w g , or A n (k x ), respectively, the incoming amplitude a in is set to one.The coefficients for the reflected waves a r (k x ) = r(ϕ(k x )) are then given by depending on the polarization of the incoming wave.For the transmitted field a t (k x ) = t(ϕ(k x ))), the coefficients are defined by Our approach also allows k z to become imaginary, so that at high angles of incidence, leading to total reflection, evanescent waves are generated in the corresponding area.Multiple reflections at the single, i.e.same, interface are not included in our method, hence the shape of the curved interfaces is restricted to that.

Examples for a single interface
We now discuss some examples using the 2D SASM as shown in Figure 3. Here, we assume a simple circular curvature function g , 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 Ψ = E y , which reaches the curved surface at an incidence angle of ϕ 0 .The Gaussian beam width is set to W c = 10µm and the structure is operating at telecommunication wavelength λ = 1.55µm.The integrals (here and in the following) required to calculate the solutions were evaluated using numerical quadrature.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 µm and incidence angle ϕ 0 .The mesh is given by a maximum element size of λ/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 ϕ 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 µm or R = 40 µm is calculated as f = 30µm or f = 60µm, respectively.The field plots in Figure 3 confirm this, even if the focal point is slightly distorted due to spherical aberration.

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 Figure 4. Now, multiple reflections take place between the interfaces, resulting in a superposition of forward and backward propagating waves.
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.

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 Figure 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 ϕ 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 (x) and g m (x) with center x n or x m and width ∆x n or ∆x m , respectively.Note, that it is not necessary to use the same discretization for both surfaces.α 1 refers to the angle of inclination for the first interface and α 2 refers to the second.Furthermore, we introduce two different local coordinate systems, (x , z ) for the first and (x, ẑ) 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 Respective wavenumbers are related as and The angular spectrum on the first interface composes now the angular spectrum of the incoming wave and a superposition of the angular spectra from all reflected fields on the second surface, thus Here, A n,0 (k x ) = A n (k x ) is defined by Eq. ( 7), describing the angular spectrum from the incoming wave, and the second term is given by a Fourier transformation (equivalent to Eq. ( 7)) where the incoming field on the first surface is calculated via Note that calculating the angular spectrum in Eq. ( 23) depends on A m,ν ( kx ), which is the angular spectrum from the second interface.
The superimposed radiated forward and backward fields from the first interface are where the negative sign refers to the transmitted waves and the positive sign to the reflected waves.The coordinate transformation from (x , z ) to (ξ n , η n ) is again given by Eq. ( 2) with the inclination angle α 1 .Furthermore, the amplitudes are 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 are where the fields on the interface are calculated via with the amplitude The radiated forward and backward fields from the second interface are then Again, the coordinate transformation Eq. ( 2) from (x , z ) to (ξ m , η m ) enters, but with the inclination angle α 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µm and the operating wavelength by λ = 1.55 µm.The incoming wave is the TE-mode with principal field component Ψ = E y , which reaches the curved surface at an incidence angle of ϕ 0 .
To keep the simulation time as low as possible, we calculated only a finite number ν < ∞ 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 ϕ 0 = 0 • .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 ϕ 0 > 0 for different curvature radii R 1 , R 2 (upper row).Again, we observe good agreement with the corresponding COMSOL solutions (bottom row).
The focus f of the outgoing beam can be calculated by the Lensmaker's equation for a thick lens [18] with For the structures considered in Figure 5, this predicts the focal points at positions z ∈ {14, 16.43, 14.62, 22.7}µm in the order from plot (a) to (d).Again, because of spherical aberrations the focus point is slightly distorted.In the same way, the values of f match also for Figure 6.  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.

Theory
The curved interface from Section 2.1 is again considered (Figure 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 (Figure 1 and Figure 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 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 ϕ 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 ϕ.More precisely, the previously presented method described for the 2D structure in Section 2.1 is applied to the configuration of Figure 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][28][29][30][31][32][33][34][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 eff 1 , n 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 Ψ 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 eff 1 and n 2 by n 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 ϕ: Here, {Φ 0 (y; k x )exp(−ik z η n ) + } =: Φ C (y, z; k x ) depicts the former quasi-2D COMSOL solution for incidence wavenumber k x , Φ 0 indicates the vectorial mode profile for oblique incidence of the incoming wave and is a remainder.Different from the 2D case, Ψ n already includes the incoming, reflected and transmitted field.The total field is then again given by summing up the n individual solutions 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 ϕ, 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.

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 (Figure 7 b).
Generally, in both slab waveguides (incoming and outgoing) the fields can be written as a superposition of all guided slab waveguide modes Ψ m with specific amplitudes a m , wavenumbers k z,m and scattered field Ψ sc where m counts the guided slab modes in either the incoming or outgoing slabs.
The amplitudes a m (k x ) can be extracted from the COMSOL solution Φ C by using a mode overlap product [38] (Ψ; Φ) defined for mode orthogonality with the transverse field components of two different modes Ψ = (E 1 , H 1 ) and Φ = (E 2 , H 2 ).Applying this product for each slab waveguide mode Ψ n (Eq.( 36)) with n ∈ [1, m] at position z = z t this results in Assuming mutual orthogonality of the guided slab modes, orthogonality of any slab mode and the scattering field (Ψ n ; Ψ sc ) y = 0 with n ∈ [1, m], and power normalization, the amplitudes can be extracted as By extending the 2D solution with its corresponding slab waveguide modes, the SASM simulation window is completely independent of the COMSOL simulation window.Moreover, the method has a substantial storage advantage in comparison to calculating the entire 2D solution for large propagation distances, and also the computational time required to compute these fields is significantly lower.In the following, we use the full 2D COMSOL solution up to the point where out-of-plane radiation is visible.Behind this area, only the semi-guided slab modes are used to assemble the field for the plots as shown.Up to this point, almost no approximations were applied when expanding the 2D field solution.

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].
Eq. ( 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 Here, Ψ 0 (y; k x ) depicts the former 1-D vectorial mode profile of the fundamental TE slab solution along the y-axis for different incidence angles.The mode profile of the rib waveguide of the fundamental TE-like mode is denoted as Φ 0 (x, y).To determine the unknown weight w 0 (k x ) the mode profile Φ 0 of the rib can be written as a superposition of the complete set of normal modes Ψ j , supported by the slab waveguide at the junction plane where rib and slab waveguide meet (here, z = z 0 ), where j identifies mode orders.By using the product (Ψ 0 , •) y (Eq.( 37)) and assuming that the modes are orthogonal and power normalized, the unknown weight for TE 0 incidence is given by For simplicity, we restrict the approach to normal incidence ϕ 0 = 0 • , at first, to obtain the new rib weight w 0 (k x ) = w 0 (k x ; 0).But the calculation can also be done analogously for oblique incidence by rotating the coordinates and fields by the appropriate angle ϕ 0 resulting in w 0 (k x ; ϕ 0 ).

Examples -Curved slab waveguide discontinuity
To present some examples for a 3D curved interface, two slabs with thicknesses d 1 = 0.05µm and d 2 = 0.22µ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 − √ 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µm, height 0.05µm (same height as the adjoining incoming slab waveguide) and etch depth 0.02µm.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 ∈ [−5, 20]µ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 ϕ 0 > 0 are shown in Figure 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.
Already for a single interface in 3D, comparably high losses occur.To get an impression of the magnitude of loss, Figure 10 shows the relative power transmittance T, the reflection R and the total outgoing power to guided modes P out depending on the incidence angle ϕ 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][28][29][30][31][32][33][34].

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µm for the incoming and outgoing slab and d 2 = 0.22µ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µm, height 0.05µm and etch depth 0.02µm.The constituting 2D COMSOL solutions are again calculated in the range z ∈ [−5, 20]µm and extended with the method from Section 4.1.1.For different radii R 1 , R 2 and angle of incidence ϕ 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.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 Figure 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.

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][40][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][28][29][30][31][32][33][34] the theory behind the oblique incidence has been explained in detail.This effect can potentially also be used as a tool for optimization.

Figure 1 :
Figure 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.

Figure 2 :
Figure 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 ϕ 0 .The dashed red line represents the discretization of the curved interface.

Figure 3 :
Figure 3: Field plots of the absolute electric field component |E y | for different curvature radii R and incidence angles ϕ 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%.

Figure 4 :
Figure 4: Configuration with two curved interfacesgiven by g(1) (x) and g(2) (x).The incoming wave reaches the first surface under an incidence angle ϕ 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.

Figure 5 :
Figure 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 ϕ 0 = 0 • .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%.

Figure 6 :
Figure 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 ϕ 0 > 0 • .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%.

Figure 7 :
Figure 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.

Figure 8
Figure 8 shows the calculated field solutions for the absolute electric field value |E x | for curvature radii R = 20µm (a) and R = 30µ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 ∈ [−5, 20]µ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 α) and for incidence angles from ϕ = −85 • to ϕ = 85 • with steps of 0.5 • .

Figure 8 :
Figure 8: Field plots of the absolute electric field component |E x | for curvature radius R = 20 µm (a) and R = 30µm (b) and normal incidence angle ϕ = 0 • ; 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%.

Figure 9 :
Figure 9: Field plots of the absolute electric field component |E x | for different curvature radii R and incidence angles ϕ 0 > 0 • ; top view (upper row) and side view (bottom row).The contour lines mark the levels of the absolute field at 2%, 5% and 10%.

Figure 10 :
Figure 10: Reflection R and transmittance T for fundamental guided TE-/TM-polarization versus the incidence angle ϕ of the 2D constitutive solutions for either incoming TE (a) or TM (b) mode.

Figure 11 :
Figure 11: Field plots of the absolute electric field component |E x | for R 1 = 30 µm, R 2 = −30µm (a) and R 1 = 40µm, R 2 = −40µm (b) and normal incidence angle ϕ = 0 • ; 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%.

Figure 12 :
Figure 12: Field plots of the absolute electric field component |E x | for R 1 = 30 µm, R 2 = −30µm and incidence angles ϕ = 15 • (a), ϕ = 30 • (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%.

Figure 11 and
Figure11and Figure12show 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 (Figure12), 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 Figure13, 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.

Figure 13 :
Figure 13: Field plots of the absolute electric field component |E x | for different radii R 1 , R 2 and oblique incidence angles ϕ 0 ; top view (upper row) and side view (bottom row).The contour lines mark the levels of the absolute field at 5% and 10%.