Methodology

Finite Difference Maxwell Solvers¶

In this work, the propagation of light in various guiding structures was determined numerically and optimized according to an objective function (see section Inverse Design with Goos). A collection of solvers were used, that rely on finite differences to discretize the partial differential Maxwell's equations, to find numerical solutions. The different methods and their application will be discussed in the following.

A commercially available FDE algorithm was used to find the eigenmodes of a translation invariant waveguide. It expresses the Maxwell's equations for the given meshed waveguide geometry as a matrix eigenvalue problem, which is solved numerically in the subsequent step https://www.lumerical.com/, undefined. The calculated eigenvalue is the propagation constant $\gamma$ from which the effective index $n_{eff}$ and attenuation $\alpha$ can be calculated, while the respective eigenvector represents the spatial distribution of complex $\vec{E}$ and $\vec{H}$ fields of the mode. The reported modal attenuation is composed of contributions by material loss (e.g. loss tangent, real conductivity), loss due to boundary conditions, and a significant additive numerical error, which can be emphasized by certain boundary conditions https://optics.ansys.com/hc/en-us/articles/360034901873-Starting-with-metal-BCs-in-FDE-simulations, undefined. The loss introduced by graphene can be modeled in multiple ways, which will be examined in section Graphene Modeling. A possible strategy to extract the modes of bend waveguides is to use a conformal coordinate transformation to map the bend geometry to an equivalent straight refractive index distribution Heiblum & Harris, 1975. In conjunction with a PML boundary condition in the radially outward direction, emulating unhindered radiation away from the simulation region, radiation loss due to tight bends can be quantified. This is however only possible for relatively small bending radii, due to the numerical inaccuracy of determining $\alpha$ quickly dominating the overall loss (see section Ring Resonator).

To describe the propagation of light in structures lacking translational invariance FDTD and FDFD solvers were used. Both use a grid similar to the one used by FDE. FDTD is able to calculate the evolution of arbitrary fields by incrementing time in small steps.

By using a broadband excitation and calculating the DTFT of the recorded fields, the device performance can be evaluated for multiple frequencies simultaneously using a single simulation Taflove et al., 2013[eq. 4.6]. FDFD on the other hand is much more similar to FDE in that it calculates the field distributions at a single frequency. Using frequency-domain-based calculations can significantly improve convergence behavior (and thus computational expense), which is pronounced in the case of resonant structures Taflove et al., 2013[figs. 20.10 and 20.11].

In order to ease the computational burden imposed by high-resolution 3D FDTD simulations it is possible under the conditions outlined below, to construct an equivalent two-dimensional geometry, which exhibits essentially the same behavior as the original device. For this method to yield valid results out-of-plane interactions have to be negligible, which is the case if no vertical coupling between different slab modes is possible. Equivalently the central assumption states that the modal fields at every point can be described with sufficient accuracy by the product of the profile evaluated along the $z$-axis $\vec{\chi}(z)$ (slab mode profile) and a complex phasor $\psi(x,y)$ dependent on the position on the $(x,y)$-plane Hammer & Ivanova, 2009[generalized from eq. 7, mind the different coordinate systems]:

In particular, this condition is not fulfilled, for arbitrarily shaped slabs of silicon on a silica BOX with an air cladding, as significant coupling between the TE and TM slab modes can be introduced due to the asymmetry of the refractive index distribution in $z$-direction, which leads to the cross-coupling coefficient derived from perturbation theory (see equation \ref{eq:perturbation_theory}) not canceling, as opposed to symmetrically clad waveguides. Mode hybridization can result from the non-zero coupling between TE and TM modes, which clearly leads to $\vec{\chi}(z)$ not being invariant with respect to $x$ and $y$. Nonetheless, a symmetric configuration can be found which exhibits the desired effective index profile, thus justifying, that the simulated behavior is the behavior of a physical device. The EIM is unsuitable however to quantitatively predict the performance of devices, that will be fabricated. In this case, full 3D simulations have to be used. The VAREIM was used to reduce the problem to two dimensions, to enable a broad parameter space exploration using the inverse design method described in section Inverse Design with Goos. This was done despite the fact, that the anisotropic behavior can not be faithfully reproduced during this reduction. However, no coupling between the TE and TM slab modes is possible. Therefore it was assumed that the orientation of the E-field in the simulated devices does not change with respect to the simulations performed to extract the effective indices. Working with only TE modes, assuming the permittivity (and conductivity) of graphene to be isotropic does not introduce significant error Sorianello et al., 2015. The statements that are derived concerning the modulation ability of EA resonant devices should persist in case the assumption stated above does not hold. Even though less common, the VAREIM can also be applied to FDTD simulations FDTD simulations.

Graphene Modeling¶

\begin{figure}[h] \centering \includegraphics[width=.6\textwidth]{sr_mobility_0_4eV.pdf} \caption{Relation between scattering rate $\gamma$ and mobility $\mu$ evaluated at a chemical potential \mu_c=\qty{0.4}{eV} according to Mohsin et al., 2015.} \label{fig:scatterin_rate_mobility} \end{figure} An integral part of describing the behavior of a graphene modulator is to model the interaction between graphene and light. Three main approaches exist in literature to model the complex conductivity of graphene at optical frequencies, these will be evaluated and compared.

Approaches to calculating the graphene absorption for a given mode profile will be presented under the assumption, that the graphene sheet does not perturb the mode significantly. It will be shown that these two approaches are equivalent given the relation between permittivity and conductivity introduced before.

Complex Sheet Conductivity¶

The interaction of graphene with light is usually described in terms of a sheet conductivity depending on the photon energy $E_ph=h\nu=\hbar\omega$ interacting with the parallel component of the electric field of the light wave. Typically the complex conductivity at optical frequencies is modeled by the RPA Falkovsky, 2008Koppens et al., 2011Emani et al., 2012Kwon, 2014Majumdar et al., 2014Hu & Wang, 2018 in the local limit, or using the Kubo formalism Hanson, 2008Vakil, 2012Mohsin et al., 2015Yang et al., 2014Koester & Li, 2014. Alternatively an experimentally verified closed form fit by Chang et al. (2014) is used Sorianello et al., 2015Romagnoli et al., 2018. \begin{figure}[hp!] \centering \includegraphics[width=\textwidth]{kubo_srs_lam1.55_Ts.pdf}

\caption{Doping dependent complex conductivity of graphene for $\lambda_0=\qty{1.55}{um}$ according to the Kubo model for different temperatures $T$ and scattering rates $\gamma$. The separate contributions of intraband $\sigma_i$ and interband $\sigma_b$  transitions, the closed form approximation to the interband $\tilde{\sigma}_b$ and the sum of intra- and interband conductivity  $\sigma$ are shown. The integral in the interband contribution is evaluated numerically. (a-b) Show the response at low temperature ($T=\qty{5}{K}$) for high ($\gamma=\qty{3e12}{\per\s}$) and low quality ($\gamma=\qty{3e13}{\per\s}$) graphene respectively. The behavior at $T=\qty{300}{K}$ is shown in (c-d) for the same graphene qualities. Note how the interband approximation fails for room temperature, predominantly in the region around \qty{0.4}{eV}, where Pauli blocking sets in.}
\label{fig:kubo_approx}

\end{figure} The Kubo formula divides the contributions to the complex conductivity into intraband and interband transitions Hanson, 2008:

With the intra- and interband contributions described by Mohsin et al., 2015:

Where $f(\xi)$ is the Fermi-Dirac statistic:

According to Hanson, 2008 the contribution of the inter band transitions can be approximated for $k_bT<<|\mu_c|, \hbar\omega$ by:

When compared to the rigorous Kubo formula it is found that the approximation does not consider $T$ leading to significant discrepancy at the onset of Pauli blocking at room temperature (see figure ). In contrast to equations \ref{eq:kubo} the fit by Chang et al., 2014 and the RPA have been found to not incorporate the effect of finite scattering rate on the width of the transition to and from Pauli blocking. As this is assumed to be inaccurate, the Kubo model will be used in this work. The chemical potential $\mu_c$ can be determined from the surface charge carrier concentration $n_s$ as Sorianello et al., 2015:

Where v_f=\qty{9.5e7}{cm/s} is the Fermi velocity in graphene. $n_s$ in turn is determined by electrostatic and chemical doping. As graphene only consists of surface the influence of adsorbates and the dielectric environment is pronounced. For a given capacitor cross-section, it can be expressed in terms of a voltage $V_D$ that has to be applied to force \mu_c=\qty{0}{eV}. $V_D$ is denoted as Dirac voltage, as it is the voltage at which $E_f$ coincides with the DP. The overdrive voltage $V_{od}$ is related to $n_s$ due to the geometrical and quantum capacitance as Sorianello et al., 2015:

Where $C_{ox}$ is the geometrical oxide capacitance per unit area and $K \in \{1,2\}$ is the number of graphene layers used as gate electrodes. The influence of the quantum capacitance is weak for large oxide thickness $d_{ox}$ and $\mu_c$ and can be disregarded commonly (see section \ref{ap:quantum_cap}).

Photonic Modeling¶

To incorporate the optical behavior of graphene as described by its sheet conductivity with numerical solvers for the Maxwell's equations a variety of approaches can be used. These include incorporating the conductivity as a boundary condition or as an effective dielectric constant for a slab of finite thickness. It is also possible to model the influence of graphene as a perturbative term after solving for the mode profile of a given waveguide structure. Two methods exist to identify such perturbative term, which are shown to be equivalent. In the following, the different approaches will be introduced and compared. The sheet conductivity can be expressed by the boundary condition Depine, 2016Koester & Li, 2014:

with $\vec{n}$ the interface normal vector pointing from domain 1 to 2, $\vec{H}^{(i)}$ the magnetic field in domain $i$ and $\vec{E}_{||}$ the in-plane part of the electric field. The most commonly used expression in literature Sorianello et al., 2015Yang et al., 2014Mohsin et al., 2015 is found by inspecting Amperes law in the frequency domain and expanding the sheet conductivity into a bulk conductivity of finite thickness $\delta_{gr}$ Vakil, 2012Depine, 2016:

Incorporating the effect of graphene at simulation time allows for a highly flexible simulation flow able to represent arbitrary geometries. When considering structures that are translation invariant in the direction of propagation it is possible to consider the relatively weak influence of graphene as a first-order perturbative term. It mainly introduces a loss term to the propagation constant, but no significant deviation of the modal field distributions compared to the unperturbed waveguide Midrio et al., 2012. The introduced power attenuation $\alpha$ can be found from perturbation theory (see appendix \ref{ap:perturbation}):

Where $k_0 = 2\pi/\lambda_0$ is the magnitude of the free space wavevector $\vec{k}$ and $\Gamma_\varepsilon$ indicates the overlap of the modal $\vec{E}$-field with the graphene. Its analytical expression can be found in equation \ref{eq:gamma_overlap}. Another expression for the graphene loss can be found by considering the ohmic loss of graphene. It is shown in appendix \ref{ap:comparison_loss_models} that this yields an expression equivalent to the one derived above. For a modulator of a given length $l$ this means that by switching $\Im(\delta \varepsilon_r)$ the power at the output of the graphene coupled waveguide can be switched (see equation \ref{eq:powerdecay}):

Note that the relative power $p\in [0,1]$. The dependence of $p$ on $l$ and $\Gamma$ is expressed by:

Inverse Design with Goos¶

As discussed in section Inverse Photonics Design the photonic inverse design methodology is a capable tool to find devices with optimized performance. It opens the possibility to investigate fundamental limitations to device performance Molesky et al., 2018[p. 667]. In this work, it will be employed to investigate, whether multiple coupled resonators or any other arbitrary distribution of dielectric can overcome the limiting relation between MD and IL described in detail in section Rectlinear Modulator. To do so, the publicly available open source framework Goos is used, which is the updated version of spins-b InvDes Su et al., 2019. For 3 dimensional simulations Goos was used in conjunction with the open source multi GPU FDFD solver maxwell-bhttps://github.com/stanfordnqp/maxwell-b, 2022.^[1] To the knowledge of the author, it is impossible to simulate anisotropic permittivities using maxwell-b, which was regarded as unproblematic for the same reasoning as stated in section Paragraph.

Device Concepts for Inverse Design¶

\begin{figure}[h] \centering \includesvg[width=\textwidth]{goos_no_arrows.svg} \caption{Modulator concepts that allow for optimization using inverse design. The dimensions are not to scale and can be adjusted. The top row shows the transmissive direct modulator (a) and Fabry-Pérot-like modulator (b), while the adjacent row shows the corresponding reflective architectures. (e) Shows a differential modulator. (f-g) rely on cascading repeating pieces of optimized waveguide, that return to the fundamental mode at the interfaces. (g) adds the possibility to optimize the graphene shape on top. A simple subset of this architecture was realized in this work by a grating-like structure of graphene. For less constraint graphene geometries connectivity and resistance checks would need to be implemented.} \label{fig:goos_architectures} \end{figure} A well-conditioned optimization problem has to be stated, including adequate constraints, to allow inverse design to effectively find good solutions. An overview of possible modulator implementations involving inverse design is given in figure . Goos allows for different constraints to be imposed on the device geometry including symmetries and grating-like structure. The gratings are described by their edge positions and can be subjected to strict minimum lengthscale constraints. Restricting the geometry to grating-like structures eliminates a large amount of design flexibility and will likely lead to less-than-ideal device geometries. It further does not represent an easily fabricable structure due to its sharp edges. Alternatively, the dielectric distribution in the design region was described by a cubic spline interpolation between control points spaced 1.5 times the minimum feature size apart, which results in geometries, that adhere to the feature size constraint almost everywhere. Additionally, a sigmoid transfer function is applied to the design to make it more binary. The steepness of the sigmoid and thus the amount of binarization was increased in the course of the optimization. The different functionality realized by the design region is determined by the implemented objective function. Table \ref{tab:obj_functions} gives an overview of the objective functions, that were tested in this work. The objective functions marked as phase-sensitive could be used to implement phase modulators, while all other objective functions only consider variations in amplitude and are thus geared to optimize amplitude modulators.

Separate Simulation of On and Off State¶

All objective functions either directly or indirectly consider the transmitted $\vec{E}$-field at the output port for both states of the modulator. Thus to evaluate the objective function two forward simulations need to be performed. Similarly to determine the gradients two separate adjoint simulations have to be run. This functionality is natively implemented in Goos and has been demonstrated for the optimization of a three-wavelength WDM (de)multiplexer Su et al., 2018. In this work, a projection of the three-dimensional geometry to an effective two-dimensional representation is used to lower the computational burden (see section Paragraph). As a result, the refractive index of the optimized geometry changes between the two modulator states. In a full three-dimensional simulation, the optimized geometry and the changing refractive index are located at different positions in space (displaced by the $z$-coordinate). The optimized geometry and switched refractive index are thus separated into different shape objects. \input{sections/Methodology/obj_functions} Goos handles the iterative optimization based on a problem graph defined prior to the optimization Su et al., 2019. The generating functions defining this problem graph do not easily permit the definition of two parameterized dielectric shapes based on the same geometry-variable but with different dielectric fore and backgrounds. Therefore the respective functions were modified to accept a Goos variable as an input, inhibiting the creation of a separate variable for every parameterized shape. The code modifications are provided in appendix \ref{ap:spins_mod}. The resulting relations between variables and shapes are compared in figure . \begin{figure}[h] \centering \includesvg[width=0.95\textwidth]{modification_param_no_arrows.svg} \caption{Modifications to the generating functions of parametrizations in Goos, to allow for multiple shapes to reference one optimization variable. The dashed arrows indicate the second execution of the generator function, which was modified to accept an existing variable.} \label{fig:goos_variables} \end{figure}

Raman Spectroscopy¶

\begin{figure}[h] \centering \includegraphics[width=0.6\textwidth]{raman_spectrum.pdf} \caption{Prototypical Raman spectrum of single-layer graphene. The notation for the characteristic peaks of graphene used in this work is indicated. The spectrum is generated from a 25x25 map by superimposing the spectra at different positions.} \label{fig:raman_spectrum} \end{figure} Raman spectroscopy is a versatile tool to noninvasively and quickly analyze multiple properties of graphene Ferrari & Basko, 2013. It relies on inelastic scattering between the incident coherent light and graphene, due to an energy exchange with phononic modes of the graphene lattice. Due to the energy exchange from the photon to the phonon (reverse), the scattered light is red (blue) shifted. This shift is called Stokes (anti-Stokes) shift. In the scope of this work, only the Stokes shift will be considered. In a Raman spectrometer, the scattered light is collected, the strong carrier peak resulting from rayleigh scattering is filtered out and the remaining light is decomposed into its spectral components, e.g. using a grating and CCD sensor. The intensities are recorded as a function of the energy shift of the scattered light, typically expressed in terms of the wavenumber shift in $cm^{-1}$ ($1/\lambda_0-1/\lambda$) as seen in figure .

The two Raman peaks of primary interest in this work are denoted as the $G$- and $2\mathrm{D}$-peaks. The $G$-peak originates from a single scattering event with a $E_{2g}$ phonon Ferrari & Basko, 2013. It is dependent on strain and doping in position and amplitude Bendiab et al., 2018.

The doubly resonant $2\mathrm{D}$-peak is caused by an inter-valley scattering process with two phonons of opposite wavevector.

Similar to the $G$-peak the position and amplitude of the $2\mathrm{D}$-peak is dependent on doping and strain, however with a different relative sensitivity, which allows deconvoluting the influence of these perturbations to a first-order. Here we will use the change of base from $G$ and $2\mathrm{D}$ positions to strain and doping given by Bendiab et al. (2018):

Where $x_i$ denotes the $i$-peak position in terms of $cm^{-1}$ relative to unstrained undoped graphene. $p_i$ and $s_i$ express the dependence of the $i$-peak position on hole doping and strain respectively. They are given as Bendiab et al., 2018:

and

It should be noted, however, that in the vicinity of the Dirac point the linear relation between doping and $G$-peak position possibly does not hold as suggested by Pisana et al. (2007), indicating that the $G$-peak position has a minimum for Fermi energies close to the phonon energy. Furthermore, it should be noted, that the $G$-peak stiffens also in the case of electron doping. A distinction of electron and hole doping from the $G$-peak position is therefore not possible.

\begin{figure}[t] \centering \includegraphics[width=\textwidth]{fit_no_fit.pdf} \caption{Comparison of fitted (a) and directly extracted (b) G-peak position. The 5x\qty{5}{\micro\m} map was recorded on a graphene sheet as transferred on $\mathrm{SiO}_2$. By fitting the peak position fine patterns become visible, that resemble the topography of a strained piece of cloth. The directly extracted map suffers strongly from quantization.} \label{fig:fit_no_fit} \end{figure} Intensities and FWHM from the measured Raman spectra Voigt lineshapes (see equation \ref{eq:ap_voigt}) were fitted to the individual peaks by minimizing the squared residuals, to allow for a reliable extraction of the peak positions.

The Voigt lineshape was selected based on Huang et al. (2014). It is noted, that the $G$-peak is experimentally fitted well by a Lorentzian Ferrari & Basko, 2013, and that a modified Lorentzian has been put forward, to accurately fit the $2\mathrm{D}$-peak with a minimal amount of free parameters Froehlicher, 2016. However fitting by Voigt lineshapes should not negatively affect the extracted peak positions or FWHM, and is thus preferred in this work due to the generality of the fitting approach for $G$ and $2\mathrm{D}$-peak.

By extracting the relevant parameters from the fitted curves instead of directly finding the peak position and intensity from the recorded intensities, a quantization of the extracted peak position to the wavenumber sampling points was avoided. Furthermore, the influence of measurement noise towards the peak intensity was reduced, as the intensity at the surrounding wavenumbers was also taken into consideration. Figure compares fitted (a) and directly extracted (b) G-peak positions of as transferred graphene. A significant quantization of the unfitted data to the wavenumber resolution of the Raman system is observed, which was overcome by fitting the resonances by a Voigt lineshape. The fitting reveals an intricate inner structure resembling wrinkles and folds of a unidirectionally strained fabric. Comparing the areas of the 2D and G peak is an alternative method to extract doping (see figure (c)) Tyagi et al., 2022Basko et al., 2009. Contrary to the method discussed previously, the extracted doping levels were sensitive to shifts in the optical focus. Furthermore, the extracted doping levels were anti-correlated to the doping levels from the previous method. The method presented by Bendiab et al., 2018 (a) will be used in the following, as it yields similar results to the method proposed by Lee et al., 2012 augmented by the relation between $\omega_G$ and $\mu_c$ extracted by Froehlicher & Berciaud, 2015 (b).

\begin{figure}[t] \centering \includegraphics[width=\textwidth]{extraction_compare.pdf} \caption{Exemplary comparison between doping extraction methods. White dots indicate areas, where the line shape fit was considered insufficient. The extraction methods are (a) the method by Bendiab et al., 2018 used here, (b) the method by Lee et al., 2012 augmented by Froehlicher & Berciaud, 2015 and (c) the method used in Tyagi et al., 2022 based on the relative area of the G and 2D-peak. (d) shows the histogram of extracted chemical potentials by method.} \label{fig:doping_extraction} \end{figure} Spatially resolved Raman spectra can be recorded when using a confocal microscope. The spatial resolution is determined by the diffraction-limited beam size of the focused laser. In this work, spatially resolved Raman measurements were recorded using the commercially available confocal microscope system Witec Alpha 300R (100X, NA=0.9, \lambda_{ex} = \qty{532}{nm}), which is advertised to have a spatial resolution of \qty{200}{nm} WITec, undefined. It is however unclear, what exactly is specified by this quantity, with an option being FWHM at the beam waist of the excitation beam. Optical characterization of a similar setup with equivalent magnification, numerical aperture and wavelength yielded a standard deviation \sigma_{beam}=\qty{138}{nm} equivalent to a FWHM of d_{beam}=\qty{325}{nm} Froehlicher, 2016[Appendix B].