1. Introduction

Photonic integration has evolved into a key technology for a wide variety of applications that range from high-speed communications [1], ultra-fast signal processing [2,3] and artificial intelligence [4] to optical metrology and sensing [5–8] and to biophotonics and life sciences [9–11]. On the technological level, photonic integrated circuits (PIC) predominantly rely on planar structures that can be fabricated with well-established microfabrication techniques based on layer deposition and 2D patterning via high-resolution electron-beam or deep-UV lithography. More recently, these techniques have been complemented by multi-photon lithography that allows for in-situ fabrication of functional 3D freeform structures that can greatly enhance the functionality and versatility of planar PIC. Examples are 3D-printed chip-chip connections, so-called photonic wire bonds [12,13] that open an attractive path towards high-performance hybrid multi-chip modules [14], 3D-printed waveguide overpasses [15], reconfigurable photonic circuits [16,17], 3D-printed waveguide interconnects for photonic neural networks [18], 3D-printed waveguide splitters [19], or 3D-printed polarization splitters and rotators [20]. However, while simulation tools for planar lightwave structures are available as part of commercial software packages [21–24], efficient modeling and design of 3D freeform waveguides still represents a challenge. This applies in particular to numerical solvers that rely on rectilinear grids within cuboid-shaped computational domains, which is, e.g., the case for most time-domain techniques. Applying such solvers to 3D freeform waveguides with strongly curved non-plane trajectories requires comparatively large computational volumes, within which the structure of interest only occupies a small fraction, thus leading to poor computational efficiency. In addition, accurate representation of curved surfaces on rectilinear grids requires local refinement of the mesh cells, which leads to reduced time steps in time-domain simulations and thus increases overall simulation time.

Here, we present a transformation-optics (TO) method for reducing the computational effort associated with freeform-waveguide simulations. Our approach relies on transforming curved freeform waveguides in the original 3D space into straight waveguides in a virtual 3D space, which can then be efficiently treated by rigorous time-domain Maxwell’s equations solvers defined on a rectilinear grid. Specifically, the transformed waveguide in the virtual space can be confined in a rectangular simulation box, whose volume is comparable to the actual freeform waveguide volume. Furthermore, in case of freeform waveguides with rectangular cores, the grid lines in the virtual space are perfectly aligned with the core surfaces thereby eliminating the need for any local mesh cell refinement. We demonstrate the viability of the concept using the commercially available time-domain solver of CST Microwave Studio (CST MWS), which is based on the finite integration technique (FIT) [25,26], reaching an acceleration by a factor of 3–6. In addition, we fabricate the simulated freeform waveguides on a silicon photonic (SiP) chip and measure the transmission losses at a vacuum wavelength of $\lambda = {1550}\,\mathrm{nm}$. We find excellent agreement between TO-based simulations in virtual space and the associated reference simulations in real space, and we confirm that the simulated transmission losses match their experimentally measured counterparts reasonably well. Although primarily aimed for use with time-domain solvers on a rectilinear grid, our method represents a general technique that allows to transform freeform waveguide-based devices into straight structures and is independent of the underlying solver.

2. TO based concept of freeform waveguide modeling

The TO concept relies on the fact that Maxwell’s equations are form-invariant with respect to coordinate transformations. In particular, if we map an original domain from an $\left (x,y,z\right )$-coordinate system, to a virtual domain in a $\left (u,v,s\right )$-coordinate system, we only need to adapt the material properties in the virtual domain, while the form of Maxwell’s equations remains unchanged [27–30]. In case of a coordinate transformation described by a differentiable function $\left (u,v,s\right )^{\textrm{T}} = \mathbf {f}\left (x,y,z\right )$, where $\left (x,y,z\right ),\left (u,v,s\right )\in \mathbb {R}^{3},$ the relationship between the material properties in the original $\left (x,y,z\right )$-space and in the virtual $\left (u,v,s\right )$-space reads [27–30]:

The TO concept has previously been used to analyze and design a variety of devices such as beam deflectors and expanders [31], polarization splitters and rotators [32], flat lenses [33], waveguide bends [34], or electromagnetic cloaks for hiding of objects [35–37] or for reshaping the perception of cloaked objects [38,39]. In general it is possible to design TO-based devices with arbitrary shapes [40] including also plasmonic structures [41], and the underlying principles of TO can be further expanded into thermodynamics and mechanics [42] or general relativity [43]. Here, we exploit TO for efficient numerical modeling of 3D freeform waveguides using well-established time-domain solvers, see Fig. 1 for an illustration of the underlying concept. Generally, time-domain solvers relying, e.g., on finite-difference-time-domain (FDTD) techniques, are perfectly suited for large-scale simulations, offering a numerical complexity that scales linearly with problem size while being amenable to efficient parallelization on large-scale computer clusters. Moreover, time-domain techniques are robust, lend themselves to broadband or transient simulations, and even offer a natural path to handling of nonlinear behavior [44]. On the other hand, time-domain techniques usually rely on rather inflexible rectilinear grids and cuboid-shaped computational domains, which severely limits the performance when applied to 3D freeform waveguides. Specifically, the rectilinear grid does not allow for efficient representation of curved surfaces, while the cuboid-shaped computational domain leads to poor computational efficiency with comparatively large computational volumes within which the structure of interest only occupies a small fraction, see Fig. 1(a).

 

Fig. 1. Transformation of a freeform optical waveguide from the original $\left (x,y,z\right )$-space to the virtual $\left (u,v,s\right )$-space. (a) Sample freeform waveguide in the original $\left (x,y,z\right )$-space. The computational domain necessary to simulate the structure with a time-domain solver on a rectilinear grid is determined by the rectangular red box. However, the region of interest, i.e., the relevant part of the computational domain, only comprises the close vicinity of the freeform waveguide — this part and its edges are depicted in blue. The rest of the computational domain is unimportant and unnecessarily consumes computing resources. Moreover, a fine discretization would be needed to correctly represent the curved surfaces of the freeform waveguide by a rectilinear computational grid. The coordinate transformation described in Eq. (3) is defined by the unit vectors $\mathbf {U}$, $\mathbf {V}$, and $\mathbf {T}$, where $\mathbf {T}$ is the local tangent unit vector $\text {d}\mathbf {r}/\text {d}s$ of the trajectory, while $\mathbf {U}$ and $\mathbf {V}$ span the transverse plane in the respective trajectory point. The vectors are chosen such that $\left ( \mathbf {U},\mathbf {V},\mathbf {T} \right )$ is a right-handed trihedron. (b) Transformed waveguide in the virtual $(u,v,s)$-space. The freeform waveguide trajectory in the original space has been mapped to a straight line, and the relevant part of the computational domain that was "twisted" in the original $\left (x,y,z\right )$-space is now represented by a rectangular blue box. The computational domain can now be reduced to the region of interest, see red dashed rectangular box, and this means a volume reduction by a factor of about 25. The coordinate transformation significantly reduces the computational domain at the cost of pre-calculating the spatial distributions of tensors of dielectric permittivity and magnetic permeability in $\left (u,v,s\right )$-space. Moreover, the rectilinear computational grid may be better adapted to the shape of the waveguide, in particular for rectangular cross sections.

Download Full Size | PDF

To overcome these problems, we map the freeform waveguide with a curved trajectory in the original $\left (x,y,z\right )$-space to an equivalent waveguide with a straight trajectory and modified permittivity and permeability tensors in a virtual $\left (u,v,s\right )$-space, see Fig. 1(b). In the virtual space, we may then use a rectangular computational domain that only encompasses the straight waveguide and its direct vicinity, along with a rectilinear grid that is oriented along the direction of the waveguide. The virtual waveguide can thus be efficiently modeled by a conventional time-domain solver, and the results are then transformed back to the $\left (x,y,z\right )$-coordinate system to obtain the field distributions in original space.

To implement this technique, we need to define a function $\left (u,v,s\right )^{\textrm{T}} = \mathbf {f}\left (x,y,z\right )$ that transforms the curved waveguide in original space into a straight path in virtual space. It is actually easier and more intuitive to analytically express the inverse function $\left (x,y,z\right )^{\textrm{T}} = \mathbf {f}^{-1}\left (u,v,s\right )$ that maps a point $\left (u,v,s\right )$ in transformed space back to original space. To arrive at a mathematical formulation of $\mathbf {f}^{-1}$, we assign the coordinate $s$ to the arc length of the waveguide trajectory $\mathbf {r}\left (s\right ) = \left (x_0\left (s\right ),y_0\left (s\right ),z_0\left (s\right )\right )$ in original space, while $u$ and $v$ are associated with the transverse position relative to the waveguide trajectory, where the direction is defined by a pair of unit vectors $\mathbf {U}$ and $\mathbf {V}$. This leads to the relation

Equation (3) still leaves the freedom of choosing the orientation of the unit vectors $\mathbf {U}$ and $\mathbf {V}$ within the transverse plane in the respective trajectory point. One obvious choice for the $\left ( \mathbf {U},\mathbf {V},\mathbf {T} \right )$ frame could be the natural Frenet-Serret frame, where $\mathbf {U}$ could be chosen as the binormal vector, and $\mathbf {V}$ could be chosen as the normal vector of the trajectory $\mathbf {r}\left (t\right )$. However, the Frenet-Serret frame is not the best choice because the binormal and the normal vectors are neither defined in points where the trajectory is straight, nor in inflection points. Furthermore, the frames immediately before and immediately after an inflection point are rotated against each other by $180^{\circ }$ about the tangent vector. We therefore use the rotation-minimizing frame (RMF) [45,46], which minimizes the frame spinning along the trajectory, and which is commonly used in computer graphics and 3D modeling [47]. To calculate the RMF, we use a simple and fast approximation method called double reflection method [48]. This method requires the trajectory sample points, the tangent vectors $\mathbf {T}$ in all sample points, and a coordinate frame in the first sample point $\left (\mathbf {U}_1,\mathbf {V}_1,\mathbf {T}_1\right )$ as an input. The coordinate frames in the remaining sample points along the trajectory are calculated recursively [48]. The coordinate frame $\left (\mathbf {U}_1,\mathbf {V}_1,\mathbf {T}_1\right )$ in the first sample point is chosen such that $\mathbf {T}_1$ is the tangent, and the two remaining mutually perpendicular vectors $\mathbf {U}_1$ and $\mathbf {V}_1$ can be chosen arbitrarily in the plane that is perpendicular to $\mathbf {T}_1$.

Note that the coordinate transformation described in Eq. (3) is generally not bijective. A trivial example is the case of self-intersecting trajectories, which must be excluded. Another obvious problem might arise if two waveguide segments pass by each other or overpass each other within a close distance, e.g., in spirals, helices, and loops. In this case, care must be taken to avoid mapping the same sub-domain of the $\left (x,y,z\right )$-space being in the vicinity of both waveguide segments twice into two distinct sub-domains of the $\left (u,v,s\right )$-space. Finally, bijectivity might be violated if the local curvature of the trajectory is too strong, such that the center of curvature falls into the transformed domain. Specifically, for a given trajectory point, the ranges for $u$ and $v$, $u\in [u_{min},u_{max}]$ and $v\in [v_{min},v_{max}]$ should be chosen not to include the center of curvature for this specific trajectory point. In other words: If, in a certain trajectory point, the center of the curvature would be represented by $\left (u_{c},v_{c}\right )$ in transformed space, then the range of $u$ and $v$ for this trajectory point must be limited not to include this point, which is ensured by the sufficient condition $\left (u_{c},v_{c}\right ) \notin \{\left (u,v \right )\vert \: u_{min} \le u \le u_{max} \; \wedge \; v_{min} \le v \le v_{max}\}$. For more details, see Appendix A.

It should be noted that the TO technique represents a very general concept, relying only on the fact that Maxwell’s equations are form-invariant with respect to coordinate transformations. The approach is hence generally independent of material properties and can be applied to a wide range of structures that can be significantly more complex than the ones investigated in our manuscript. Examples comprise, but are not limited to, multi-layer waveguide structures made from doped semiconductor materials with non-vanishing imaginary part of the complex electric permittivity, metal waveguides and plasmonic structures [41], perfect electric conductors (PEC) [49], or refractive optical elements such as lenses. In the latter case, curved lens surfaces in the original space may be transformed to rectangular cuboids in the transformed space. This would not lead to a reduction of the computational volume as for the case of freeform waveguides considered in our work but may still avoid local mesh refinement in the vicinity of curved surfaces and thereby allow to keep larger time steps. The overall reduction of computational complexity will strongly depend on the respective structure.

3. Implementation of TO modeling of freeform waveguides

To prove the viability of the proposed approach, we implement and use the TO concept to simulate freeform waveguides with an invariant rectangular cross section of the core. For the coordinate transform, we use MATLAB, and the resulting waveguide in virtual space is then treated with CST MWS, which relies on a time-domain finite-integration technique (FIT) [25,26]. We restricted the implementation to freeform waveguides that are made from isotropic material and that feature plane trajectories, such that the tensors $\boldsymbol{\mathrm{\epsilon}}^{\prime }$ and $\boldsymbol{\mathrm{\mu}}^{\prime }$ are diagonal in the transformed space, see Appendix B. Anisotropic materials and/or a non-plane trajectory would lead to non-diagonal tensors $\boldsymbol{\mathrm{\epsilon}}^{\prime }$ and $\boldsymbol{\mathrm{\mu}}^{\prime }$, which cannot be treated by the FIT solver of CST MWS (CST Studio Suite version 2019). Note that the treatment of such non-diagonal tensors associated with non-plane trajectories increases the computational effort by a factor of approximately two compared to diagonal tensors associated with plane trajectories, while offering a larger reduction of the computational volume, e.g., by a factor of 25 in the example illustrated in Fig. 1. A more detailed discussion of the impact of non-plane trajectories on computational complexity can be found in Appendix D.

The input data for our MATLAB code consist of the trajectory sample points $\left (x_0,y_0,z_0\right )$ in the original $\left (x,y,z\right )$-space, the orientation of the RMF in the initial point, the cross-sectional shape of the waveguide core along with the material properties of the core and the cladding ($\varepsilon _\mathrm {core}$, $\mu _\mathrm {core}$, $\varepsilon _\mathrm {cladding}$, and $\mu _\mathrm {cladding}$), and the ranges of $u$- and $v$-coordinates that define the space to be transformed. Without loss of generality, we can assume that the plane trajectory lies in the $(y,z)$-plane, i.e., that the $x$-coordinates of all trajectory points are equal to zero. For simplicity, we chose the right-handed RMF $\left (\mathbf {U}_1,\mathbf {V}_1,\mathbf {T}_1\right )$ in the first sample point of the trajectory such that $\mathbf {U}_1$ is parallel to the $x$-axis and $\mathbf {V}_1$ is perpendicular to $\mathbf {U}_1$ and to the tangent vector $\mathbf {T}_1$. We further assume a rectangular core cross section with width $w_{\text {WG}}$ and height $h_{\text {WG}}$, measured along the $\mathbf {U}$- and $\mathbf {V}$-directions of the RMF. In a first step of the transformation, we numerically calculate the $s$-coordinate sample points from the trajectory sample points. The RMF $\left (\mathbf {U}\left (s\right ),\mathbf {V}\left (s\right ),\mathbf {T}\left (s\right )\right )$ in the remaining trajectory points is then calculated numerically by the aforementioned double reflection method [48]. In the special case of a plane trajectory in the $(y,z)$-plane with $\mathbf {U}_1$ chosen to be parallel to the $x$-axis, $\mathbf {U}\left (s\right )$ is parallel to the $x$-axis in all trajectory points. Having calculated the $s$-coordinates with given $u$- and $v$-coordinate ranges, we have the rectangular box defining the computational domain in $(u,v,s)$-space. The ranges of $u$ and $v$ determine how much space around the trajectory will be taken into account for the TO-based simulations. These ranges should be large enough to include a sufficient portion of the evanescent field outside the core, but not too large, to ensure bijectivity of the function $\left (u,v,s\right )^{\textrm{T}}=\mathbf {f}\left (x,y,z\right )$, see Appendix A. In the next step, the transformed material properties in virtual space need to be calculated and fed to the numerical solver. Generally, this can be done by adapting $\boldsymbol{\mathrm{\epsilon}}^{\prime }$ and $\boldsymbol{\mathrm{\mu}}^{\prime }$ in each point of space according to Eq. (1). For commercially available simulation programs with CAD-type user interfaces, however, setting the spatial variations of the materials properties on the level of individual spatial grid cells is not very efficient. We therefore approximate the transformed structure with continuously varying material properties by a multitude of small bricks with constant material properties, which can be fed to the solver as cuboid CAD elements by a standard scripting interface. We partition the range of $u$-coordinates into $N_u$ steps, the range of $v$-coordinates into $N_v$ steps, and the range of $s$-coordinates into $N_s$ steps. Then the total number of bricks is given by

Assuming a simplified waveguide structure with a plane trajectory that entirely lies in the $(y,z)$-plane allows to greatly simplify the Jacobian according to Eq. (2), see Eq. (9) of Appendix B. In this case, regions for which the material properties of the original waveguide in original space do not change along $x$ do not need to be subdivided into bricks along $u$. Since we consider a waveguide with homogeneous core and cladding region, we may thus reduce the transformed structure to $N_u^{\text {core}} = 3$ bricks along $u$ in the core region and $N_u^{\text {cladding}} = 1$ brick along $u$ in the cladding region below and above the core, see Fig. 2. The range of $v$-coordinates is divided into $N_v = N_v^{\text {core}} + N_v^{\text {cladding}}$ steps, which comprises $N_v^{\text {core}}$ steps in the core region, and $N_v^{\text {cladding}}$ steps in the cladding region below and above the core. Therefore, the number of bricks per slice in this case is $N_u^{\text {core}}N_v^{\text {core}} + N_u^{\text {cladding}}N_v^{\text {cladding}}$, and Eq. (4) becomes

 

Fig. 2. Representation of a waveguide section with a trajectory lying in the $(y,z)$-plane by rectangular bricks with constant material properties in virtual $\left (u,v,s\right )$-space. The freeform waveguide is discretized into sections with approximately constant curvature of the waveguide trajectory (slices) along $s$, defined by coordinates $s_i$ and $s_{j}$. Green: Core region of the transformed freeform waveguide. Cyan: Cladding region of the transformed freeform waveguide. For a freeform waveguide with a plane trajectory in the $(y,z)$-plane, the axis $\mathbf {U}$ of the RMF is parallel to the $x$-axis in all trajectory sample points, which allows us to completely omit the discretization along $u$ above and below the freeform waveguide core and to reduce the representation along $u$ to three segments in the region of the waveguide core. The illustrated freeform waveguide slice is represented by an overall of 28 bricks — 5 each above and below the core, 6 on the each right and the left of the core, and 6 within the core.

Download Full Size | PDF

Finally, to calculate the material properties of the various bricks in $\left (u,v,s\right )$-space, we first need to know the material properties of the corresponding bricks in $\left (x,y,z\right )$-space and then transform them using Eq. (1). After calculating the material properties of all bricks, the MATLAB code then writes a script which is loaded and run in CST MWS to generate the transformed structure for simulation by the time-domain solver of CST MWS. After the simulation is done, scattering parameters (S-parameters) of the freeform waveguide can be read directly, while the field distribution must additionally undergo the inverse coordinate transformation given by Eq. (3), in order to be represented in the original $\left (x,y,z\right )$-space. This is done by exporting the field distribution from the CST MWS simulation, and performing the inverse transformation by an additional MATLAB script. Note that our technique leaves the underlying solver unchanged and may also allow to exploit existing FDTD implementations [21–24] to efficiently treat 3D freeform waveguides that comprise nonlinear or dispersive media [50,51].

4. Freeform waveguide simulations and experimental benchmarking

To verify our TO approach, we apply it to a series of freeform waveguide structures that connect two SiP waveguides, see Fig. 3(b). We calculate the transmission and the electric field distribution for these structures at a wavelength of $\lambda = {1550}\,\mathrm{nm}$ using the TO concept with CST MWS as a FIT time-domain solver for the waveguide in the virtual $\left (u,v,s\right )$-space, and we compare the results to conventional simulations of the original structure in $\left (x,y,z\right )$-space. Moreover, we experimentally realized and characterized the structures — Fig. 3 shows the experimental setup, some freeform waveguide illustrations, scanning electron microscope (SEM) pictures of the fabricated structures, as well as a comparison of the simulated and measured transmission for different waveguide trajectories, that are described by the height $h$ of the apex above the substrate. As a further reference, we also calculated the transmission using our previously reported fundamental-mode approximation (FMA) method [52]. The FMA is based on a look-up table of pre-calculated eigenmodes (propagation constants and modal field distributions) of waveguide segments with constant radii of curvature, on their tabulated bending losses and on transition losses due to modal field mismatch between two adjacent segments with different radii of curvature. The FMA subdivides freeform waveguides in segments with constant radii of curvature and calculates the transmission losses based on the look-up table. While exceptionally fast, this method disregards the existence of higher-order modes.

 

Fig. 3. Benchmarking of the TO approach with respect to different simulation techniques and to measurements. (a) Experimental setup for measuring the losses of freeform waveguides with different trajectories: Continuous wave (CW) light at a wavelength of $\lambda = {1550}\,\mathrm{nm}$ is launched to the device under test (DUT) consisting of 3D freeform waveguides that are connected to on-chip access waveguides. The light is coupled to the chip by single-mode fibers (SMF) and grating couplers (GC). The polarization of the incoming light is adjusted by a polarization controller (PC), and the power of the transmitted signal is measured by an optical power meter (OPM). (b) Artist’s impression of the test structures on silicon photonic (SiP) chip. The 3D freeform waveguides feature different apex heights $h$, measured between the trajectory in the center of the waveguide and the top surface of the buried oxide (BOX) layer. (c) Side view and (d) top view of the tapered transition between a freeform waveguide and a SiP waveguide. (e) Scanning electron microscope (SEM) image of a series of freeform waveguides fabricated on a SiP chip. The apex height $h$ of the trajectory is swept between $h_{\text {min}}={6.2}\,\mathrm{\mu} \mathrm{m}$ and $h_{\text {max}}={16.2}\,\mathrm{\mu} \mathrm{m}$. Bottom inset: Close-up of a freeform waveguide having an apex height of $h \approx {13}\,\mathrm{\mu} \mathrm{m}$. The central part (red) is 55 µm long. The coupling sections and the adjacent sections (cyan) are kept the same for all freeform waveguides. (f) Simulated and measured transmission of the freeform waveguides for different apex heights $h$. The deviations of the transmissions predicted by the TO technique and the reference simulation (Ref.) range from $-0.5$ dB to $+0.6$ dB, with the biggest differences occurring at the extreme values $h_{\text {min}}$ and $h_{\text {max}}$ of the apex height $h$ where the waveguide curvature is strongest. Methods: TO — Transformation optics, Ref. — Reference simulation in original space, FMA — Fundamental-mode approximation, Exp. — Experimental. (g) Freeform waveguide trajectories with indicated radii of curvature $R_0$, $R_1$, and $R_2$.

Download Full Size | PDF

For the experimental benchmark, the freeform waveguides were 3D-printed on a SiP chip, Figs. 3(b)–3(e), that was fabricated through a commercial foundry using standard CMOS process and ${248}\,\mathrm{nm}$ deep-UV lithography. The waveguides are normally covered by a $\mathrm {SiO_{2}}$ top cladding layer, which was locally removed to make the SiP waveguides accessible to the 3D printing system. For fabrication of the 3D freeform waveguides, a negative-tone photoresist is deposited on the chip and structures are then 3D-printed by two-photon polymerization [53]. Subsequently, the unexposed resist is removed in a separate development step, and the freeform waveguides are covered by a low-index polymer (not shown in Fig. 3) that serves as a cladding and a protection against environmental influences. The refractive index of the freeform waveguide core at $\lambda = {1550}\,\mathrm{nm}$ amounts to $n_{\mathrm {core}} \approx 1.53$, and the cladding refractive index is $n_{\mathrm {cladding}} \approx 1.36$. Both materials are assumed to be lossless dielectrics (relative magnetic permeability $\mu _r=1$ in the original $(x,y,z)$-space), and the corresponding values of the dielectric permittivity are calculated by squaring the refractive indices. Each freeform waveguide core has a $w_{\text {WG}} \times h_{\text {WG}} = {2}\,\mathrm{\mu} \mathrm{m} \times {1.8}\,\mathrm{\mu} \mathrm{m}$ rectangular cross-section and bridges a $l_\text {g} = {100}\,\mathrm{\mu} \mathrm{m}$ gap between a pair of SiP strip waveguides (width $w_{\mathrm {Si}} = {500}\,\mathrm{nm}$, height $h_{\mathrm {Si}} = {220}\,\mathrm{nm}$) that lie on a buried oxide (BOX) $\mathrm {SiO_2}$ layer with a height of $h_{\mathrm {BOX}} = {3}\,\mathrm{\mu} \mathrm{m}$. Each SiP waveguide is connected to a grating coupler (GC) on one side and to a freeform waveguide on the other side. The connections between SiP waveguides and freeform waveguides are made through linear inverse tapers with a length of $l_{\mathrm {taper}} = {60}\,\mathrm{\mu} \mathrm{m}$ on both connected waveguides. The chosen length ensures an adiabatic transition between the different mode fields in the two waveguides and thus improves the coupling, see Figs. 3(c) and 3(d). On the SiP waveguide side, the linear taper converts the initial $w_{\mathrm {Si}} \times h_{\mathrm {Si}} = {500}\,\mathrm{nm} \times {220}\,\mathrm{nm}$ cross section to a $w_{\mathrm {Si,taper}} \times h_{\mathrm {Si}} = {130}\,\mathrm{nm} \times {220}\,\mathrm{nm}$ cross section at the taper tip, see Figs. 3(c) and 3(d). On the freeform waveguide side, the initial $w_{\text {WG}} \times h_{\text {WG}} = {2}\,\mathrm{\mu} \mathrm{m} \times {1.8}\,\mathrm{\mu} \mathrm{m}$ cross section is linearly tapered to a $w_{\mathrm {WG,taper}} \times h_{\mathrm {WG,taper}} = {0.76}\,\mathrm{\mu} \mathrm{m} \times {1.8}\,\mathrm{\mu} \mathrm{m}$ cross section. All freeform waveguide trajectories are plane curves with an apex height $h$ swept between $h_{\mathrm {min}} = {6.2}\,\mathrm{\mu} \mathrm{m}$ and $h_{\mathrm {max}} = {16.2}\,\mathrm{\mu} \mathrm{m}$ with a step of about 670 nm. For these experiments, we designed the beginning and the ending of each freeform waveguide trajectory (length of 22.5 µm on each side) as a circular arc, bending upwards with a radius of $R_0 = {55}\,\mathrm{\mu} \mathrm{m}$, while the central freeform waveguide trajectory part (55 µm length) is variable, generated by a parameterized B-spline, see Figs. 3(e) and 3(g). For small values of $h$, there are two sharp bends with radius $R_1$ at the two connections between the central freeform waveguide section and the two parts with the constant bend radius, see Fig. 3(g). This bend radius $R_1$ becomes larger with increasing apex height $h$, such that the loss contribution of these bends decreases. For large values of $h$, there is a sharp bend with a radius $R_2$ at the trajectory apex, which again increases the overall loss. It is therefore expected that the transmission increases with $h$ near $h_{\mathrm {min}}$ and decreases with $h$ near $h_{\mathrm {max}}$, with a maximum between the two extreme values of $h$.

To measure the transmission loss of a 3D printed freeform waveguide, continuous-wave (CW) light from a laser source at a vacuum wavelength of $\lambda = {1550}\,\mathrm{nm}$ is launched to the SiP input waveguide through a standard single-mode fiber (SMF) and a GC, Figs. 3(a) and 3(b). The light propagates through the 3D-printed freeform waveguide, and is finally coupled out through another SiP waveguide and probed by another SMF through the corresponding GC. The transmitted optical power is measured by an optical power meter (OPM). In order to exclude the coupling losses between the chip and the fibers, we measure the transmission loss of a reference structure (not shown) comprising two GC connected by a short SiP waveguide, and we refer all other transmission measurements to this value. The corresponding result is plotted on a logarithmic scale in Fig. 3(f) (marked Exp., black). To the best of our knowledge, these measurements represent the first experimental study of trajectory-dependent losses of 3D-printed freeform waveguides. The small variations in transmission for adjacent heights demonstrate the precision of the 3D printing system and the reproducibility of the overall fabrication process. Note that the cross section of 3D-printed freeform waveguide cores ($2\,\mathrm{\mu} \mathrm{m} \times 1.8\,\mathrm{\mu} \mathrm{m}$) in combination with the index difference between the core ($n_{\mathrm {core}} \approx 1.53$) and the cladding ($n_{\mathrm {cladding}} \approx 1.36$) permits propagation not only of the fundamental mode, but also of some higher-order modes. The local transmission minimum at an apex height of $h \approx {10}\,\mathrm{\mu} \mathrm{m}$, see Fig. 3(f), stems from the excitation of higher-order modes in the 3D freeform waveguide and from the resulting multimode interference at the transition to the single-mode SiP waveguide, see also Fig. 4.

 

Fig. 4. Simulated normalized magnitude $|\underline {\mathbf {E}}|$ of the complex electric field vectors $\underline {\mathbf {E}}$ for the freeform waveguide with apex height $h={9.54}\,\mathrm{\mu} \mathrm{m}$. The complex electric field $\underline {\mathbf {E}}$ is obtained from the time-domain simulation results through a Fourier transform at the target frequency, corresponding to a vacuum wavelength of $\lambda = {1550}\,\mathrm{nm}$. We compare results in the virtual $\left (u,v,s\right )$-space and in the original $\left (x,y,z\right )$-space for a TE polarized light, having a dominant electric field oriented along $x$ in $\left (x,y,z\right )$-space and along $u$ in $\left (u,v,s\right )$-space. The freeform waveguide has a plane trajectory in the $\left (y,z\right )$-plane ($x = 0$). The white contour lines in (a)–(d) are added for a better visualization of the freeform waveguide core. (a) Field distribution in $\left (u,v,s\right )$-space in plane $u = 0$, which corresponds to plane $x = 0$ in $\left (x,y,z\right )$-space. The freeform waveguide in $\left (u,v,s\right )$-space is straight, which allows to reduce the computational domain to a minimum-size rectangular volume. (b) Field distribution in $\left (x,y,z\right )$-space in the plane $x = 0$ obtained by applying the inverse space transformation $\left (x,y,z\right )^{\textrm{T}} = \mathbf {f}^{-1}\left (u,v,s\right )$ to the distribution shown in (a). (c) Field distribution in $\left (x,y,z\right )$-space in plane $x = 0$, as obtained from a reference simulation in $\left (x,y,z\right )$-space. The rectangular computational volume in $\left (x,y,z\right )$-space is not optimal since it comprises much space far away from the freeform waveguide, where the field is close to zero. The field distribution shows a good match to the distribution from (b). (d) Relative deviation between the field distributions obtained from the TO simulation in (b), and the reference simulation in (c), calculated by normalizing the magnitude of the difference of the respective electric field magnitudes to the maximum of the field magnitude found for the reference simulation. Referring to Fig. 3(f) and the field in (a)–(c), we see that multimode interference effects could explain the local minimum in transmission for apex heights $h$ around ${10}\,\mathrm{\mu} \mathrm{m}$.

Download Full Size | PDF

We then simulate the freeform waveguide transmission losses with the FIT time domain solver of CST MWS using the TO approach (in the virtual $\left (u,v,s\right )$-space as described in Section 3), and with the conventional approach (in the original $\left (x,y,z\right )$-space) as a reference. We set up the ranges of $u$ and $v$ as $u\in \left [-{2.5}\,\mathrm{\mu} \mathrm{m},{2.5}\,\mathrm{\mu} \mathrm{m}\right ]$ and $v\in \left [-{4}\,\mathrm{\mu} \mathrm{m},{4}\,\mathrm{\mu} \mathrm{m}\right ]$ in the TO-based simulations, while in the reference simulations these ranges correspond to $x\in \left [-{2.5}\,\mathrm{\mu} \mathrm{m},{2.5}\,\mathrm{\mu} \mathrm{m}\right ]$ and $y\in \left [-{4}\,\mathrm{\mu} \mathrm{m},h + {4}\,\mathrm{\mu} \mathrm{m}\right ]$. Since we are predominantly interested in the impact of the waveguide trajectory on the transmission behavior, these simulations do not take into account the coupling of the SiP waveguide to the 3D-printed freeform waveguide as detailed in Figs. 3(c) and 3(d). Instead, we assume the freeform waveguide to be embedded into a homogeneous cladding material, and we cut it at the end of linear taper of the SiP waveguide, as indicated by a dotted red line in Fig. 3(c). The tapered structure to the right of this line is then replaced by a straight waveguide section of length $l_\text {s}={1}\,\mathrm{\mu} \mathrm{m}$, in which we define the ports for the CST MWS simulation. Since a straight waveguide in $\left (x,y,z\right )$-space maps to an identical straight waveguide in $\left (u,v,s\right )$-space without any change of material properties, the modes of the straight waveguides in both spaces are identical. The length of the computational domain in the virtual $\left (u,v,s\right )$-space is thus $s_{\text {tot}} + 2l_\text {s}$, where $s_{\text {tot}}$ denotes the total length of the freeform waveguide trajectory. For the simulation in the original $\left (x,y,z\right )$-space, the length of the computational domain is only $\left (l_\text {g}+2l_\text {s}\right )$. All simulations were done with the same settings and on the same simulation machine. The discretization of the freeform waveguide into bricks for performing the TO simulations was done as explained in Section 3, with 28 bricks per slice as illustrated in Fig. 2. The total number of bricks for different freeform waveguides ranges from $1288$ to $3080$. TO-simulated transmission values in Fig. 3(f) are displayed in red, the results of the reference simulations by the conventional approach are displayed in blue. The curves obtained by the TO and the reference simulations are nearly identical, with minimal differences that can be explained by the discretization of the freeform waveguide model into finite bricks with constant material properties. We analyzed these differences and find that deviations of the transmissions predicted by both techniques range from $-0.5$ dB to $+0.6$ dB, with the biggest differences occurring at the extreme values $h_{\text {min}}$ and $h_{\text {max}}$ of the apex height $h$ where the waveguide curvature is strongest. Both curves, similarly to the experimentally obtained curve, exhibit a local minimum near $h = {10}\,\mathrm{\mu} \mathrm{m}$, and agree qualitatively very well with the measurement. The deviations between measurements and simulations can be explained by the fact that the latter do not account for the coupling of the SiP waveguide to the 3D-printed freeform waveguide as detailed in Figs. 3(c) and 3(d). In addition, the higher measured transmission for $h \approx {6}\,\mathrm{\mu} \mathrm{m}$ can be explained by a possible shrinking of the freeform waveguides during the development process, which smoothens the two sharp bends designed to have a curvature radius $R_1 \approx {4.3}\,\mathrm{\mu} \mathrm{m}$, at the connections of the center freeform waveguide section to the initial and the final sections having a constant bend radius $R_0 = {55}\,\mathrm{\mu} \mathrm{m}$, see Fig. 3(g).

We also simulate the transmission losses using the FMA method [52], which subdivides the waveguide into sections of constant curvature and calculates the transmission losses based on the propagation of the fundamental modes in these waveguide segments — the results are displayed in green in Fig. 3(f). The transmission obtained by the FMA method, does not show the local minimum near $h = {10}\,\mathrm{\mu} \mathrm{m}$, which we attribute to the fact that the FMA method disregards possible excitation of higher-order modes, such that it cannot take into account any multi-mode interference effects. Still, the optimum configurations with the least losses around a height of ${13}\,\mathrm{\mu} \mathrm{m}$ are reasonably well reproduced. Within certain limits, the FMA can therefore be used for a real-time trajectory optimization, which is important for, e.g., photonic wire bonds [12–14,16] or waveguide overpass fabrication [15] in an industrial setting.

For comparing the electric field distribution obtained by the TO and by the conventional approach, we plot the magnitude of the complex electric field vectors in the plane $u = 0$ of the virtual $\left (u,v,s\right )$-space, Fig. 4(a), and in the corresponding plane $x = 0$ of the original $\left (x,y,z\right )$-space after back-transformation, Fig. 4(b). The apex height of the depicted waveguide amounts to $h={9.54}\,\mathrm{\mu} \mathrm{m}$, which corresponds to the local minimum of the simulated transmission, caused by higher-order mode excitation. Figure 4(c) depicts the results of a CST MWS reference simulation in $\left (x,y,z\right )$-space, and Fig. 4(d) displays the magnitude of the relative deviation between the field distributions obtained from the TO simulation, Fig. 4(b), and the reference simulation, Fig. 4(c), normalized to the maximum of the field magnitude found for the reference simulation. It can be seen that the two field distributions match well. The biggest differences occur in the regions with a small bend radius — these differences are attributed to the discretization of the waveguide into bricks in the transformed space, see Section 3.

5. Computational complexity

To quantify the advantages of the TO approach, we compare the associated computational effort to that of the conventional approach. The results are displayed in Fig. 5. Figure 5(a) shows the ratios of the volumes of the computational domains for different apex heights $h$ of the freeform waveguide, see Fig. 3(b), as well as ratio of the numbers of grid cells. Both ratios indicate a reduction of the computational effort by more than a factor of two when using the TO technique. This ratio increases nearly linearly with the apex height $h$, because the computational volume of the conventional approach increases more strongly with $h$ than the one of the TO approach. Note that, for the TO-based simulations, the volume depends predominantly on the total arc length $s_{\text {tot}}$ of the freeform waveguide, which increases only slightly with $h$ as long as $h\ll s_{\text {tot}}$. The ratio of the number of grid cells for the reference simulation and the TO-simulation follows about the same proportionality as the ratio of the volumes of the computational domains. Note that the ratio of the grid cell numbers is even slightly larger than the ratio of the computational-domain volumes. We attribute this finding to the fact that, in virtual space the side walls of the straight waveguides are perfectly aligned to the rectilinear grid. Local mesh refinement as needed to accurately represent the curved waveguide surfaces in real space is hence unnecessary for the TO simulation.

 

Fig. 5. Comparison of the computational effort for the TO technique and for conventional reference simulations of freeform waveguides with different apex heights $h$ of the trajectories, see Fig. 3(b). We compare the computational volume and the number of grid cells, the time steps, and the overall simulation time of the simulations in the transformed $\left (u,v,s\right )$-space (TO) to that of the direct reference simulations in the original $\left (x,y,z\right )$-space (Ref.). (a) Ratio of the computational domain volumes and of the corresponding numbers of grid cells for the reference simulations (Ref.) and the associated TO simulations (TO). Both curves follow a (nearly) straight line as a function of $h$. (b) Comparison of time steps, determined by the smallest grid cell size and by the material properties. In the conventional reference simulations, smaller grid cell sizes $\Delta x$, $\Delta y$, and $\Delta z$, are needed to correctly represent the curved surfaces of the freeform waveguides, thus leading to smaller time steps according to Eq. (7). The large variations of the time steps in the TO-based simulations originate from different material properties. Specifically, in case of freeform waveguides with sharp bends, the tensors of material properties in the transformed $\left (u,v,s\right )$-space may have elements with values close to zero, see Appendix C, which leads to large phase velocity and thus small time steps in the TO-based model. (c) Comparison of total simulation times. Both the different numbers of grid cells and different time steps contribute to different simulation times. Overall, for the structures simulated here, the TO approach is 3–6 times more efficient than the conventional simulations, with significant potential for further improvement. Note that the trajectories here are plane curves in the $(y,z)$-plane that start and end at the same height $y$. The advantages of the TO approach related to the volume reduction of the computational domain and thus to the total simulation time reduction, would be even more pronounced for waveguides with non-plane trajectories as, e.g., shown in Fig. 1.

Download Full Size | PDF

Another parameter that influences the total simulation time is the time step used in the FIT simulation, for which the Courant-Friedrichs-Lewy stability condition for solving partial differential equations [26,54] dictates an upper limit,

The waveguide trajectories in the presented examples are still rather simple, consisting of plane curves in the $(y,z)$-plane that start and end at the same height $y$. For waveguides with non-plane trajectories, see, e.g., Fig. 1, the reduction of the computational volume through the TO approach will be even more pronounced, but the numerical treatment also becomes more complicated. Specifically, for waveguides made from isotropic materials and having plane trajectories, the tensors $\boldsymbol{\mathrm{\epsilon}}^{\prime }$ and $\boldsymbol{\mathrm{\mu}}^{\prime }$ maintain their diagonal shape in $\left (u,v,s\right )$-space, see Appendix B, while waveguides with non-plane trajectories require the consideration of off-diagonal elements in the transformed material tensors $\boldsymbol{\mathrm{\epsilon}}^{\prime }$ and $\boldsymbol{\mathrm{\mu}}^{\prime }$. A more detailed discussion of the impact of non-plane trajectories on computational complexity can be found in Appendix D. For commonly used leap-frog schemes, the FDTD update equations for the electric and the magnetic field in generally anisotropic media can be written as

In this context, it should also be noted that non-plane trajectories would also lead to an increase of the computation time that is needed to generate the transformed model from the original structure, since we also need to discretize along the $u$-axis in transformed space. For a quantitative estimate, we start from the structure shown in Fig. 2 and assume, for simplicity, that the discretization along the $u$-axis also relies on 16 bricks as chosen along the $v$-axis. The total number of bricks per $s$-slice would thus increase from presently $28$ to $16 \times 16 = 256$, i.e., by approximately an order of magnitude. This would increase the time to generate the corresponding 3D models from the current value of approx. 120 s (2 min) to approx. 1200 s (20 min) in the worst case. This extra time is significantly smaller than the more than 20 000 s of the overall simulation time that was needed to treat even the rather simple waveguides with plane trajectories in real space, see Fig. 5(c). We hence conclude that the overhead needed for generating the transformed model is significantly overcompensated by the savings of the simulation time used for the Maxwell’s equations solver also for non-plane trajectories.

Note that the simulations shown in the previous sections refer to waveguides with rectangular cross sections that are invariant along the waveguide trajectory, which can be accurately represented by a rather simple rectilinear grid. In case of arbitrary cross sections that vary along the waveguide trajectory, this advantage of the TO approach might be maintained by adapting the transformation to not only map the curved trajectory into a straight one, but to also map an arbitrary waveguide cross section into a rectangular one [59] that is invariant along the propagation direction.

6. Summary

We introduced a transformation-optics (TO) approach for simulating freeform optical waveguides, which is applicable to commercially available time-domain Maxwell’s equations solvers. The method reduces the computational volume by transforming a curved freeform waveguide in the original 3D space into a straight structure with a modified permittivity and permeability profile. Furthermore, the method allows for a better alignment of grid lines and surfaces of rectangular freeform waveguides, thus avoiding local mesh cell refinement. We verify the viability of our technique by simulating freeform waveguides with plane trajectories, thereby demonstrating a significant reduction in computational effort compared to the simulations in the original 3D space while maintaining the same level of accuracy. For experimental benchmarking, we realized a series of freeform waveguides and measured their transmission losses, finding good qualitative agreement. To the best of our knowledge, these measurements represent the first experimental study of trajectory-dependent losses of 3D-printed freeform waveguides. We believe that our TO-based simulation approach has the potential to greatly facilitate design and prototyping of optical devices based on 3D freeform waveguides.

Appendix A: on the bijectivity of the coordinate transformation function

A necessary condition for the bijectivity of the coordinate transformation function $\left (u,v,s\right )^{\text {T}} = \mathbf {f}\left (x,y,z\right )$ described by its inverse function $\left (x,y,z\right )^{\text {T}} = \mathbf {f}^{-1}\left (u,v,s\right )$ in Eq. (3) is that the computational domain does not contain local center of curvature of the freeform waveguide trajectory. This can be achieved by appropriately choosing the ranges of $u$- and $v$-coordinates. As a simple illustrative example we may think of a freeform waveguide with a plane trajectory in the $\left (y,z\right )$-plane consisting of two straight sections that are connected by a $90^{\circ }$ bend with a constant bend radius $R$. The local center of curvature is point $C$ in the $\left (y,z\right )$-plane, see Fig. 6(a). As explained in Section 3, the RMF in all points of the trajectory is chosen such that vector $\mathbf {U}$, pointing out of the drawing plane, is parallel to the $x$-axis and perpendicular to the $\left (y,z\right )$-plane. As a consequence, vector $\mathbf {V}$ is parallel to the $\left (y,z\right )$-plane. Since the vector $\mathbf {V}$ is parallel to the $v$-axis in $\left (u,v,s\right )$-space, the range $v \in \left [-d,d\right ]$ determines whether the point $C$ is inside or outside the computational domain.

 

Fig. 6. Waveguide bends for illustrating a necessary condition for the bijectivity of the coordinate transformation function $\left (u,v,s\right )^{\text {T}} = \mathbf {f}\left (x,y,z\right )$, which is defined by its inverse function $\left (x,y,z\right )^{\text {T}} = \mathbf {f}^{-1}\left (u,v,s\right )$ in Eq. (3). The different computational domains are limited by the outer boundaries of the hatched areas. We consider a plane freeform waveguide trajectory (dashed lines) in the $\left (y,z\right )$-plane, consisting of two straight waveguide segments connected by a $90^{\circ }$ circular bend with a bend radius $R$. The $\left (y,z\right )$ plane in $\left (x,y,z\right )$-space (left drawings) corresponds to the $\left (u,v\right )$-plane in $\left (u,v,s\right )$-space (right drawings). The local center of curvature corresponds to the center $C$ of the $90^{\circ }$ bend. We consider a symmetric range of $v$-coordinates, $v \in \left [-d,d\right ]$. The RMF is shown only in case (a). (a) If $d<R$, the point $C$ is outside the computational domain, ensuring bijective mapping between the two spaces. (b) If $d=R$, the point $C$ is on the border of the computational domain. In this critical case, the space transformation is not anymore a bijection, since point $C$ is mapped onto a line segment $C_1 C_2$ in $\left (u,v,s\right )$-space. (c) If $d>R$, the subdomains $A$ and $B$ in $\left (x,y,z\right )$-space are mapped to multiple sub-domains $A_1, A_2, A_3$ and $B_1, B_2$, respectively, in $\left (u,v,s\right )$-space, and the one-to-one correspondence between the two spaces is further violated.

Download Full Size | PDF

Without going into mathematical details, we give a qualitative analysis of three cases: $d<R$, $d=R$, and $d>R$, Figs. 6(a)–6(c) here — a more detailed mathematical analysis of the spatial transformation of $90^{\circ }$ bends can be found in Appendix C. Dashed lines represent waveguide trajectories, the white part in the middle represents the waveguide core, and the hatched parts to both sides represent the portion of cladding that is included in the computational domain. For all three cases we provide two drawings: One for the original freeform waveguide in the $\left (y,z\right )$-plane of $\left (x,y,z\right )$-space, and one for the corresponding straight waveguide in the $\left (v,s\right )$-plane of $\left (u,v,s\right )$-space. If $d<R$, the computational domain does not include the center point $C$ of the bend, and the mapping between the two spaces is a bijection, see Fig. 6(a). In case $d=R$, point $C$ is on the border of the computational domain and is mapped to a line segment $C_1 C_2$ in $\left (u,v,s\right )$-space, the length of which is equal to the arc-length of the $90^{\circ }$ bend of the trajectory, $\overline {C_1 C_2} = R\pi /2$, see Fig. 6(b). This is the critical case when the spatial transformation function is not a bijection anymore. In case $d>R$, the one-to-one correspondence between the two spaces is further violated. Not only the point $C$ is again mapped to the line segment $C_1C_2$, but also areas marked with $A$ and $B$ in the $\left (y,z\right )$-plane are mapped into multiple areas in the $\left (v,s\right )$-plane — see Fig. 6(c). In particular, area $A$ is mapped to areas $A_1$, $A_2$, and $A_3$, while area $B$ is mapped to areas $B_1$ and $B_2$. Based on this simple example, it is clear that a necessary condition for ensuring bijectivity is that $C$ in Fig. 6 is outside of the computational domain. On the other hand, for a physically correct representation of the waveguide, the essential part of the evanescent fields at the outside of the bend must be inside the computational domain. This might be accomplished by choosing a computational domain that comprises an asymmetric range of $v$-coordinates instead of the symmetric range $\left [-d,d\right ]$ used in our example.

This simple example is also very illustrative for the general case of freeform waveguides. Bends need not necessarily be $90^{\circ }$ bends, and the bend radius can continuously change along the trajectory. In case of non-plane true 3D trajectories, the local center of the bend lies in the local osculating plane, and the same reasoning provided above can be applied. As a matter of fact, the osculating plane for our example is the $\left (y,z\right )$-plane — the only difference to 3D trajectories would be that the local vectors of the RMF are not necessarily perpendicular and parallel to the local osculating plane for each point on the trajectory.

Appendix B: tensors of $\boldsymbol{\mathrm{\epsilon}}^{\prime }$ and $\boldsymbol{\mathrm{\mu}}^{\prime }$ in $\mathbf {\left (u,v,s\right )}$-space for freeform waveguides with plane trajectories

In case of freeform waveguides with plane trajectories and isotropic material properties $\varepsilon$ and $\mu$ in the original $\left (x,y,z\right )$-space, the material properties $\boldsymbol{\mathrm{\epsilon}}^{\prime }$ and $\boldsymbol{\mathrm{\mu}}^{\prime }$ in the virtual $\left (u,v,s\right )$-space are diagonal tensors. This can be shown by assuming that, without loss of generality, the plane trajectory lies in $\left (y,z\right )$-plane, as explained in Section 3. and Appendix A, such that the vector $\mathbf {U}$ of the RMF in each trajectory point is oriented parallel to the $x$-axis, while the remaining two vectors $\mathbf {V}$ and $\mathbf {T}$ of the RMF lie in the $\left (y,z\right )$-plane, with the vector $\mathbf {T}$ forming an angle $\theta$ with the positive $z$-axis (see Fig. 7). Since the vectors $\mathbf {V}$ and $\mathbf {T}$ are parallel to the $v$- and $s$-axes in the $\left (u,v,s\right )$-space, respectively, the axes of the 2D $\left (v,s\right )$-coordinate system are rotated by the same angle $\theta$ with respect to the axes of the 2D $\left (y,z\right )$-coordinate system. This allows us to simplify the calculation of the Jacobian matrix given by Eq. (2). Since the vector $\mathbf {U}$ and the $x$-axis are parallel to each other, it follows that $\partial u/\partial x = 1$. This also implies that the partial derivatives of $u$ with respect to $y$ and $z$ must be zero, $\partial u/\partial y = 0$, and $\partial u/\partial z = 0$. Furthermore, the vectors $\mathbf {V}$ and $\mathbf {T}$ are perpendicular to the $x$-axis, which implies that the derivatives of $v$ and $s$ with respect to $x$ must be zero, too: $\partial v/\partial x = 0$ and $\partial s/\partial x = 0$. The Jacobian matrix of the coordinate transformation function thus reads

(9)$$\mathbf{J} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \frac{\partial v}{\partial y} & \frac{\partial v}{\partial z} \\ 0 & \frac{\partial s}{\partial y} & \frac{\partial s}{\partial z} \\ \end{bmatrix}\,.$$

Assuming the material properties are isotropic in $\left (x,y,z\right )$-space, Eq. (1) can be simplified,

(10)$$ \begin{aligned} &\boldsymbol{\mathrm{\epsilon}}^{\prime}(u, v, s)=\varepsilon(x, y, z) \frac{\mathbf{J} \cdot \mathbf{J}^{\mathbf{T}}}{\operatorname{det}(\mathbf{J})} \\ &\boldsymbol{\mathrm{\mu}}^{\prime}(u, v, s)=\mu(x, y, z) \frac{\mathbf{J} \cdot \mathbf{J}^{\mathbf{T}}}{\operatorname{det}(\mathbf{J})} \end{aligned} $$

where the product of the Jacobian matrix and its transposed reads

(11)$$\mathbf{J}\cdot\mathbf{J^{T}} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \left(\frac{\partial v}{\partial y}\right)^{2} + \left(\frac{\partial v}{\partial z}\right)^{2} & \left(\frac{\partial y}{\partial v}\frac{\partial y}{\partial s}\right)^{{-}1} + \left(\frac{\partial z}{\partial v}\frac{\partial z}{\partial s}\right)^{{-}1} \\ 0 & \left(\frac{\partial y}{\partial v}\frac{\partial y}{\partial s}\right)^{{-}1} + \left(\frac{\partial z}{\partial v}\frac{\partial z}{\partial s}\right)^{{-}1} & \left(\frac{\partial s}{\partial y}\right)^{2} + \left(\frac{\partial s}{\partial z}\right)^{2} \end{bmatrix}.$$

Equation (11) is a diagonal matrix if its two off-diagonal elements are equal to zero, which reduces to

(12)$$\frac{\partial y}{\partial v}\frac{\partial y}{\partial s} + \frac{\partial z}{\partial v}\frac{\partial z}{\partial s} = 0\,.$$

Since $\varepsilon \left (x,y,z\right )$, $\mu \left (x,y,z\right )$, and $\det \left (\mathbf {J}\right )$ are scalars, it follows from Eqs. (10) and (11) that Eq. (12) is a sufficient condition that ensures that $\boldsymbol{\mathrm{\epsilon}}^{\prime }\left (u,v,s\right )$ and $\boldsymbol{\mathrm{\mu}}^{\prime }\left (u,v,s\right )$ are diagonal tensors.

 

Fig. 7. Relationship between coordinates in $\left (x,y,z\right )$- and $\left (u,v,s\right )$-space in a point $A$ with coordinates $\left (y_0\left (s\right ),z_0\left (s\right )\right )$ on the trajectory for the case of a freeform waveguide with a plane trajectory. The trajectory lies in the $\left (y,z\right )$-plane, and the RMF is oriented such that the vector $\mathbf {U}$ and the associated $u$-axis are parallel to the $x$-axis in all points of the trajectory. The remaining two axes $\mathbf {V}$ and $\mathbf {T}$ of the RMF define a 2D frame that lies in the $\left (y,z\right )$-plane. The axes of the local $\left (v,s\right )$-coordinate system are parallel to the $\left (\mathbf {V},\mathbf {T}\right )$ frame and rotated by an angle $\theta$ with respect to the axes of the $\left (y,z\right )$-coordinate system.

Download Full Size | PDF

The coordinate transformation from the $\left (v,s\right )$-coordinate system to the $\left (y,z\right )$-coordinate system can be extracted from the sketches in Fig. 7,

(13)$$\begin{bmatrix} y\left(s,v\right)\\ z\left(s,v\right)\\ \end{bmatrix} = \begin{bmatrix} y_0\left(s\right)\\ z_0\left(s\right)\\ \end{bmatrix} + \begin{bmatrix} v\cos{\theta\left(s\right)}\\ v\sin{\theta\left(s\right)}\\ \end{bmatrix} ,$$

where $\left (y_0\left (s\right ),(z_0\left (s\right )\right )$ are coordinates of the point $A$ on the trajectory that is defined by arc-length coordinate $s$. Equation (12) can be expressed as a dot-product of two $(2,1)$-vectors,

(14)$$\frac{\partial}{\partial v} \begin{bmatrix} y\left(s,v\right)\\ z\left(s,v\right)\\ \end{bmatrix} \cdot \frac{\partial}{\partial s} \begin{bmatrix} y\left(s,v\right)\\ z\left(s,v\right)\\ \end{bmatrix} = 0\,.$$

Inserting Eq. (13) into Eq. (14), we obtain

(15)$$\begin{bmatrix} \cos{\theta\left(s\right)}\\ \sin{\theta\left(s\right)}\\ \end{bmatrix} \cdot \left(\mathbf{T}\left(s\right) + \begin{bmatrix} -v\sin{\theta\left(s\right)}\frac{\partial \theta\left(s\right)}{\partial s}\\ v\cos{\theta\left(s\right)}\frac{\partial \theta\left(s\right)}{\partial s}\\ \end{bmatrix} \right) = \mathbf{N}\left(s\right) \cdot \mathbf{T}\left(s\right)\left(1 + v\frac{\partial \theta\left(s\right)}{\partial s}\right) = 0,$$

where $\mathbf {N}\left (s\right )$ denotes the unit normal vector in point $A$ on the trajectory, which in case of plane trajectories is parallel to vector $\mathbf {V}\left (s\right )$. Since the dot product of the normal and the tangent vector is always equal to zero, Eq. (15) is always fulfilled, and $\boldsymbol{\mathrm{\mu}} ^{\prime }\left (u,v,s\right )$ are thus diagonal tensors.

Appendix C: time-stepping

We have already shown in Appendix A that no local center of curvature of the trajectory is allowed to be within the computational domain to maintain the bijectivity of the space transformation defined by Eq. (3). In addition, if the computational domain border is too close to a local center of curvature, some elements of the tensors $\boldsymbol{\mathrm{\epsilon}} ^{\prime }\left (u,v,s\right )$ and $\boldsymbol{\mathrm{\mu}} ^{\prime }\left (u,v,s\right )$ can assume values close to zero. The Courant-Friedrichs-Lewy stability condition, Eq. (7), then leads to small time steps of the corresponding time-domain simulation and greatly increases the overall simulation time. For illustration and similarly to Appendix A, we discuss a circular $90^{\circ }$-bend with radius $R$ in the $\left (y,z\right )$-plane, see Fig. 8. Assuming, without loss of generality, that the center of the circular arc is in the origin of the $\left (x,y,z\right )$ coordinate system, the trajectory is given by

(16)$$\begin{aligned} & x_0 = 0, \\ & y_0 ={-}R\cos{\left(t\right)}, \\ & z_0 = R\sin{\left(t\right)}, \end{aligned}$$

with $t\in \left [0,\pi /2\right ]$. The coordinate $s = R t$ is the arc length of the trajectory. Since the $u$-axis is parallel to the $x$-axis, and the $v$-axis is in the plane of the trajectory, we can write the one-to-one correspondence between the $\left (x,y,z\right )$- and $\left (u,v,s\right )$-spaces, see Fig. 8,

(17)$$\begin{aligned} & x\left(u,v,s\right) = u, \\ & y\left(u,v,s\right) ={-}\left(R-v\right)\cos{\left(\frac{s}{R}\right)}, \\ & z\left(u,v,s\right) = \left(R-v\right)\sin{\left(\frac{s}{R}\right)}. \end{aligned}$$

The Jacobian $\mathbf {J}$ of the function $\left (u,v,s\right )^{\text {T}} = \mathbf {f}\left (x,y,z\right )$ can be found as the inverse of the Jacobian of the function $\left (x,y,z\right )^{\text {T}} = \mathbf {f}^{-1}\left (u,v,s\right )$,

(18)$$\mathbf{J} = \begin{bmatrix} \frac{\partial x}{\partial u} & \frac{\partial x}{\partial v} & \frac{\partial x}{\partial s} \\ \frac{\partial y}{\partial u} & \frac{\partial y}{\partial v} & \frac{\partial y}{\partial s} \\ \frac{\partial z}{\partial u} & \frac{\partial z}{\partial v} & \frac{\partial z}{\partial s} \\ \end{bmatrix}^{{-}1} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos{\left(\frac{s}{R}\right)} & -\sin{\left(\frac{s}{R}\right)} \\ 0 & \frac{R}{R-v}\sin{\left(\frac{s}{R}\right)} & \frac{R}{R-v}\cos{\left(\frac{s}{R}\right)}\\ \end{bmatrix}.$$

Inserting this result into Eq. (10) leads to

(19)$$\boldsymbol{\mathrm{\epsilon}}^{\prime}\left(u,v,s\right) = \varepsilon\left(x,y,z\right) \begin{bmatrix} \frac{R-v}{R} & 0 & 0 \\ 0 & \frac{R-v}{R} & 0 \\ 0 & 0 & \frac{R}{R-v}\\ \end{bmatrix}= \begin{bmatrix} \varepsilon_{1,1} & 0 & 0 \\ 0 & \varepsilon_{2,2} & 0 \\ 0 & 0 & \varepsilon_{3,3}\\ \end{bmatrix}.$$

A similar result can be obtained for $\boldsymbol{\mathrm{\mu}}^{\prime }$. Hence, when $v$ approaches $R$, $\varepsilon _{1,1}$ and $\varepsilon _{2,2}$ ($\mu _{1,1}$ and $\mu _{2,2}$) tend to 0, causing the maximum phase velocity to tend to infinity. According to the Courant-Friedrichs-Lewy stability condition Eq. (7), the maximal time step then tends to zero, and the total simulation time approaches infinity. On the outer side of the bend, for $v<0$, tensor elements $\varepsilon _{3,3}$ and $\mu _{3,3}$ are more problematic, since they tend to zero for $v \rightarrow -\infty$. This is, however, not of practical relevance since the space far away from the trajectory is commonly not of interest for modeling of 3D freeform waveguides.

 

Fig. 8. Illustration of the space transformation for a $90^{\circ }$-bend. (a) $90^{\circ }$-bend waveguide in the $\left (y,z\right )$-plane of $\left (x,y,z\right )$-space. (b) Same waveguide in the $\left (v,s\right )$-plane of $\left (u,v,s\right )$-space. In $\left (u,v,s\right )$-space, the $90^{\circ }$-bend is straightened. It is divided into squares by taking equidistant divisions along the $s$-axis (slicing) and along the $v$-axis. The squares in $\left (u,v,s\right )$-space correspond to sections bounded by two line segments and two concentric arcs in $\left (x,y,z\right )$-space. The areas of these sections differ along the radial coordinate and are smaller on the inner side and larger on the outer side of the bend. Areas close to the origin tend to zero, and mapping these sections to finite-size squares in $\left (u,v,s\right )$-space causes two elements of tensors of material properties to tend to zero, see Eq. (19). This leads to a small time step of the corresponding FDTD simulation due to the Courant-Friedrichs-Lewy stability condition, Eq. (7), and thus to long simulation times. Areas of sections far from the origin tend to infinity, and mapping infinite sections to finite size squares in $\left (u,v,s\right )$-space also causes one of the tensor elements of material properties to tend to zero. These cases are, however, not critical, since the space far away from the trajectory is commonly not of interest for modeling of 3D freeform waveguides.

Download Full Size | PDF

Appendix D: impact of non-plane trajectories on computational complexity

In general, freeform waveguide trajectories are non-plane, which leads to non-diagonal tensors $\boldsymbol{\mathrm{\mu}}^{\prime }$ and $\boldsymbol{\mathrm{\epsilon}}^{\prime }$ in the transformed $\left (u,v,s\right )$-space. This increases the computational complexity, as more additions and multiplications need to be done in each time step of the Maxwell’s equations solver compared to the case with diagonal tensors $\boldsymbol{\mathrm{\mu}}^{\prime }$ and $\boldsymbol{\mathrm{\epsilon}}^{\prime }$. As an example, we will consider an FDTD solver, which is based on the leap-frog algorithm, where in each time step, the electric field $\underline {\mathbf {E}}^{\prime }$ and the magnetic field $\underline {\mathbf {H}}^{\prime }$ are updated successively. The update relation is given by Eq. (8), that we repeat here for convenience:

(20)$$\begin{aligned} & \underline{\mathbf{H}}^{\prime}(t + \Delta t/2) = \underline{\mathbf{H}}^{\prime}(t - \Delta t/2) - \Delta t \left(\boldsymbol{\mathrm{\mu}}^{\prime}\left(u,v,s\right)\right)^{{-}1} \cdot \left(\mathbf{\nabla}\times{\underline{\mathbf{E}}^{\prime}(t)}\right),\\ & \underline{\mathbf{E}}^{\prime}(t + \Delta t) = \underline{\mathbf{E}}^{\prime}(t) + \Delta t \left(\boldsymbol{\mathrm{\epsilon}}^{\prime}\left(u,v,s\right)\right)^{{-}1} \cdot \left(\mathbf{\nabla}\times{\underline{\mathbf{H}}^{\prime}(t + \Delta t/2)}\right).\\ \end{aligned}$$

For updating the vector of the electric (magnetic) field in Eq. (20), we need to calculate a curl of the magnetic (electric) field vector, multiply the $\left (3,3\right )$ matrix $\Delta t \left (\boldsymbol{\mathrm{\epsilon}} ^{\prime }\right )^{-1}$ ($\Delta t \left (\boldsymbol{\mathrm{\mu}} ^{\prime }\right )^{-1}$) by the calculated curl, and add the result to the vector of the electric (magnetic) field in the previous time point. To find a curl of the $\left (3,1\right )$ vector, we need to calculate six partial derivatives. This requires nine additions — six additions ($6A$) to calculate the finite differences between field components in spatial grid points and three more additions ($3A$) for calculating the difference between the two spatial derivatives in each component of the curl — and six multiplications ($6M$) for dividing the finite differences of the fields by the corresponding size of the spatial grid cell. In case of plane trajectories and isotropic materials, the transformed permittivity and permeability tensors $\boldsymbol{\mathrm{\epsilon}} ^{\prime }$ and $\boldsymbol{\mathrm{\mu}} ^{\prime }$ are diagonal, and the same applies to the $\left (3,3\right )$ matrices $\Delta t \left (\boldsymbol{\mathrm{\epsilon}} ^{\prime }\right )^{-1}$ and $\Delta t \left (\boldsymbol{\mathrm{\mu}} ^{\prime }\right )^{-1}$. In this case, multiplication of these matrices by the curl vector requires three additional scalar multiplications ($3M$). In contrast to that, non-plane trajectories lead to non-diagonal tensors, for which a full vector-matrix multiplication with nine multiplications and six additions is needed ($6A+9M$). Finally, in both cases, adding the resulting vector to the field vector in the previous time point requires three more additions ($3A$). Hence, updating the electric or the magnetic field in one spatial grid point requires twelve additions and nine multiplications $12A + 9M$ in case of plane trajectories and diagonal permittivity and permeability tensors, and eighteen additions and fifteen multiplications $18A + 15M$ in case of non-diagonal tensors. The ratio of necessary arithmetic operations is therefore $\frac {18A + 15M}{12A + 9M}<\frac {5}{3}$. Based on this rather simple model for computational effort, we estimate that moving from plane trajectories to non-plane trajectories increases the computational effort by a factor of roughly 2. Note that additional complications arise from the fact that the standard Yee grid cannot be used for non-diagonal permittivity and permeability tensors [55]. On the other hand, the reduction of simulation time for non-plane trajectories should even exceed the factor of 3–6 that we found for the case of a plane trajectory, see Fig. 5(c). We thus expect that the overall computational effort should still be greatly reduced by using the proposed TO approach

Funding

Deutsche Forschungsgemeinschaft (258734477, 390761711, EXC-2082/1); Erasmus Mundus Joint Doctorate Program Europhotonics (159224-1-2009-1-FR-ERA MUNDUS-EMJD); European Research Council (773248); Bundesministerium für Bildung und Forschung (16ES0948); Karlsruhe School of Optics and Photonics (KSOP); Alfried Krupp von Bohlen und Halbach-Stiftung.

Acknowledgments

This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the Collaborative Research Center “Wave Phenomena: Analysis and Numerics” (CRC $1173$), Project C4 $(\#258734477)$ as well as under Germany’s Excellence Strategy via the Excellence Cluster 3D Matter Made to Order (3DMM2O, EXC-2082/1, #390761711), by the European Research Council (ERC) Consolidator Grant “TeraSHAPE” (#773248), by the Bundesministerium für Bildung und Forschung (BMBF) via the project DiFeMiS (#16ES0948), by the Karlsruhe School of Optics and Photonics (KSOP), and by the Alfried Krupp von Bohlen und Halbach Foundation. A.N. was supported by the Erasmus Mundus doctorate programme Europhotonics (grant number 159224-1-2009-1-FR-ERA MUNDUS-EMJD).

Disclosures

C.K. is a co-founder and a shareholder of Vanguard Photonics GmbH and Vanguard Automation GmbH, and A.N. is a former employee of Vanguard Automation GmbH. Vanguard Photonics GmbH and Vanguard Automation GmbH are start-up companies engaged in commercializing 3D-printing technologies in the field of photonic integration. M.B. is employed by Nanoscribe GmbH, a company developing 3D printing tools.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. P. Dong, Y.-K. Chen, G.-H. Duan, and D. T. Neilson, “Silicon photonic devices and integrated circuits,” Nanophotonics 3(4-5), 215–228 (2014). [CrossRef]  

2. T. Harter, S. Muehlbrandt, S. Ummethala, A. Schmid, S. Nellen, L. Hahn, W. Freude, and C. Koos, “Silicon–plasmonic integrated circuits for terahertz signal generation and coherent detection,” Nat. Photonics 12(10), 625–633 (2018). [CrossRef]  

3. S. Ummethala, T. Harter, K. Koehnle, Z. Li, S. Muehlbrandt, Y. Kutuvantavida, J. Kemal, P. Marin-Palomo, J. Schaefer, A. Tessmann, S. K. Garlapati, A. Bacher, L. Hahn, M. Walther, T. Zwick, S. Randel, W. Freude, and C. Koos, “THz-to-optical conversion in wireless communications using an ultra-broadband plasmonic modulator,” Nat. Photonics 13(8), 519–524 (2019). [CrossRef]  

4. Y. Shen, N. C. Harris, S. Skirlo, M. Prabhu, T. Baehr-Jones, M. Hochberg, X. Sun, S. Zhao, H. Larochelle, D. Englund, and M. Soljačić, “Deep learning with coherent nanophotonic circuits,” Nat. Photonics 11(7), 441–446 (2017). [CrossRef]  

5. S. Schneider, M. Lauermann, P.-I. Dietrich, C. Weimann, W. Freude, and C. Koos, “Optical coherence tomography system mass-producible on a silicon photonic chip,” Opt. Express 24(2), 1573–1586 (2016). [CrossRef]  

6. P. Trocha, M. Karpov, D. Ganin, M. H. P. Pfeiffer, A. Kordts, S. Wolf, J. Krockenberger, P. Marin-Palomo, C. Weimann, S. Randel, W. Freude, T. J. Kippenberg, and C. Koos, “Ultrafast optical ranging using microresonator soliton frequency combs,” Science 359(6378), 887–891 (2018). [CrossRef]  

7. J. Sun, E. Timurdogan, A. Yaacobi, E. S. Hosseini, and M. R. Watts, “Large-scale nanophotonic phased array,” Nature 493(7431), 195–199 (2013). [CrossRef]  

8. C. Ranacher, C. Consani, A. Tortschanoff, R. Jannesari, M. Bergmeister, T. Grille, and B. Jakoby, “Mid-infrared absorption gas sensing using a silicon strip waveguide,” Sensors Actuators A: Phys. 277, 117–123 (2018). [CrossRef]  

9. A. F. Gavela, D. G. García, J. Ramirez, and L. Lechuga, “Last advances in silicon-based optical biosensors,” Sensors 16(3), 285 (2016). [CrossRef]  

10. P. Steglich, M. Hülsemann, B. Dietzel, and A. Mai, “Optical biosensors based on silicon-on-insulator ring resonators: A review,” Molecules 24(3), 519 (2019). [CrossRef]  

11. D. Kohler, G. Schindler, L. Hahn, J. Milvich, A. Hofmann, K. Länge, W. Freude, and C. Koos, “Biophotonic sensors with integrated Si₃N₄-organic hybrid (SiNOH) lasers for point-of-care diagnostics,” Light: Sci. Appl. 10(1), 64 (2021). [CrossRef]  

12. N. Lindenmann, S. Dottermusch, M. L. Goedecke, T. Hoose, M. R. Billah, T. P. Onanuga, A. Hofmann, W. Freude, and C. Koos, “Connecting silicon photonic circuits to multicore fibers by photonic wire bonding,” J. Lightwave Technol. 33(4), 755–760 (2015). [CrossRef]  

13. M. R. Billah, M. Blaicher, T. Hoose, P.-I. Dietrich, P. Marin-Palomo, N. Lindenmann, A. Nesic, A. Hofmann, U. Troppenz, M. Moehrle, S. Randel, W. Freude, and C. Koos, “Hybrid integration of silicon photonics circuits and InP lasers by photonic wire bonding,” Optica 5(7), 876–883 (2018). [CrossRef]  

14. M. Blaicher, M. R. Billah, J. Kemal, T. Hoose, P. Marin-Palomo, A. Hofmann, Y. Kutuvantavida, C. Kieninger, P.-I. Dietrich, M. Lauermann, S. Wolf, U. Troppenz, M. Moehrle, F. Merget, S. Skacel, J. Witzens, S. Randel, W. Freude, and C. Koos, “Hybrid multi-chip assembly of optical communication engines by in situ 3D nano-lithography,” Light: Sci. Appl. 9(1), 71 (2020). [CrossRef]  

15. A. Nesic, M. Blaicher, T. Hoose, A. Hofmann, M. Lauermann, Y. Kutuvantavida, M. Nöllenburg, S. Randel, W. Freude, and C. Koos, “Photonic-integrated circuits with non-planar topologies realized by 3D-printed waveguide overpasses,” Opt. Express 27(12), 17402–17425 (2019). [CrossRef]  

16. T. Hoose, M. Blaicher, J. N. Kemal, H. Zwickel, M. R. Billah, P.-I. Dietrich, A. Hofmann, W. Freude, S. Randel, and C. Koos, “Hardwire-configurable photonic integrated circuits enabled by 3D nano-printing,” arXiv:1912.09942 (2019).

17. H. Gehring, M. Blaicher, T. Grottke, and W. H. P. Pernice, “Reconfigurable nanophotonic circuitry enabled by direct-laser-writing,” IEEE J. Sel. Top. Quantum Electron. 26(5), 1–5 (2020). [CrossRef]  

18. J. Moughames, X. Porte, M. Thiel, G. Ulliac, L. Larger, M. Jacquot, M. Kadic, and D. Brunner, “Three-dimensional waveguide interconnects for scalable integration of photonic neural networks,” Optica 7(6), 640–646 (2020). [CrossRef]  

19. P. Gašo, D. Pudiš, D. Seyringer, A. Kuzma, L. Gajdošová, T. Mizera, and M. Goraus, “3D polymer based 1 × 4 beam splitter,” J. Lightwave Technol. 39(1), 154–161 (2021). [CrossRef]  

20. A. Nesic, M. Blaicher, P. Marin-Palomo, C. Füllner, S. Randel, W. Freude, and C. Koos, “Ultra-broadband polarization beam splitter and rotator based on 3D-printed waveguides,” arXiv:2107.12282 (2021).

21. Ansys-Lumerical, https://www.lumerical.com. Accessed: 2022-04-04.

22. Synopsis-RSoft Photonic Device Tools, https://www.synopsys.com/photonic-solutions/rsoft-photonic-device-tools.html. Accessed: 2022-04-04.

23. Optiwave Photonic Software, https://optiwave.com/. Accessed: 2022-04-04.

24. MIT Electromagnetic Equation Propagation — MEEP, https://meep.readthedocs.io/en/latest/. Accessed: 2022-04-04.

25. M. Clemens and T. Weil, “Discrete electromagnetism with the finite integration technique,” Prog. Electromagn. Res. 32, 65–87 (2001). [CrossRef]  

26. Computer Simulation Technology, Understanding Time Domain Meshing in CST Microwave Studio® (2010).

27. L. S. Dolin, “To the possibility of comparison of three-dimensional electromagnetic systems with nonuniform anisotropic filling,” Izv. VUZov, Radiofizika 4(5), 964–967 (1961).

28. E. Post, Formal Structure of Electromagnetics: General Covariance and Electromagnetics, North-Holland Series in Physics (North-Holland Publishing Company, 1962).

29. I. Lindell, Methods for Electromagnetic Field Analysis (Clarendon Press, 1992).

30. A. J. Ward and J. B. Pendry, “Refraction and geometry in Maxwell’s equations,” J. Mod. Opt. 43(4), 773–793 (1996). [CrossRef]  

31. M. Rahm, D. A. Roberts, J. B. Pendry, and D. R. Smith, “Transformation-optical design of adaptive beam bends and beam expanders,” Opt. Express 16(15), 11555–11567 (2008). [CrossRef]  

32. D.-H. Kwon and D. H. Werner, “Polarization splitter and polarization rotator designs based on transformation optics,” Opt. Express 16(23), 18731–18738 (2008). [CrossRef]  

33. D.-H. Kwon and D. H. Werner, “Flat focusing lens designs having minimized reflection based on coordinate transformation techniques,” Opt. Express 17(10), 7807–7817 (2009). [CrossRef]  

34. L. H. Gabrielli, D. Liu, S. G. Johnson, and M. Lipson, “On-chip transformation optics for multimode waveguide bends,” Nat. Commun. 3(1), 1217 (2012). [CrossRef]  

35. J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic fields,” Science 312(5781), 1780–1782 (2006). [CrossRef]  

36. D. Schurig, J. J. Mock, B. J. Justice, S. A. Cummer, J. B. Pendry, A. F. Starr, and D. R. Smith, “Metamaterial electromagnetic cloak at microwave frequencies,” Science 314(5801), 977–980 (2006). [CrossRef]  

37. U. Leonhardt and T. Tyc, “Broadband invisibility by non-Euclidean cloaking,” Science 323(5910), 110–112 (2009). [CrossRef]  

38. Y. Lai, J. Ng, H. Chen, D. Han, J. Xiao, Z.-Q. Zhang, and C. T. Chan, “Illusion optics: The optical transformation of an object into another object,” Phys. Rev. Lett. 102(25), 253902 (2009). [CrossRef]  

39. O. Ozgun and M. Kuzuoglu, “Form invariance of Maxwell’s equations: The pathway to novel metamaterial specifications for electromagnetic reshaping,” IEEE Trans. Antennas Propag. 52(3), 51–65 (2010). [CrossRef]  

40. E. A. Berry, J. Gutierrez, and R. Rumpf, “Design and simulation of arbitrarily-shaped transformation optic devices using a simple finite-difference method,” Prog. Electromagn. Res. B 68, 1–16 (2016). [CrossRef]  

41. M. Kadic, S. Guenneau, S. Enoch, P. A. Huidobro, L. Martín-Moreno, F. J. García-Vidal, J. Renger, and R. Quidant, “Transformation plasmonics,” Nanophotonics 1(1), 51–64 (2012). [CrossRef]  

42. M. Kadic, T. Bückmann, R. Schittny, and M. Wegener, “Experiments on cloaking in optics, thermodynamics and mechanics,” Phil. Trans. R. Soc. A. 373(2049), 20140357 (2015). [CrossRef]  

43. M. McCall, J. B. Pendry, V. Galdi, Y. Lai, S. A. R. Horsley, J. Li, J. Zhu, R. C. Mitchell-Thomas, O. Quevedo-Teruel, P. Tassin, V. Ginis, E. Martini, G. Minatti, S. Maci, M. Ebrahimpouri, Y. Hao, P. Kinsler, J. Gratus, J. M. Lukens, A. M. Weiner, U. Leonhardt, I. I. Smolyaninov, V. N. Smolyaninova, R. T. Thompson, M. Wegener, M. Kadic, and S. A. Cummer, “Roadmap on transformation optics,” J. Opt. 20(6), 063001 (2018). [CrossRef]  

44. R. Joseph and A. Taflove, “FDTD Maxwell's equations models for nonlinear electrodynamics and optics,” IEEE Trans. Antennas Propag. 45(3), 364–374 (1997). [CrossRef]  

45. R. L. Bishop, “There is more than one way to frame a curve,” The Am. Math. Mon. 82(3), 246–251 (1975). [CrossRef]  

46. F. Klok, “Two moving coordinate frames for sweeping along a 3D trajectory,” Comput. Aided Geom. Des. 3(3), 217–229 (1986). [CrossRef]  

47. D. Bechmann and D. Gerber, “Arbitrary shaped deformations with DOGME,” Visual Comp 19(2), 175–186 (2003). [CrossRef]  

48. W. Wang, B. Jüttler, D. Zheng, and Y. Liu, “Computation of rotation minimizing frames,” ACM Trans. Graph. 27(1), 1–18 (2008). [CrossRef]  

49. J. Guan, W. Li, W. Wang, and Z. Fu, “General boundary mapping method and its application in designing an arbitrarily shaped perfect electric conductor reshaper,” Opt. Express 19(20), 19740–19751 (2011). [CrossRef]  

50. M. Fujii, C. Koos, C. Poulton, I. Sakagami, J. Leuthold, and W. Freude, “A simple and rigorous verification technique for nonlinear FDTD algorithms by optical parametric four-wave mixing,” Microw. Opt. Technol. Lett. 48(1), 88–91 (2005). [CrossRef]  

51. C. Koos, M. Fujii, C. G. Poulton, R. Steingrueber, J. Leuthold, and W. Freude, “FDTD-modelling of dispersive nonlinear ring resonators: Accuracy studies and experiments,” IEEE J. Quantum Electron. 42(12), 1215–1223 (2006). [CrossRef]  

52. F. Negredo, M. Blaicher, A. Nesic, P. Kraft, J. Ott, W. Dörfler, C. Koos, and C. Rockstuhl, “Fast and reliable method to estimate losses of single-mode waveguides with an arbitrary 2D trajectory,” J. Opt. Soc. Am. A 35(6), 1063–1073 (2018). [CrossRef]  

53. M. Deubel, G. von Freymann, M. Wegener, S. Pereira, K. Busch, and C. M. Soukoulis, “Direct laser writing of three-dimensional photonic-crystal templates for telecommunications,” Nat. Mater. 3(7), 444–447 (2004). [CrossRef]  

54. R. Courant, K. Friedrichs, and H. Lewy, “Über die partiellen Differenzengleichungen der mathematischen Physik,” Math. Ann. 100(1), 32–74 (1928). [CrossRef]  

55. J. Schneider and S. Hudson, “A finite-difference time-domain method applied to anisotropic material,” IEEE Trans. Antennas Propag. 41(7), 994–999 (1993). [CrossRef]  

56. S. G. García, T. M. Hung-Bao, R. G. Martín, and B. G. Olmedo, “On the application of finite methods in time domain to anisotropic dielectric waveguides,” IEEE Trans. Microwave Theory Tech. 44(12), 2195–2206 (1996). [CrossRef]  

57. M. Nauta, M. Okoniewski, and M. Potter, “FDTD method on a Lebedev grid for anisotropic materials,” IEEE Trans. Antennas Propag. 61(6), 3161–3171 (2013). [CrossRef]  

58. M. Salmasi, M. E. Potter, and M. Okoniewski, “Implementation of general dispersive anisotropic materials in Lebedev FDTD,” IEEE Trans. Antennas Propag. 66(12), 7171–7179 (2018). [CrossRef]  

59. C. Fong, “Elliptification of rectangular imagery,” arXiv:1709.07875 (2019).