## Abstract

Extreme skin depth engineering (*e-skid*)
can be applied to integrated photonics to manipulate the evanescent
field of a waveguide. Here we demonstrate that *e-skid* can be implemented in two directions in order to
deterministically engineer the evanescent wave allowing for dense
integration with enhanced functionalities. In particular, by
increasing the skin depth, we enable the creation of two-dimensional
(2D) *e-skid* directional couplers with
large gaps and operational bandwidth. Here we experimentally validate
2D *e-skid* for integrated photonics in a
complementary metal–oxide semiconductor (CMOS) photonics foundry and
demonstrate strong coupling with a gap of 1.44 µm.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. EVANESCENT WAVES IN SILICON PHOTONICS

Evanescent waves in silicon photonic waveguides have the propensity to cause parasitic optical cross talk. In traditional photonic circuits, design strategies must consider minimum separation distances between any closely spaced waveguides to prevent unwanted coupling [1]. This problem is inhibiting for many photonic circuits due to cost and size constraints. Many efforts have been made to overcome these issues, battling size constraints by employing inverse design [2,3] or implementing metamaterials to increase performance [4–6].

Recent work introduced a new, metamaterial paradigm for waveguiding that
fundamentally suppresses coupling between waveguides [7]. In this approach, a subwavelength,
multi-layer cladding is placed in plane and in parallel with the
waveguide, decreasing the skin depth of the fundamental
transverse-electric (TE) mode’s evanescent field. The concept is called
*extreme skin depth engineering*, or *e-skid*. *E-skid* has
been employed as cross talk suppression [7,8], and for high
performance polarization splitting [9,10]. The *e-skid* features are created in the same processing
step as the waveguide itself, allowing this to be an innate no-cost
addition to any design. The addition of these features can reduce the
cross talk between waveguides by more than 3 orders of magnitude, which
will dramatically reduce the photonic design footprint [7].

Here we expand on this work by using *e-skid*
engineering in *two directions*, within the
same plane as the waveguide. Using both a parallel (as in [7]) and perpendicular *e-skid* cladding, we engineer the coupling between waveguides
throughout a photonic circuit. Specifically, the perpendicular *e-skid* cladding (with features orthogonal to the
waveguide) can dramatically increase the evanescent tail of the mode [or
equivalently decrease its decay constant—as shown in Fig. 1(c)]. This then enables arbitrary
enhancement in coupling between nearby waveguides anywhere in the circuit.
While parallel *e-skid* cladding can be used
to suppress coupling in the routing of the rest of the circuit. Here we
explore this with the design of a two-dimensional (2D) *e-skid* directional coupler with a large gap ($\ge\!
1.4\,\,\unicode{x00B5}{\rm m}$) between the waveguides and large
operational bandwidth ($\ge\! 40\;{\rm nm}
$). We also demonstrate these 2D *e-skid* directional couplers in a complementary
metal–oxide semiconductor (CMOS) photonic platform, thereby affirming the
manufacturability of *e-skid* components and
integration with foundry offerings. Ultimately, this work is motivated by
the idea that we can employ 2D *e-skid*
techniques to ensure low cross talk in the routing phase (with parallel
*e-skid*) and, without the use of bends,
efficient coupling between waveguides with larger-than-normal gaps (with
perpendicular *e-skid*)—fully leveraging both
directions of *e-skid*. Consequently, 2D
*e-skid* allows for dense integration within
the constraints of CMOS manufacturing.

## 2. THEORY

#### A. *E-Skid* in Two Directions

Consider two media with index of refraction ${n_1},{n_2}$. When an incoming wave from ${n_1}$ meets the boundary at ${n_2}$ and the angle is greater than the
critical angle, ${\theta _i} \lt {\theta
_c} = {{\sin}^{- 1}}({n_2}/{n_1})$, an evanescent wave is formed in the
second medium. This wave does not carry power across the boundary; it
exponentially decays into the second medium [11]. *E-skid* allows us
to tune the decay constant of this evanescent wave by introducing
subwavelength, periodic structures that transform a wave’s
(specifically, a polarized wave’s) momentum [12,13]. These
features change the second medium from an isotropic material to an
anisotropic metamaterial. The anisotropy here refers to the
permittivity values of the dielectric tensor of the material (where we
are assuming that the permittivity can be defined ${\epsilon _r} =
{n^2}$). For deep-subwavelength features,
these component values are defined by the Rytov relations [14,15]:

The key result of the *e-skid* derivation
leverages this anisotropy for the evanescent wave, which is
characterized by the decay constant, $\beta$:

In Fig. 1(a), we show the
dielectric tensor for the *e-skid*
structure, where the periodicity of the subwavelength features is
parallel to the boundary ($y = 0$). In this orientation, the diagonal
components of the second material become $[n_{2x}^2,n_{2y}^2,n_{2z}^2] = [n_\parallel ^2,n_ \bot
^2,n_\parallel ^2]$ in accordance with the Rytov
relations [Eqs. (1) and
(2)]. This structure
will increase the decay constant of the evanescent wave, thereby
decreasing the skin depth [7].
Without loss of generality, we recognize that we can rotate the
optical axis by rotating the subwavelength features and realize
*e-skid* in a *second* direction. Due to the direction dependency outlined
by the Rytov relations, when we rotate the periodicity of the
features, we effectively swap the $xx$ and $yy$ components of the dielectric tensor
of the parallel cladding such that we now see $[n_{2x}^2,n_{2y}^2,n_{2z}^2] = [n_ \bot ^2,n_\parallel
^2,n_\parallel ^2]$ [Fig. 1(b)]. The values of ${n_{2x}},{n_{2y}}$ in Eq. (3) control the decay constant, and by rotating the
periodic structure, we are able to dictate an increase or
decrease.

We populated Eq. (3) with
the new dielectric tensor values outlined in Figs. 1(a) and 1(b) such that we show in Fig. 1(c) the full range of decay constant tunability of
*e-skid* in two directions. Figure 1(c) shows clearly that both
decreasing and increasing skin depth can be achieved by the parallel
features; however, applying this to CMOS photonics manufacturing, we
generally omit the higher and lower fill factors due to resolution
constraints [16]. We used a
material platform consistent with CMOS photonics in Fig. 1(c), such that material one is
silicon (Si) and material two is silicon dioxide (${{\rm
SiO}_2}$); however, this is true for any
optical material combination as long as ${n_1} \gt
{n_2}$.

#### B. *E-Skid* in Two Directions in Waveguides

Optical waveguiding is not fully described by the simple
electromagnetic wave-at-a-boundary example above. While it lends
intuition, we must find the electromagnetic mode of the entire
structure to get a clear picture of this effect. We used a commercial
finite difference eigenmode (FDE) solver to simulate the TE modes of
three specific types of waveguides to demonstrate *e-skid* in two directions [17]. Figure 2(a) shows
a top view of a single-mode strip waveguide, where the propagation is
in the $\hat x$ direction. Next to the strip
waveguide, Fig. 2(d) shows a 2D
mode profile of the fundamental TE propagating mode. We introduce the
wave supressing *e-skid* features on one
side of the waveguide in Fig. 2(b) and show the corresponding 2D mode profile in Fig. 2(e). Finally, we introduce the wave
enhancing *e-skid* features in Fig. 2(c) and the corresponding 2D mode
profile in Fig. 2(f). We
compiled the cross sections of all three modes in Fig. 2(g) to demonstrate the effect of the
features on the evanescent wave of the mode. Figure 2(g) clearly demonstrates, with a
log-scale in $y$, that the decaying wave outside of
the center of the waveguide is suppressed by the parallel features,
and greatly enhanced by the perpendicular features.

We extend our analysis to include a full wave three-dimensional (3D) finite-difference time-domain (FDTD) solver with Bloch-periodic boundary conditions [17]. Figure 3 displays the decay constant as extracted for varying fill factors and periods of subwavelength features. The decay constant extraction fit for both Figs. 3(a) and 3(b) had average correlation coefficients of ${R^2} = 0.98$. In Fig. 3(a), we show the results of the parallel cladding, and in Fig. 3(b), we show the perpendicular, both of which follow the general shape outlined by the analytical expression in Eq. (3) at the paraxial limit, or simply, a guided wave. As the period increases, we notice a shift in behavior of the peak toward the lower fill factors; however, this is not significant enough to alter the design parameters. We also note the similarity of Figs. 3 and 1(c), where the parallel features are the only features that increase the decay constant and the perpendicular features allow for deterministic reduction of the decay constant, thereby enhancing cross talk. Figure 3 stops at 270 nm to allow each simulation to fall below the Bragg limit. We note that the operational window for subwavelength devices is dependent on the minimally reliable feature (or hole) size for the manufacturing process. Here, we see at the highest period that only 30–70% of fill factors are expected to resolve for an 80 nm minimum.

## 3. 2D *E-SKID* DIRECTIONAL COUPLER
DESIGN

We propose and demonstrate a new coupler that leverages *e-skid* in two directions to create coupling in
desired regions. In a traditional integrated photonic platform, a
directional coupler is created by bending two waveguides close to each
other [Fig. 4(a)]. The waveguides
must otherwise be kept far apart in other parts of a circuit in order to
avoid unwanted coupling, limiting the circuit density. By using *e-skid* with parallel subwavelength features
[Fig. 4(b)] that suppress coupling,
we overcome this limitation and keep two waveguides within close
proximity, with negligible coupling. Furthermore, when the design with
*e-skid* needs coupling, we show that the
introduction of perpendicular subwavelength features in the coupling
region, as are seen in Fig. 4(c),
will significantly enhance coupling. These features have tunable variables
(i) period (${\Lambda _ \bot}$) and (ii) fill factor (${\rho _ \bot}$), which directly tune the amount of
coupling experienced. We introduce 2D *e-skid*
as a way to leverage the size reduction offered by the parallel features
with the addition of the perpendicular features to create practical
circuits as seen in Fig. 4(d).

#### A. Coupled Modes for *E-Skid*

The evanescent wave of the mode, even though it carries no power across the boundary, causes coupling between parallel guides if the overlap between the evanescent waves of supported modes is large enough [18]. From coupled mode theory [18], we define the power in the bar and cross ports as

This approach allows for an intuitive understanding of the device. The crossover length is given as

where $\lambda$ is the free space wavelength and ${n_{{\rm even},\lambda}},{n_{{\rm odd},\lambda}}$ are the effective indices of the even and odd modes, respectively [1]. The field of the odd mode is antisymmetric across the coupling region, and it remains generally unaffected by symmetric features there [19]. However, by introducing features into the coupling region, the even mode is affected, thereby enabling dispersion engineering of the directional coupler—specifically, controlling the directional couplers’ optical bandwidth [19]. By crafting ${L_x}$, we can dictate how the device performs according to Eqs. (4) and (5). Essentially, if we make the slope of ${L_x}$ as flat as possible over a span of $\lambda$, we ensure a useful operating bandwidth (e.g., a 3 dB coupler) is preserved for that span. Our design is fundamentally different than [19] due to the structural asymmetry, which encourages coupling, and the higher fill factor. These parameters allow us to create a directional coupler with more than an order of magnitude shorter crossover length in comparison, at the penalty of reduced operating bandwidth.In order to demonstrate the effect of dispersion engineering, we begin
by simulating the photonic band structure of the directional coupler
in a full wave 3D FDTD solver with Bloch-periodic boundary conditions
[17]. Figures 5(a)–5(c) show the photonic band structure of a traditional,
one-dimensional (1D) *e-skid* and 2D
*e-skid* directional coupler. These
directional couplers are fundamentally different from photonic
crystals as they are not designed to work in the photonic bandgap;
instead, these subwavelength features allow for low loss propagation
through the periodic structures below the bandgap [15]. The traditional and *e-skid* couplers exhibit similar band structures,
but the 2D *e-skid* directional coupler’s
even mode is approaching the band edge just above 200 THz (1500 nm).
Because the even mode is near the band edge, dispersion is increased,
which allows for flexibility in tuning the behavior.

From the band structures, we extract the dispersive effective index of
both the fundamental even and odd supermodes [Figs. 5(d)–5(f)]. The effective indices exhibit a similar characteristic
shape to their corresponding band structures. The *e-skid* coupler brings the even mode effective index much
closer to the odd mode in comparison with the traditional coupler, and
the 2D *e-skid* directional coupler has
increased the difference between the two effective indices.
Figures 5(g)–5(i) show the crossover length given
the corresponding effective indices. The traditional coupler and
*e-skid* behave as expected, with an
increase in crossover length for the latter. The 2D *e-skid* directional coupler exhibits a
dramatically reduced crossover length and, in relation to dispersion
engineering, a completely different shape. It is important to note
that the 2D *e-skid* directional coupler
does support a higher-order mode [Fig. 5(c)]; however, the coupling efficiency extracted from a modal
overlap integral between the fundamental and the first higher-order
mode is ${\lt}8\%$ over the wavelength span for our
designs, according to our 3D FDTD simulations [17]. While 8% is not insignificant, tweaking our
parameters (specifically ${\rho _
\bot}$ and ${\Lambda _
\bot}$) can reduce this coupling efficiency
into higher-order modes, increasing device performance [19]. For this experiment, the primary
design choices were dictated from a perspective of manufacturability.
In the future, small decreases to ${\rho _
\bot}$ and ${\Lambda _
\bot}$ will result in higher performing
couplers verified by our 3D FDTD simulations [17] and prior work [19]. Additionally, careful
consideration of tapering, which is not investigated in this work, can
also mitigate the excitation of higher-order modes [15,20].

The 2D *e-skid* directional coupler is
fundamentally different to a multimode interferometer (MMI). MMIs
offer compact power splitting via a self-imaging effect that manifests
from the excitation of higher-ordered modes, creating an interference
pattern that repeats at a length specific to the physical parameters
[21,22]. The 2D *e-skid*
directional coupler specifically operates as a *directional coupler*, where the characteristics of the device
arise from the coupled supermodes supported by the two waveguides
[18]. The perpendicular,
metamaterial features act to control the decay constant of the two
modes’ evanescent tails by shaping the material between the
waveguides, enabling 2D *e-skid* and
dispersion engineering, as demonstrated by the simulation results in
Fig. 3.

#### B. Device Design

We investigate the effect of different parameter variations of the 2D
*e-skid* directional coupler. Figure 6 shows the results of varying the
fill factor, period, and gap between waveguides. These parameter
variations indicate the substantial tunability offered by
two-directional *e-skid*. First, we
selected the operating gap between the two waveguides to be 1.44 µm in
order to stay consistent with the parallel *e-skid* features (${\rho _\parallel} = 60\%
,{\Lambda _\parallel} = 225\;{\rm nm} $, six layers deep results in a 1.44 µm
gap). We designed these devices for manufacturing with the American
Institute of Manufacturing (AIM) Photonics CMOS foundary Multi-Project
Wafer (MPW) offering. For a photolithographic process like this one,
we must take in to account the limitations of the processing, like
feature size. For example, many prior work designs with features
smaller than 60 nm would not resolve with CMOS processing compared to
electron beam lithography. We chose to design our devices with ${\Lambda _ \bot} =
275\;{\rm nm} $ to remain beneath the Bragg limit but
maintain high manufacturing quality. The parallel features were
designed with ${\Lambda _\parallel} =
225\;{\rm nm} $. It should be noted that the parallel
cladding structures are less challenging for a lithographic system
because they are lines, not holes [16]. We targeted ${\rho _ \bot} =
60\%$ for the majority of our devices
because we assumed that the features would be over etched, a common
practice in silicon-on-insulator fabrication, so that the fill factor
would decrease [23].

## 4. 2D *E-SKID* DIRECTIONAL COUPLER DEVICE
MEASUREMENTS AND PARAMETER EXTRACTION

#### A. Experiment

Our experimental setup is shown in Fig. 7(a). We placed the chip on a mount in between two 3-axis
stages with bare fiber on either side for coupling in and out. For the
input, we connected the fiber to a tunable laser source (TLS) and a
polarization controller (PC) to ensure TE polarization, and we
measured the output signal with an optical powermeter (OPM). The
fibers were edge coupled to the chip, which routed the light through
strip/wire waveguides to the devices. To transition from the strip
waveguide mode to the *e-skid*’s, we
slowly introduced the parallel and perpendicular claddings as seen in
Fig. 7(c) and depicted in the
schematic in Fig. 7(a). We were
careful to design simple tapers because when periodic, asymmetric
features are introduced, there is a chance for radiative losses [15]. We measured total device loss at $\le\! 2\;{\rm
dB}$ per device, which can be imrpoved
with longer tapers or carefully designed width-varying tapers to
ensure smooth modal transition.

We collected the transmission spectra [an example is shown in
Fig. 7(b)]. The spectrum shows
a characteristic “chirp”-like behavior, which is due to the dispersive
nature of the crossover length. This is seen in Fig. 6, where, based on the perpendicular
*e-skid* parameters of the device, the
crossover length can exhibit a significant change with wavelength.
This will result in a quickly oscillating output [given by Eqs. (4) and (5)]. With that said, using the
experimental measurements in conjunction with the previously described
models, we were able to extract the parametric dependence of the
coupler designs.

#### B. Parameter Extraction Method

After extracting the transmission spectra from our fabricated devices, we used a dispersive model for extracting the crossover length from the data. Because inverse sine functions are multi-valued, we cannot obtain ${L_x}$ directly from Eq. (4) or Eq. (5). We instead employed the behavioral model for characterizing directional couplers [24]. For the crossover length, we are interested in the wavelength dependence, so we prepared the data by filtering the noise and normalizing the bar and cross measurements. We used a polynomial expansion of the coupling coefficient coupled with a non-linear least squares (NLS) optimization algorithm to find the best fit for ${L_x}$ [17,25,26]. Because many of our devices exhibits a strong “chirp-like” behavior [Fig. 7(b)], we used a third-order polynomial expression for ${L_x}$ to determine the best fit, such that

where the curve is characterized by fitting parameters ${L_{x,0}},{L_{x,1}},{L_{x,2}},{L_{x,3}}$ pertaining to the wavelength, $\lambda$. The NLS optimization aimed to minimize the difference between the measured and theoretical spectra by adjusting the fit parameters in Eq. 7.#### C. Experimental Results

In Fig. 8(a), we show the
device’s dependence on fill factor variation. Fill factors up to 60%
were successfully fabricated in the CMOS process and will be reported
below, while fabrication specific optimization needs to go into higher
fill factor devices. The fill factor variations resemble those from
the simulations [Fig. 6(a)],
where higher fill factors increased ${L_x}$. In Fig. 8(b), we investigated coupler-length variation for a
fixed fill factor (60%) and gap (1.44 µm). We fit nine different
directional couplers with the exact same parameters changing only the
coupling length. According to Eqs. (4) and (5), the length is independent of the coupling coefficient, $\kappa$_{,} and therefore these
couplers should exhibit identical ${L_x}$ measurements. There are manufacturing
variations, measurement errors, and fitting errors that reveal
themselves in Fig. 8(b). The
inset shows a 95% confidence interval for these nine couplers’ ${L_x}$ extraction and showcases expected
similar behavior for all of these couplers. This truly highlights the
utility of these devices, as they are able to couple fully in $\le\!
50\,\,{\unicode{x00B5}{\rm m}}$, even with a large gap of 1.44 µm.
Finally, Fig. 8(c) shows the ${L_x}$ extraction for varying gaps. Even
though there is a qualitative match between the simulations and
experimental results, manufacturing variations account for the
quantitative differences.

## 5. FUTURE WORK

In this work, we performed a comprehensive study of the parameter space of
our proposed directional coupler using 2D *e-skid*. In the future, it will be desirable to realize devices
with particular performance characteristics. From our experimental data
(Fig. 8), we extracted that we are
able to realize a 2D *e-skid* directional
coupler that achieves 100% coupling (100/0 splitting ratio) with a
coupling length of $L =
50\,\,{\unicode{x00B5}{\rm m}}$ as seen in Fig. 9(a), using ${\rho _ \bot} = 60\%
,\;{\Lambda _ \bot} = 275\;{\rm nm} $, and coupling gap of 1.44 µm. We can take
this to design a 50/50 directional coupler. Figure 9(b) shows the theoretical transmission spectrum of this
device. We set the coupling length $L = {\rm avg}({L_x})/2 =
23.5\,{\unicode{x00B5}{\rm m}}$, and we see broadband behavior of nearly
40 nm. Additionally, we can slightly vary parameters to tune the device to
a more desirable center wavelength. For example, by reducing the period of
the 2D *e-skid* directional coupler to ${\Lambda _ \bot} = 255\;{\rm
nm} $, ${\rho _ \bot} =
50\%$, and coupling gap of 1.44 µm, the
device’s operating band can now be centered closer to 1.55 µm and achieves
an even larger operating bandwidth of ${\gt}40\;{\rm nm}
$ [Fig. 9(c)], which is sufficient for many applications.

## 6. CONCLUSION

We introduced the deterministic, targeted control of the evanescent wave in
the TE mode of strip waveguides by employing *e-skid* features in two directions. We designed and demonstrated
2D *e-skid* directional couplers on a CMOS
photonic chip fabricated by AIM Photonics. All of the results were
compared to simulations by extracting design parmaters, using a NLS
optimization technique coupled with a behavioral model of the directional
coupler [24]. With the parameter
extraction, we show experimentally that *e-skid* waveguides, and the 2D *e-skid* directional coupler in particular, are possible to
realize in a CMOS platform. Moving forward, we can design full 2D *e-skid* circuits that achieve high densities by both
suppressing and enhancing coupling at will, while operating with large
bandwidths.

## APPENDIX A: VERIFYING MODE CONVERSION EFFICIENCY OF 1D *E-SKID* FOR CMOS PHOTONICS

We verify the conversion efficiency between strip waveguides and
*e-skid* waveguides by comparing the loss
coefficients at the resonances of a racetrack resonator [27,28]. We designed three strip waveguide racetrack resonators,
Fig. 10(a), with varying gaps,
and then the exact same racetrack resonators in which we add two
*e-skid* features into the ring,
Fig. 10(c). The intent of this
experiment is to quantify the additional loss created by these
features, and thereby measure the mode conversion efficiency between
strip and *e-skid* waveguides [7]. A simple ring resonator [shown in
Fig. 10(a)] can be
paramaterized by two coefficients, $\tau$, the self-coupling coefficient that
indicates how much light goes through the coupler, and $\alpha$, the loss coefficient that indicates
how much light is lost into the ring. We extract $\alpha$ and $\tau$ according to the method described by
[27], such that

The finesse, ${\cal F}$, is defined as the ratio between the free spectral range, $\Delta {\lambda _{{\rm FSR}}}$, and the full width at half-maximum, $\Delta {\lambda _{{\rm FWHM}}}$, of each resonance. The exctinction ratio, ${\cal E}$, is defined as the ratio between the transmission maximum, ${T_{{\rm MAX}}}$, off resonance, and the minimum, ${T_{{\rm MIN}}}$, at each resonance. We can decouple $\alpha ,\tau$ in Eq. (A5) using the method further discussed in [27]. When we determine $\alpha$, we know that ${\alpha ^2}$ indicates the percentage lost into the ring, which allows us to compare these two different resonators. This resulted in average values of ${\alpha ^2} = 97.2\%$ and ${\alpha ^2} = 95.7\%$ for the two ring types [Figs. 10(a) and 10(c), respectively]. This leads to a mode conversion efficiency of 99.6%.

## Funding

National Science Foundation (1810282); Air Force Research Laboratory (FA8650-15-2-5220, FA8750-16-2-0140).

## Acknowledgment

Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of Air Force Research Laboratory or the U.S. Government. M. v. N. and S. F. P. would like to acknowledge Navin B. Lingaraju for sparking a collaboration between RIT and Purdue.

## Disclosures

The authors declare that there are no conflicts of interest related to this article.

## REFERENCES

**1. **L. Chrostowski and M. Hochberg, *Silicon Photonics Design: From Devices
to Systems* (Cambridge
University, 2015).

**2. **A. Y. Piggott, J. Lu, K. G. Lagoudakis, J. Petykiewicz, T. M. Babinec, and J. Vučković, “Inverse design and
demonstration of a compact and broadband on-chip wavelength
demultiplexer,” Nat. Photonics **9**, 374–377
(2015). [CrossRef]

**3. **B. Shen, P. Wang, R. Polson, and R. Menon, “An integrated-nanophotonics
polarization beamsplitter with 2.4 × 2.4 µm 2
footprint,” Nat. Photonics **9**, 378–382
(2015). [CrossRef]

**4. **S. Jahani and Z. Jacob, “All-dielectric
metamaterials,” Nat. Nanotechnol. **11**, 23–36
(2016). [CrossRef]

**5. **J. M. Luque-González, A. Herrero-Bermello, A. Ortega-Moñux, Í. Molina-Fernández, A. V. Velasco, P. Cheben, J. H. Schmid, S. Wang, and R. Halir, “Tilted subwavelength gratings:
controlling anisotropy in metamaterial nanophotonic
waveguides,” Opt. Lett. **43**, 4691–4694
(2018). [CrossRef]

**6. **I. Staude and J. Schilling, “Metamaterial-inspired silicon
nanophotonics,” Nat. Photonics **11**, 274–284
(2017). [CrossRef]

**7. **S. Jahani, S. Kim, J. Atkinson, J. C. Wirth, F. Kalhor, A. Al Noman, W. D. Newman, P. Shekhar, K. Han, V. Van, R. G. DeCorby, L. Chrostowski, M. Qi, and Z. Jacob, “Controlling evanescent waves
using silicon photonic all-dielectric metamaterials for dense
integration,” Nat. Commun. **9**, 1893 (2018). [CrossRef]

**8. **M. B. Mia, S. Z. Ahmed, I. Ahmed, Y. J. Lee, M. Qi, and S. Kim, “Exceptional coupling in
photonic anisotropic metamaterials for extremely low waveguide
crosstalk,” Optica **7**, 881–887
(2020). [CrossRef]

**9. **H. Xu, D. Dai, and Y. Shi, “Anisotropic
metamaterial-assisted all-silicon polarizer with 415-nm
bandwidth,” Photon. Res. **7**, 1432–1439
(2019). [CrossRef]

**10. **K. Chen, K. Yu, and S. He, “High performance polarization
beam splitter based on cascaded directional couplers assisted by
effectively anisotropic structures,” IEEE
Photon. J. **11**,
1–9 (2019). [CrossRef]

**11. **E. Hecht, *Optics*
(Pearson Education/Addison-Wesley,
2002).

**12. **S. Jahani and Z. Jacob, “Transparent subdiffraction
optics: nanoscale light confinement without metal,”
Optica **1**,
96–100 (2014). [CrossRef]

**13. **S. Jahani and Z. Jacob, “Photonic skin-depth
engineering,” J. Opt. Soc. Am. B **32**, 1346–1353
(2015). [CrossRef]

**14. **S. Rytov, “Electromagnetic properties of
a finely stratified medium,” Sov. Phys.
JETP **2**,
466–475 (1956).

**15. **P. Cheben, R. Halir, J. H. Schmid, H. A. Atwater, and D. R. Smith, “Subwavelength integrated
photonics,” Nature **560**, 565–572
(2018). [CrossRef]

**16. **B. Smith, K. Suzuki, and J. Sheats, *Microlithography: Science and
Technology* (Taylor &
Francis, 1998).

**17. **Lumerical Inc., https://www.lumerical.com/products/.

**18. **W.-P. Huang, “Coupled-mode theory for
optical waveguides: an overview,” J. Opt. Soc.
Am. A **11**,
963–983 (1994). [CrossRef]

**19. **R. Halir, A. Maese-Novo, A. Ortega-Moñux, I. Molina-Fernández, J. Wangüemert-Pérez, P. Cheben, D.-X. Xu, J. Schmid, and S. Janz, “Colorless directional coupler
with dispersion engineered sub-wavelength structure,”
Opt. Express **20**,
13470–13477 (2012). [CrossRef]

**20. **R. Halir, P. J. Bock, P. Cheben, A. Ortega-Moñux, C. Alonso-Ramos, J. H. Schmid, J. Lapointe, D.-X. Xu, J. G. Wangüemert-Pérez, Í. Molina-Fernández, and S. Janz, “Waveguide sub-wavelength
structures: a review of principles and applications,”
Laser Photon. Rev. **9**,
25–49 (2015). [CrossRef]

**21. **W. J. Westerveld and H. P. Urbach, *Silicon Photonics*
(IOP Publishing,
2017).

**22. **R. Halir, P. Cheben, J. M. Luque-González, J. D. Sarmiento-Merenguel, J. H. Schmid, G. Wangüemert-Pérez, D.-X. Xu, S. Wang, A. Ortega-Moñux, and Í. Molina-Fernández, “Ultra-broadband nanophotonic
beamsplitter using an anisotropic sub-wavelength
metamaterial,” Laser Photon. Rev. **10**, 1039–1046
(2016). [CrossRef]

**23. **W. Bogaerts, D. Taillaert, B. Luyssaert, P. Dumon, J. Van Campenhout, P. Bienstman, D. Van Thourhout, R. Baets, V. Wiaux, and S. Beckx, “Basic structures for photonic
integrated circuits in silicon-on-insulator,”
Opt. Express **12**,
1583–1591 (2004). [CrossRef]

**24. **Y. Xing, U. Khan, A. R. Alves Júnior, and W. Bogaerts, “Behavior model for directional
coupler,” in *Proceedings Symposium IEEE
Photonics Society Benelux* (2017),
pp. 128–131.

**25. **P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. Jarrod Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0
Contributors, “SciPy 1.0: fundamental
algorithms for scientific computing in Python,”
Nat. Methods **17**,
261–272
(2020). [CrossRef]

**26. **S. van der Walt, S. C. Colbert, and G. Varoquaux, “The NumPy array: a structure
for efficient numerical computation,” Comput.
Sci. Eng. **13**,
22–30 (2011). [CrossRef]

**27. **W. McKinnon, D.-X. Xu, C. Storey, E. Post, A. Densmore, A. Delâge, P. Waldron, J. Schmid, and S. Janz, “Extracting coupling and loss
coefficients from a ring resonator,” Opt.
Express **17**,
18971–18982 (2009). [CrossRef]

**28. **K. Han, S. Kim, J. Wirth, M. Teng, Y. Xuan, B. Niu, and M. Qi, “Strip-slot direct mode
coupler,” Opt. Express **24**, 6532–6541
(2016). [CrossRef]