Supplementary information for
Foundry-Fabricated Grating Coupler Demultiplexer
Inverse-Designed via Fast Integral Methods
Constantine Sideris,
1
∗
Aroutin Khachaturian,
2
Alexander D. White
3
Oscar P. Bruno,
4
†
Ali Hajimiri,
2
†
1
University of Southern California
Ming Hsieh Department of Electrical and Computer Engineering,
Los Angeles, CA 90089, USA
2
California Institute of Technology
Electrical Engineering, Pasadena, CA 91125, USA
3
Stanford University
Department of Electrical Engineering and Ginzton Laboratory
Stanford, California 94305, USA
4
California Institute of Technology
Computational and Mathematical Sciences
Pasadena, CA 91125, USA
†
These authors contributed equally.
∗
To whom correspondence should be addressed; E-mail: csideris@usc.edu.
1
Supplementary Notes
The following sections provide additional details concerning the computational design meth-
ods used as well as further analysis of the character of the proposed GCWD device. Thus,
Section Supplementary Note 1 presents details concerning the Maxwell solver and geometry
representation methods used (Sections Supplementary Note 1.1 and Supplementary Note 1.2
respectively) as well as comparisons with a state-of-the-art commercial solver (Section Sup-
plementary Note 1.3). Then, Section Supplementary Note 2 presents the actual geometric pa-
rameters that characterize the fabricated design, and Section Supplementary Note 3 presents a
simulation-based study of the robustness of the GCWD device with respect to variations in the
fabrication process and perturbations in the orientation of the incident illumination. For refer-
ence, Section Supplementary Note 4 presents and analyzes an alternate design, for which the
minimum feature-size constraints were significantly relaxed, leading to a GCWD structure with
50nm minimum-size features. Section Supplementary Note 5, finally, considers the potential of
the proposed GCWD device as a wavelength multiplexer or combiner, and presents simulation
results in support of this use case scenario.
Supplementary Note 1
Green function-based simulation and optimization for photonic
waveguide structures
This section describes several recent innovations, introduced in refs. [S1, S2], which have
facilitated efficient modeling and optimization of photonic grating structures by means of Green
function-based methods. Section Supplementary Note 1.1 thus presents, in the context of the
relevant literature, one of the main enabling elements, namely, the windowed Green function
method (WGF) for accurate termination of infinite integration domains, along with some de-
tails concerning the underlying quadrature problem. The geometry representation strategy used,
which is well suited for the photonic optimization problem at hand, is then described in Sec-
tion Supplementary Note 1.2, together with the device-optimization algorithm, including, in
particular, a description of a matrix update approach utilized for efficient implementation of
the adjoint method. Additionally, details on the application of the L-BFGS quasi-Newton algo-
rithm [S3] are provided. The efficiency and accuracy of the Green-function simulation approach
is demonstrated in Section Supplementary Note 1.3, which includes, in particular, detailed com-
parisons vis-
`
a-vis alternative methods.
Supplementary Note 1.1
Windowed Green function method (WGF)
Most available numerical methods for waveguide problems rely on use of either finite-
difference or finite-element volumetric discretizations. These solvers truncate the computational
domain, including the unbounded waveguides, on the basis of absorbing boundary conditions
such as the Perfectly Matched Layer technique (PML). As is well known, at the very least a few
points per wavelength must be used in the discretization since, clearly, it would be impossible
to resolve the details in the electromagnetic field with discretizations that contain one or fewer
points per wavelength. While necessary, the use of a fixed number of points per wavelength
in this context is generally not sufficient to maintain a prescribed accuracy in general solutions
of the Maxwell equations. To appreciate this fact, we consider an experiment in which a fixed
number of points per wavelength is used for the solution of problems for smaller and smaller
wavelengths, for a fixed overall structure. At a fixed number of points per wavelength, the
error in approximation of continuous derivatives by discrete derivatives remains constant, but
the number of times such approximations are used grows as the wavelength decreases and the
number of wavelengths spanned by the domain grows. This leads to so-called “dispersion er-
rors” that grow without bound as the number of wavelengths spanned by the domain grows.
Equivalently, to maintain a certain accuracy, the number of points used per wavelength must be
increased as the electrical size of the problem grows. In all, in view of the fine meshes needed
and the PML regions necessary to effectively eliminate reflections, volumetric methods give
rise to high computing costs. Still, finite-difference and finite-element methods can be success-
fully implemented [S4–S6], and they remain powerful, widely used numerical techniques for
the types of problems under consideration.
Boundary Integral methods can significantly reduce the required numbers of unknowns on
two main counts: as a result of the reduced dimensionality, on one hand—since they rely on use
of surface, rather than volume, mesh discretizations—and in view of their virtually dispersion-
less character—as the Green’s function analytically propagates oscillatory fields, without errors,
for arbitrarily long distances. The full system of integral equations used for electromagnetic so-
lution in the context of this paper are detailed in [S1]; to provide an illustrative description in
subsequent sections, however, we also consider here the simpler integral equation
1
2
ψ
(
r
) +
∫
Γ
∂G
(
r
,
r
′
)
∂n
(
r
)
ψ
(
r
′
)
ds
r
′
=
−
∂u
inc
(
r
)
∂n
(
r
)
(SE1)
over a single
unbounded
curve
Γ
which, outside of some bounded region, coincides with the
x
axis in the two-dimensional
r
= (
x,y
)
plane. Here
G
(
r
,
r
′
) =
i
4
H
(1)
0
(
k
|
r
−
r
′
|
)
denotes the
two-dimensional Green function for the spatial frequency
k
. This integral equation corresponds
to the problem of electromagnetic scattering of a TM polarized incident electromagnetic wave
u
inc
by a perfectly conducting surface that is invariant along the
z
axis, and whose cross-section
with the plane
z
= 0
is, precisely, the curve
Γ
. While different from the integral equation
problems presented in [S1] for a dielectric structure of the type depicted in Figure 1 of the
main text, which amount to systems of several equations similar in form to (SE1), this single
equation provides a valuable framework to demonstrate some of the main principles underlying
the numerical method used for the full integral equation system.
Prior to the recent work [S1, S2], BIE methods had almost exclusively been relegated to
applications concerning bounded scattering structures such as radiating antennas, aircraft and
spacecraft, etc., due to the lack of a suitable condition for termination of the infinite com-
putational domain arising from a semi-infinite waveguide. The challenge is that a truncation
of the infinite integration domain in SE1 creates artificial edges which cause the Green func-
tion to produce strong artificial edge reflections, thereby fatally compromising the accuracy of
the solution. The windowed Green function method (WGF) [S1] eliminates this difficulty by
multiplying the integrand in the integral of the form (SE1) over the unbounded curve
Γ
by a
smooth and
slowly-rising
windowing function
w
=
w
(
x
)
(similar to the Hann window utilized
for other purposes in signal processing) which smoothly truncates the integral to the bounded
curve
Γ
∩{
w
6
= 0
}
, and, thus reduces the numerical errors to levels that decay extremely fast
as the window-size
A
—faster, in fact, than
O
(1
/A
p
)
for all positive integers
p
[S1, S2, S7]. (An
additional consideration is necessary concerning the incoming excitation, which is omitted in
the present description for brevity, but which is detailed in the aforementioned contributions for
various types of incident field configurations. In all cases these considerations provide certain
correct prescriptions for the right-hand side of the equation for which the aforementioned fast
convergence to the true Maxwell solutions takes place as
A
→∞
.)
Two alternatives to the WGF method have been proposed. The previous physically-motivated
Surface Conductive Absorber (SCA) approach [S8], on one hand, utilizes a matching layer char-
acterized by a surface conductivity, whose absorption, as in the volumetric version of the PML
approach, results in field decay. Unfortunately, although the SCA approach can absorb waveg-
uide modes, this method generally requires large propagation lengths for adequate absorption
of scattered fields—since it does not absorb radiating waves that are generated from scattering
events other than conduction along waveguides. As a result, the SCA methodology has not
found application in practice. The PML contribution [S9] for integral methods, in turn, refor-
mulates the WGF methodology, incorporating a windowing effect by means of a complexified
distance in the Green-function exponential, and thus it provides a variant of the WGF method.
In addition to the termination problem, the novel integral solver [S1] incorporates a num-
ber of specialized techniques, including 1) An effective methodology, described in Section 3.5
of [S10], for integration of the Green-function singularity with high accuracy on the basis of
small number of discretization points (and, thus, small numbers of unknowns); 2) A change of
variables around corners, described on pp. 83-84 in [S10], whose Jacobian eliminates infinite
currents near edges; and 3) A methodology that addresses a structural characteristic, namely
nearly touching boundaries that result for small inter-tooth widths and/or spacings, as well as
for large etch depths for which the bottom of one or more spacings closely approaches the
buried oxide layer, that is specific to the simulation of the grating couplers of the type consid-
ered in this paper: instead of utilizing large numbers of unknowns to guarantee a sufficiently
fine discretization of the Green function singularity, the method continues to utilize coarse dis-
cretizations, as in other regions, but it then interpolates the unknown integral-equation solution,
on the basis of certain Chebyshev expansions that underly the integration algorithm, to a suffi-
ciently fine mesh—as needed for accurate integration. As illustrated in Section Supplementary
Note 1.3, the innovations outlined above collectively enable the WGF methodology to tackle
the waveguide problems at hand with unprecedented accuracy and efficiency.
Supplementary Note 1.2
Geometry representation and efficient adjoint method implementation
As is the case in connection with the FDTD simulation method, the calculation of the gra-
dient of the objective function with respect to parameters in the integral-equation context can
be greatly accelerated by means of the adjoint method. Although the general principles for the
adjoint method can be applied in connection with Green functions and integral equations, for
efficiency, the methodology requires a number of innovations, as described in what follows (see
also Sec. SV in [S1]).
As is in the FDTD case, the adjoint method for integral equation problems proceeds by sub-
tracting from the expression for the gradient of the objective function an adequate linear com-
bination of the equations for the gradient components—which can themselves be obtained by
directly differentiating the system matrix. Upon such subtraction the resulting gradient expres-
sion does not contain derivatives of the electromagnetic solution, and can be computed without
requiring a system solution per gradient component. The coefficients of the aforementioned lin-
ear combination are obtained as the solution of a linear system whose matrix is the adjoint of the
system matrix for the electromagnetic problem, and whose right-hand sides involve derivatives
of the Green-function matrix itself.
Due to the complicated integral expressions involved, these right-hand side quantities must
be computed numerically using finite differences with respect to each optimization parameter,
which requires calculating changes in the adjoint system matrix. The geometrical parametriza-
tion approach that we use provides a significant advantage, as it gives rise to highly localized
variations in this matrix as each parameter is varied—so that computing the variation for all
possible parameters requires a computational cost comparable to the computation of a single
full system matrix instead of
2
n
+ 1
matrices, one per each of the
2
n
+ 1
parameters used.
In detail, for the particular type of grating geometry considered in this paper, the interfaces
can be easily parametrized by a sequence of parametrizations, one for each straight segment in
the structure. The straight segments are characterized by the positions
x
j
of the vertical sides of
the individual teeth along the horizontal
x
axis (which can easily be made to correspond with
the
n
widths
w
i
and
(
n
+ 1)
spacings
s
i
depicted in Figure 1 in the main text). In this manner,
a given variation of one the
x
j
positions does not cause a global modification of the Green-
function matrix. As mentioned above, this approach, from which actual widths and spacings
can then be easily recovered to enable enforcement of the box constraints, is advantageous
since such variations only give rise to small changes in the system matrix, which allows for an
efficient matrix update.
On the basis of this geometric setup, the BFGS gradient-descent algorithm [S3] with simple
box constraints
(termed L-BFGS-B in that reference) was used to optimize the design under
the given parametrization. All grating widths and spacings were restricted on the basis of the
L-BFGS-B box constraints to a minimum width of 160 nm. Typically, the L-BFGS-B optimizer
converged on the prescribed problems in 50 to 200 iterations.
Supplementary Note 1.3
Comparison vis-a-vis FDTD solvers
The WGF formulation described above provides a very significant advantage in the context
of device optimization. Indeed, since only boundaries are discretized in the WGF approach, the
discretization points are carried along with the boundary as part of the optimization process,
and, thus, the method avoids difficulties inherent in volumetric methods, as discussed in the
introductory section in [S1]. Additional details concerning the geometry character of the WGF
method and associated device-optimization approach are presented in Section Supplementary
Note 1.2.
The WGF-BIE also provides significantly improved solution capabilities over leading meth-
ods in practical engineering design and optimization contexts. To demonstrate the advantages
derived from WGF-BIE in this regard we consider the accuracy levels that are achievable by
a leading commercial FDTD-based solver. A comparison of the accuracy of the field solu-
tion of the final optimized GCWD structure in the region where the mode overlap integrals are
computed for optimization of the device is provided in Supplementary Figure 1, which uses a
very finely discretized solution for each solver as a reference. We see that, in order to achieve
an accuracy of approximately
1%
(which is necessary to avoid degradation of the typical
10%
experimental accuracy
1
) requires computing times of 10 sec and 2000 sec for the BIE and
1
Clearly, a notional
10%
solution error could lead to device sub-optimality by up to
20%
when combined with
(a)
(b)
Supplementary Figure 1. Proposed and commercial solvers: Errors and computing times.
Relative errors resulting from the WGF-BIE and commercial FDTD solvers, evaluated via com-
parison with finely resolved solutions, versus (a) the number of points per wavelength, and
(b) required computing time. Single core runs. The same computer (2.9GHz Intel Core i9) was
used for both the WGF-BIE and commercial-solver runs. For a
10
−
2
relative error level the
required BIE and FDTD solution times are approximately 10 sec and 2000 sec, respectively.
FDTD solvers respectively. In particular, clear advantages are provided by the proposed WGF-
BIE solver the context of device optimization, wherein hundreds of solutions must obtained to
achieve an optimized design.
Supplementary Note 2
Grating Coupler Wavelength Demultiplexer design dimensions
Supplementary Table 1 contains the complete set of widths and spacings
w
j
and
s
j
(see
Figure 1 in the text) that were used in the fabricated GCWD design presented in this work. The
thicknesses of the various layers are as follows:
• Passivation layer:
2
.
78
μ
m
• Si device (upper Si layer):
0
.
22
μ
m
• Buried oxide layer:
2
.
0
μ
m
an underlying
10%
fabrication error.
j
w
j
s
j
0
—
1.00
1
0.46
0.16
2
0.16
0.41
3
0.16
0.16
4
0.45
0.16
5
0.16
0.42
6
0.16
0.65
7
0.16
0.23
8
0.45
0.26
9
0.18
0.36
10
0.16
0.55
11
0.33
0.16
12
0.39
0.37
13
0.16
0.32
14
0.28
0.49
15
0.16
0.27
16
0.37
0.35
17
0.16
0.62
18
0.16
0.16
19
0.16
2.87
Supplementary Table 1. Proposed GCWD design parameters.
Widths and spacings
w
j
and
s
j
that characterize the fabricated GCWD design (in
μ
m).
• Si substrate (lower Si layer):
725
μ
m (computationally modeled as infinitely thick)
Refractive indices of 3.48, 1.44, and 1.0 for silicon, oxide, and air, respectively, were used
in the numerical simulations of the device. A TM-polarized vertical Gaussian beam with a
10
μ
m waist diameter at a height of 0.22
μ
m above the GCWD passivation layer was used as an
approximation of the incident excitation provided by the fundamental mode of the illuminating
SMF-28 optical fiber.
Supplementary Note 3
Grating Splitter Robustness Evaluation
This section discusses the variation in performance experienced by the proposed GCWD
device in response to fabrication defects as well as perturbations in the orientation of the incident
illumination.
Supplementary Note 3.1
Impact of fabrication process variations
We consider results of a Monte Carlo-based statistical analysis on the performance vari-
ability of the final optimized GCWD with regards to geometrical perturbations, demonstrating
the device’s robustness with regard to fabrication defects. To do this we evaluated the perfor-
mance of
855
,
893
different random variations of the proposed GCWD obtained by adding to
each of the vertical wall position
x
j
(as described above in Section Supplementary Note 1.2)
an independent zero-mean normally distributed random number with 10nm standard deviation.
(A standard deviation of
10
nm was chosen as a conservative over-estimate of typical foundry
process variation.) Supplementary Figure 2 displays histograms of the resulting coupling effi-
ciencies and splitting ratios. As can be seen, both parameters are fairly insensitive to random
variations in the design parameters. For example, for 2
σ
variation, which includes
95%
of all
the samples, the coupling efficiencies for the right and left outputs at 1.55
μ
m and 1.3
μ
m respec-
tively vary by at most
±
3
.
3%
and
±
2
.
7%
; the corresponding splitting ratios exhibit a maximum
deviation of
±
3
.
8
dB for both wavelengths. Thus, the proposed GCWD design is expected to be
robust to process variation—a fact that is corroborated by the agreement observed between the
simulated and measurement performance results.
(a)
(b)
Supplementary Figure 2. Monte Carlo simulations of GCWD process variation.
His-
togram depiction of the variation observed in the coupling efficiency over a total of
855
,
893
samples obtained as random geometric perturbations of the proposed GCWD. The random de-
vice samples were generated by adding normally distributed zero-mean random perturbations,
with 10nm standard deviation, to the GCWD wall positions determined by the widths and spac-
ings presented in Supplementary Table 1. Panels (a) and (b) present the resulting histogram
plots of coupling efficiencies and splitting ratios, respectively.
Supplementary Note 3.2
Sensitivity to incident fiber angle
Although the GCWD was optimized for and intended to be excited via a perfectly vertical
(
θ
= 0
) incident fiber, in a practical experimental setup or packaged module a small tilt may
occur in the fiber angle which could affect the performance metrics. (Here
θ
is defined as the
angle between the fiber and the direction normal to the propagation axis.) In order to study
the sensitivity of the GCWD structure’s insertion loss and splitting ratio parameters (defined
in the main text), we performed a set of simulations where we swept the incident fiber angle
θ
from -10 to 10 degrees in the plane that contains the fiber and the GCWD propagation axis.
Supplementary Figure 3 displays the insertion loss and splitting ratio versus coupling angle for
both wavelengths,
λ
1
.
3
μ
m and
λ
1
.
5
μ
m.
(a)
(b)
Supplementary Figure 3. GCWD performance sensitivity to incident fiber angle.
Simu-
lated GCWD insertion loss (a) and splitting ratio (b) as functions of the fiber angle
θ
relative to
the normal.
According to Figure 3, even up to
θ
=
±
5
◦
, a reduction of insertion losses by no more than
approximately 2dB is observed. The splitting ratios are even less sensitive to negative angle
perturbations (only 1dB worse at
−
5
◦
for the 1.55
μ
m port and, in fact, better in some cases
for the 1.33
μ
m port); however, the splitting ratio for the
1
.
3
μ
m wavelength at a
+5
◦
angle is
degraded by up to 4.5dB. This suggests that the fiber should be deliberately tilted by a small
angle in the negative
θ
direction in an experimental setup in order to ensure that any error falls
on the side where the GCWD is the least sensitive to deviations from
0
◦
incidence.
Supplementary Note 4
Relaxing minimum feature-size constraints
Unlike UV lithography, E-beam lithography is time consuming, expensive, and does not
scale for mass production. However, E-beam lithography allows etching of much smaller fea-
tures, which could potentially lead to improvements in the device performance. Since the in-
verse design methodology we use is capable of designing grating splitters with any minimum
and maximum feature size constraints, as dictated by the fabrication process, a study was con-
ducted seeking to estimate the potential performance gains that might result from relaxing the
feature-size constraint to 50nm (from 160nm which was implemented in our proposed GCWD
design, in accordance with our foundry design rules). We thus used our algorithm to design a
new grating splitter with the same objective function as the one presented in the paper. We in-
creased the number of widths to 79 and number of spacings to 80 for a total of 159 optimization
parameters so that we could keep the total size dimension of the grating the same as that of the
proposed device while allowing for smaller features to span the length of the design.
j
w
j
s
j
0
—
0.050
1
0.063
0.087
2
0.050
0.050
3
0.109
0.166
4
0.052
0.050
5
0.054
0.125
6
0.050
0.050
7
0.050
0.147
8
0.129
0.050
9
0.050
0.099
10
0.050
0.050
11
0.077
0.079
12
0.050
0.050
13
0.127
0.148
14
0.147
0.162
15
0.058
0.050
16
0.050
0.092
17
0.050
0.050
18
0.084
0.156
19
0.050
0.050
j
w
j
s
j
20
0.096
0.115
21
0.050
0.050
22
0.050
0.169
23
0.165
0.051
24
0.050
0.050
25
0.078
0.110
26
0.050
0.064
27
0.162
0.050
28
0.050
0.122
29
0.062
0.050
30
0.050
0.160
31
0.050
0.050
32
0.068
0.050
33
0.054
0.050
34
0.050
0.206
35
0.107
0.050
36
0.050
0.111
37
0.071
0.199
38
0.050
0.050
39
0.097
0.105
j
w
j
s
j
40
0.093
0.108
41
0.050
0.050
42
0.150
0.050
43
0.077
0.162
44
0.050
0.050
45
0.050
0.050
46
0.050
0.095
47
0.122
0.050
48
0.065
0.050
49
0.050
0.167
50
0.079
0.107
51
0.079
0.050
52
0.050
0.132
53
0.053
0.066
54
0.050
0.050
55
0.129
0.133
56
0.073
0.050
57
0.050
0.057
58
0.202
0.095
59
0.065
0.127
j
w
j
s
j
60
0.050
0.050
61
0.050
0.148
62
0.137
0.050
63
0.050
0.050
64
0.081
0.112
65
0.050
0.050
66
0.050
0.085
67
0.157
0.116
68
0.082
0.050
69
0.050
0.083
70
0.147
0.137
71
0.050
0.050
72
0.050
0.050
73
0.050
0.050
74
0.050
0.084
75
0.163
0.134
76
0.152
0.105
77
0.050
0.050
78
0.050
0.050
79
0.050
1.860
Supplementary Table 2. Reduced minimum feature-size GCWD design (for reference).
Widths and spacings
w
j
and
s
j
that characterize the 50nm minimum feature-size GCWD design
(in
μ
m).
As expected, the resulting simulated performance improved appreciably. For the new de-
sign, the coupling efficiencies into the right and left output waveguides at 1.55
μ
m and 1.3
μ
m
equal 52.1% and 66.7% respectively, compared to 37.6% and 27.5% for the grating obtained un-
der the 160nm minimum feature size constraints. The splitting ratios of the new 50nm relaxed-
constraint coupler, in turn, equal 23.5dB and 21.6dB at 1.55
μ
m and 1.3
μ
m respectively, which
also represents an improvement, albeit a lesser one, over the 160nm constrained design—for
which splitting ratios of 19dB and 21dB at 1.55
μ
m and 1.3
μ
m, respectively, were obtained. As
mentioned previously, although the relaxed design constraints result in improved simulated per-
formance, such a coupler is not foundry compatible and may also be more susceptible to fabri-
cation error than a grating with minimum feature size of 160nm. For reference Supplementary
Table 2 lists the optimized width and spacing parameters of the 50nm minimum feature size
grating obtained. All stack-up and excitation parameters were kept unchanged.