Appendix

Silicon Photonics¶

Splitters and Couplers\label{ap:splitters}¶

Directional Coupler¶

In many cases it is necessary to couple a defined amount of power from one waveguide. One versatile approach to accomplish this coupling is to use the overlap between the evanescent field of one waveguide mode with the mode of a closely spaced waveguide. The change of base used to decouple the propagation equations in CMT can also be understood as the orthogonal modes of the perturbed system of two waveguides in close proximity. As these so called supermodes propagate with different phase velocities they alternatingly interfere constructively in one or the other waveguide. The first constructive interference in the non-excited waveguide is located at the so called cross-over length $L_x$ defined by Chrostowski & Hochberg, 2015:

Where $\Delta n$ is the difference in effective refractive index of the supermodes.

Y-junction and Multi Mode Interferometer (MMI)¶

In a directional coupler the ration of coupled and transmitted power is sensitive to deviations from the modeled device geometry and refractive index distribution (i.e. the instantaneous coupling coefficient $\kappa^\prime$ changes exponentially with the gap between the waveguides) Romero-García et al., 2013. In case that specific splitting ratios (e.g. $50:50$) are desired more reliable alternatives exist in the form of Y-junctions and MMIs. Y-junctions consist of a single waveguide splitting in two. The area around the junction is typically optimized in shape to introduce little excess loss. Due to the devices symmetry close to symmetric power coupling is guaranteed Zhang et al. (2013)De la Cruz-Coronado et al. (2020).

Resonators¶

Photonic Crystal and Nanobeam Cavities\label{ap:crystal}¶

Bragg reflectors are the one-dimensional variant of photonic crystals, which can prohibit the propagation of specific wavelengths by a periodic arrangement of objects with different refractive indices. Similar to actual crystals bandgaps can form, and by carefully introducing defects to the lattice, light can be guided and resonators can be formed Joannopoulos, 2008. A special type of photonic crystal cavity is formed by etching a shifted periodic pattern of holes into a rectangular beam of material. The light confinement in such nanobeam cavities can be enhanced by introducing a slot in the elongated central section. The resulting slotted nanobeam cavities exhibit simultaneously a relatively high quality factor and ultralow mode volume Yu et al. (2011)Ryckman & Weiss (2012)Seidler et al. (2013).

Optical Bandwidth¶

In high-$Q$ resonant modulators the photon lifetime in the cavity can be long, limiting the maximum EOBW of such systems Lin et al., 2013. Schemes exist to mitigate this limitation, e.g. by modulating the resonator coupling Sacher et al., 2013. As the introduction of graphene to the resonator implies additional cavity loss, such modulators typically do not exhibit very high intrinsic $Q$ Midrio et al., 2012. As discussed the resonant modulator maintains a significant advantage in ROMA only for a fraction of approximately \qty{20}{%} of the resonance FWHM. Thus an amplified FWHM is required to maintain the modulation performance over the full channel bandwidth. The FWHM is evaluated analytically for ring resonators of different radius $R$ and coupler transmission $\tau$ in the critical coupling state (figure Figure 1). It is observed that FWHM of more than \qty{100}{GHz} are reached for relevant $\tau$ and $R$ (see section Minimizing Graphene Footprint).

Inverse Design¶

Adjoint Method\label{ap:adjoint}¶

Navigating large parameter spaces to find optimal solutions benefits from the possibility to evaluate the gradient of the objective function with respect to all parameters Dababneh et al., 2018. The adjoint method is capable to provide such gradients from only two Maxwell simulation runs independent of the number of free parameters. Here only a brief introduction to the adjoint method will be given. Please find a more detailed description in Piggott, 2018 and Hammond, 2022. The Maxwell's equations for a given dielectric distribution and excitation are a system of linear PDEs and can be expressed in matrix notation:

With $A$ the system matrix depending on the parametrization $p$. $\vec{E}$ is the sought \efield and b an excitation term depending on the sources present. Solving for $\vec{E}$ is denoted the forward pass. Let $o(\vec{E})$ be an objective function which eventually should be minimized over $p$. Separating $p$ into its individual parameters $p_i$ one finds using the sum and product rule Piggott, 2018:

with $\vec{E}_{a d j}$ defined by:

These adjoint fields can be found by solving a matrix problem analogous to the forward solve. Using this method the gradients towards all parameters are found using a single matrix solve, which makes it highly efficient compared to a finite differences approach for example Piggott, 2018.

Parametrization¶

The index distribution in the design region has to be described by a finite set of parameters, to use the inverse design approach. Multiple such parametrizations have been put forward over the last decade for planar integrated photonic devices. These include levelset- Piggott et al., 2020 and density based Hammond, 2022 parametrizations. Recently a parametrization based on a conditional generator was proposed, which allows meeting the fabrication constraints in every step of the optimization Schubert et al., 2022.

Bandwidth Scaling with Graphene Length\label{ap:bw_experimental}¶

The bandwidth limited by the RC-time constant can in theory be determined from a cross-sectional model disregarding the modulator length Cheng et al., 2020. Experimental realizations of graphene modulators have however shown a length dependence of the bandwidth Hu et al. (2016)Wu et al. (2021). A possible explanation for this discrepancy might be the method by which the devices are driven. Typically the modulation voltage is applied via a transmission line or waveguide matched to \qty{50}{\ohm}. These are in the order of the expected total series resistance for the modulators. The RC time constant can be decomposed into a constant contribution from the internal resistance, which is length invariant, and a contribution from the wave impedance $Z_0$. This effect is illustrated in figure Figure 2 for a modulator with the following specifications:

Note, that if the given explanation for the bandwidth scaling is appropriate, a suitable low-impedance driver could be integrated to mitigate this effect. The integration of specialty drivers poses additional challenges, especially when considering back-end-of-line integration with any given chip technology. It is thus beneficial to scale down the graphene capacitance to reach high bandwidth without the need for specialty drivers.

Switching Energy\label{ap:switching_energy}¶

It should be noted, that under most circumstances the relation given by equation \ref{eq:e_per_bit} holds even with an applied bias voltage contrary to the formulation used in Midrio et al., 2012:

With $Q_{on}$ ($Q_{off}$) the total charge present on one capacitor plate in the on (off) state. Using the defining equation of capacitance, one finds $E_{bit,midrio}$ to be:

This contradiction to equation \ref{eq:e_per_bit} stems from the fact, that the energy transferred back from the capacitor to the supply Miller, 2012 is not considered by Midrio et al., 2012.

Raman Fitting - Voigt lineshape¶

Where $G$ is a Gaussian and $L$ is a Lorentzian lineshape, and $\sigma$ determines the width of the Gaussian, while $\gamma$ is the HWHM of the Lorentzian. $*$ denotes the convolution operation. To perform the numerical fit the following expression was used https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.voigt_profile.html, undefined:

Where $w(z)$ is the Faddeeva-function, for which a numerical approximation exists https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.wofz.html, undefined, and z is defined as:

The FWHM of the Voigt fit $\Gamma_\mathrm{V}$ can be determined from the FWHM of the Gaussian ($\Gamma_\mathrm{G}$) and the Lorentzian ($\Gamma_\mathrm{L}$) contributions Olivero & Longbothum, 1977:

With

And

Graphene Modelling¶

Quantum Capacitance \label{ap:quantum_cap}¶

Graphene has a linear dispersion around the $K$ and $K^\prime$ points, denoted here as DP Fang et al., 2007:

Describing the energy of the electron in the conduction (CB) and valence bands (VB) relative to the DP, having a relative wavevector $\vec{k}$ with respect to the DP. v_f\approx\qty{1e6}{\m \per \s} denotes the Fermi velocity, which can be described as an effective speed of light experienced by the carriers in graphene, having zero mass at rest Novoselov et al., 2005. A potentially important quantity to describe the electrostatic doping control of a low dimensional system is the quantum capacitance defined as:

Where $e$ is the elementary charge, $Q$ is the net charge in the graphene layer and $\mu_c$ is the chemical potential. The size of graphenes quantum capacitance can be found by considering its linear dispersion relation, in conjunction with degeneracy in spin and valley to find the DOS. The density of electrons and holes can then be found by integrating the DOS filled by carriers according to Fermi-Dirac statistics. Using the defining equation \ref{eq:quantum_cap} one finds Fang et al., 2007:

Which can be approximated in the regime where the chemical potential is much larger than the thermal energyFang et al., 2007:

Not that $C_Q$ should stay positive for all $\mu_c$. In the following $|\cdot|$ will be introduced to ensure that. Starting from this equation for the quantum capacitance of graphene an equation directly relating the surface charge density $n_s$ and chemical potential will be derived. The resulting formula (equation \ref{eq:n_s_mu_c}) is broadly used in literature Wang et al. (2008)Das et al. (2008)Sorianello et al. (2015), however an original source was not cited and could not be determined by the author. The quantum capacitance can be understood to describe the additional energy that needs to be put into a gate insulator graphene capacitor to store some amount of charge besides the energy stored in the displacement field $\vec{D}$. This additional energy arises due to the gradual filling of states of the graphene DOS. In an equivalent electrical representation the quantum capacitance is placed in series to the geometrical capacitance, with the channel voltage $V_{ch} = \mu_c/e$ present at the intermediate node.

It is important to note, that the usual relation between voltage, capacitance, and charge does not hold here due to the voltage-dependent quantum capacitance $C_Q = C_Q(V_{ch})$. As a consequence to obtain the charge the integral expression has to be used:

With the surface charge $Q=e n_s$ and the relation $V_{ch} = \mu_c/e$ this yields:

Perturbative Graphene Modeling¶

Two approaches to perturbative graphene modeling will be introduced and compared.

Perturbation Theory\label{ap:perturbation}¶

The change in $n_{eff}$ due to a small (and possibly local) change in permittivity $\delta \varepsilon_r$ can be described using perturbation theory by Yariv & Yeh, 2007:

Where $Z_0=\mu_0 c$ is the wave impedance in vacuum, $\vec{E}$ and $\vec{H}$ are the modal electric and magnetic fields respectively. $\delta \varepsilon_{r}$ is the dielectric perturbation tensor introduced by the graphene sheet. Assuming a homogeneous $\delta \varepsilon_{r}$ in the directions parallel to the graphene sheet (0 otherwise) across the graphene sheet(s) one finds $\delta n_{eff} = \Gamma_\varepsilon \delta \varepsilon_{r}$ with a proportionality factor of:

The graphene absorption is assumed to be the dominant contribution to the waveguide loss. The real part of the propagation constant $\gamma=-\alpha/2+j\beta$, describing the field decay along the direction of propagation, is thus determined by:

Ohmic Dissipation¶

The time-averaged work done by an electric field applied to an infinitesimal sheet conductor of area $dA$ and complex conductivity $\sigma$ is expressed as Depine, 2016:

Due to translation invariance in the direction of propagation one finds $\alpha$ from the relation between the total local power dissipation $W_{tot}(z)$ and the overall average power propagating in the waveguide $P(z)$:

For $\sigma$ constant across the graphene sheet on ends with:

Similar to $\Gamma_\varepsilon$ the proportionality constant $\alpha = \Gamma_\sigma \Re(\sigma)$ is defined as:

Mesh Convergence\label{ap:fde_convergence}¶

To verify the different techniques presented in Photonic Modeling, they were used to calculate the absorption loss in graphene-covered silicon waveguides of different widths. It was observed, that the different methods yielded significantly different results, under the condition, that the mesh resolution was not sufficient to resolve the \qty{0.35}{nm} thick graphene layers (see fig. Figure 3). For a mesh with less than \num{50} cells in the \qty{15}{nm} high region surrounding the graphene sheets, the effective permittivity method fails rapidly as indicated in figure Figure 4, while no significant deviation was observed for the other methods, as long as the global mode profile was mostly unaffected.

Comparison\label{ap:comparison_loss_models}¶

In the following, it will be shown, that the two expressions found above are equivalent. Starting with the observation, that the denominator in equation \ref{eq:gamma_overlap} contains the surface integral over the time-averaged pointing vector $\vec{S}=1/2\vec{E} \times \vec{H}^{*}$. The integral thus represents the power flow perpendicular through integration plane $2P(z)$. Thus one finds:

Comparison to equation \ref{eq:ohmic} considering that the thickness of the integration region is fixed to the effective graphene thickness $d_{gr}$:

With the expressions for $k_0 = \omega \sqrt{\varepsilon_0 \mu_0}$ and $Z_0 = \sqrt{\mu_0/\varepsilon_0}$ one finds:

Which is consistent with the expression in equation \ref{eq:eff_epsilon}. The two perturbative expressions for graphene absorption have thus been shown to be equivalent.

Locality of Doping and Strain\label{sec:acf}¶

Spatial inhomogeneity in the local doping of graphene can lead to a reduction of the effective absorption switching Mohsin et al., 2014 (more details in section Device Size with locally correlated Doping). It is thus important to quantify these variations and their topography. The variations in doping were extracted from the $G$- and $2D$-peak positions, as described in Raman Spectroscopy, for each position on a grid rastered by the confocal microscope. From these maps the mean $\mu_{dop}$ and standard deviation $\sigma_{dop}$ in doping, and its PDF were determined. To describe the topography of the doping surfaces the two-dimensional discrete (normalized) ACF of the doping/Fermi energy surfaces ($A$, $B$) was determined by:

Where $\mathrm{CCF}_{A,B}$ is the crosscorrelation function of $A$ and $B$ defined as (generalized from Zhao et al., 2006):

Where $N=(2w+1)^2$ is the total number of pixels, and $w$ is the number of pixels from the center of the map to the edge. Direct evaluation of equation \ref{eq:crosscorr} for maps with large ($>1025\times1025$) pixel count is computationally intensive. It can however be sped up by considering the similarities to the convolution and leveraging the frequency domain formulation of the latter. Surface topographic features that exhibit a property analogous to the Markov property have been shown to yield an exponential autocorrelation function decaying radially with Agterberg, 1970:

Where $\rho$ is the radial displacement from the origin and $L_{corr}$ is defined as the correlation length of the surface. To obtain $\mathrm{ACF}_{A}(\rho)$ in this work, $\mathrm{ACF}_{A}(x,y)$ was numerically averaged radially. Surfaces with exponential ACF have been observed and analyzed extensively in literature Zadehgol, 2021. It is of interest, how the influence of the doping variations changes for different device sizes with respect to the correlation length. To evaluate the switching degradation due to local doping variations independent of other degrading factors a Heavyside switching behavior was assumed. The influence of doping variations, measured at a certain gate voltage (e.g. V_g=\qty{0}{V}) was modeled by the following first-order approximation: The doping was expressed in terms of a local Dirac voltage $V_D$ which was assumed to be constant with changes in chemical potential $\mu_c$ for every point on the graphene sheet. It has been comprehensively shown in literature, that this approximation does not fully hold i.e. due to border traps charging and discharging, which alters the potential landscape seen by the graphene Singh & Gupta (2018)Knobloch et al. (2022). However, modeling such effects would introduce a large amount of complexity and is out of the scope of this work. From the local $V_D$ and a global overdrive voltage the local $\mu_c$ can be found according to equation \ref{eq:v_od}, disregarding the quantum capacitance:

and with equation \ref{eq:mu_c}:

To generate random signals which exhibit a specified autocorrelative behavior the Wiener-Khinchin theorem was leveraged, stating that for a given random signal $s(t)$ the Fourier transform of the ACF is the PSD of $s(t)$ Ohm & Lüke, 2014. Where the PSD is the squared magnitude of the FT of $s(t)$ ($S(f)$)^[1]. A spectrum $S(f)$ can be generated to have a given PSD by taking the square root of the PSD and adding a uniformly sampled random phase $\phi(f)$:

$S(f)$ can then be transformed to $s(t)$ using the inverse FT. This approach can be generalized to discrete time signals using the DFT and inverse DFT, commonly implemented numerically as a FFT. For one-dimensional signals with an exponential ACF, the Wiener-Khinchin theorem yields a Lorentzian PSD, as is typically observed in Wiener processes (e.g. Brownian motion, laser phase noise)

Switching Degradation due to Doping Variations\label{ap:dop_var}¶

The following assumptions are made: The doping variations are normally distributed with a mean $\mu_{dop}$ and standard deviation $\sigma_{dop}$. The shape of the doping distribution is unaffected by the electrostatic modulation, which only shifts $\mu_{dop}$:

Where $g$ is the PDF of the doping distribution. Moreover, the ideal modulation behavior is modeled by a shifted and mirrored Heavyside function $\Theta$:

The absorption $\alpha$ is used here to avoid confusion between the standard deviation and surface conductivity. It is noted that the derivation also holds for the conductivity at optical frequencies, which is proportional to $\alpha$. $\alpha_0$ denotes the absorption without doping. The expectation value for $\alpha$ evaluates to:

The approximation holds for $\mu_c$ being significantly positive. For negative $\mu_c$ the derived result holds due to symmetry. As the absolute positions of $\mu_c$ are irrelevant in this analysis the evolution of $\bar{\alpha}$ with $\mu_{dop}$ can be described by

Evaluating the positions at which $\tilde{\alpha}/\alpha_0$ reaches $p$ and $(1-p)$ on finds:

Resulting in $\Delta\mu_{c,80\

Noise Performance Modulator¶

Here the expressions for $\delta P_{l}$ will be derived for a link limited by thermal noise, shot noise, and a combination of RIN and shot noise.

Thermal Noise\label{ap:thermal_noise}¶

As discussed in section Sources of noise in optical links thermal noise is independent of the optical power at the detector:

Thus equation \ref{eq:q-link} becomes:

Where $Q_{link}$ is abbreviated as $Q$ in this section. Which can be expressed in terms of the input power to the modulator $P_{in}$ and the scalar $md=P_1/P_0$ and $il=P_1/P_{in}$:

Solving for $P_{in}$ and comparing it to the power required for an ideal modulator ($il=1$ and $md\rightarrow\infty$) one finds:

Defining the ROMA. $\delta P_{l}$ indicates the factor by which the laser power has to be increased due to the nonideal modulator $P_{l,nonideal}=\delta P_{l} \cdot P_{l,ideal}$.

Using equation \ref{eq:md_il_f} in its scalar form one finds $il/md = il^f$ and thus:

Which is shown for different $f$ in figure . The position of the maximum is found as:

the maximum achievable ROMA thus evaluates to:

Shot Noise\label{ap:shot_noise}¶

Similarly in the shot noise limit, one finds from equation \ref{eq:q-link} that:

Again solving for $P_{in}$ and comparing the result to the ideal modulator one thus finds:

RIN and thermal noise\label{ap:rin_power}¶

For a combination of thermal noise and RIN, the link quality factor is determined by:

Evaluating \ref{eq:rin_p1} for $\delta P_{laser}$ is performed as follows. The equation is brought into the form of equation \ref{eq:overview_eq} by the following substitutions:

Leading to:

This equation is solved for $x$ by squaring both sides collecting all terms but the remaining square root on one side and squaring again. By eliminating trivial solutions, that were introduced by the squaring one finds:

Introducing the values from before one finds:

from there the comparison to an ideal modulator with infinite modulation depth yields:

Singly Coupled Resonator¶

HWHM \label{ap:hwhm}¶

The HWHM in terms of phase can be found from \ref{eq:t_allpass} to be:

which was found to be approximately $(1-\tau^2)/\tau$ for $\tau$ close to 1 (see fig. ). This approximation is used to normalize $\Delta \phi$ in the following. According to \ref{eq:hwhm} the HWHM can only be determined for $\tau>\tau_{crit} = \sqrt{3-\sqrt{8}}$, as $T=0.5$ is not reached else. \begin{figure}[h] \centering \begin{subfigure}[b]{.5\linewidth} \includesvg[width=\textwidth]{HWHM_phase.svg} \caption{} \end{subfigure} \begin{subfigure}[b]{.3\linewidth} \includesvg[width=\textwidth]{HWHM_tau.svg} \caption{} \end{subfigure} \caption{The HWHM of the resonance dip in terms of $\phi$ for different coupler transmission $\tau$ (a) and its dependence on $\tau$ (b). An approximation to the dependence is indicated by the dashed line.} \label{fig:HWHM} \end{figure}

Offresonant IL¶

Figure Figure 5 shows the IL in the offresonant case.

Q/V \label{ap:qv_eta}¶

Motivated by section Arbitrarily shaped Cavities lets assume:

Where $l_gr$ is the required graphene length and $\eta$ is some proportionality constant. Solving for eta and inserting the expressions for $Q$, $V$ and $l_gr$ one finds:

Undefined control sequence: \unit at position 69: …alpha_{trans} [\̲u̲n̲i̲t̲{dB \per \m}]} …

\eta = \frac{-20 \pi n_g}{A \lambda_{res}} \frac{d}{\alpha_{trans} [\unit{dB \per \m}]} \frac{\log_{10}(a\tau)\sqrt{a\tau}}{1-a\tau}

In the limit of large round trip transmission the latter part converges to:

The expression for $\eta$ thus becomes:

Undefined control sequence: \unit at position 85: …\alpha_{trans}[\̲u̲n̲i̲t̲{dB \per \m}]} …

\lim_{a\tau\to1} \eta = \frac{20 \pi n_g \log_{10}(e)}{A\lambda_{res}\alpha_{trans}[\unit{dB \per \m}]} \left[ \unit{m^{-2}}\right]

Using the approximate values presented in table \ref{tab:eta_parameters} eta is evaluated to be \eta\approx\qty{2e17}{m^{-2}}. \begin{table}[] \centering \begin{tabular}{c|cc} Property & Quantity & Unit \ \hline $d$ & 5 & - \ $n_g$ & 4.4 & - \ $A$ & 0.2 & \unit{\micro\m^2} \ $\lambda_{res}$ & 1550 & nm \ \multicolumn{1}{l|}{$\alpha_{trans}$} & \multicolumn{1}{l}{\num{1e4}} & \multicolumn{1}{l}{\unit{dB\per\m}} \end{tabular} \caption{Approximate values used to evaluate \eta\approx\qty{2e17}{m^-2}. $n_g$ and $A$ are based on FDE simulations of a \qtyproduct{220x500}{nm} SOI waveguide. The other quantities are chosen based on analysis in the main body of this work.} \label{tab:eta_parameters} \end{table}

Domination factor¶

\label{ap:domination} Given

and

note that

The reduced switching factor $f^\prime$ is described by:

with

and thus

and with

one finds

Reflective Modulator Converter \label{ap:refl_conversion}¶

The photonic circuit shown in Figure 6 is suitable to extract the output light from a reflective modulator. It relies on a $\pi/2$ phase shift applied to the light that couples from one waveguide to the other in a 50:50 DC. As the light transits the DC twice light ending up in the upper left waveguide has experienced this phase shift exactly once regardless of the path it took (upper or lower reflector). The light exiting at the lower left has either experienced the phase shift twice or not at all. The resulting $\pi$ phase difference leads to destructive interference in the lower arm, with all light interfering constructively in the upper arm. The configuration is essentially a Michelson interferometer that has arms of equal length.

Goos Modification\label{ap:spins_mod}¶

The following example demonstrates how the functions generating the parametrizations were modified. The original function definition can be found here.

from discretize import TensorMesh

def compat_param(...
                 var = None, 
                 #added to allow using existing variable
                 **kwargs) 
    -> Tuple[goos.Variable, OldStyleParametrization]:
    '''
    ...
    '''
    _, var_shape, lower_bound, upper_bound = create_old_param(
        param, extents, pixel_size)
        
    if not var: #if exists use given variable
        var = goos.Variable(initializer(var_shape),
                            name=var_name,
                            lower_bounds=lower_bound,
                            upper_bounds=upper_bound)
                            
    return var, OldStyleParametrization(var=var,
                                        ...
                                        **kwargs)

Modifications for BEOL Integration\label{ap:beol}¶

The planned BEOL integration poses multiple challenges to resonant graphene modulators. Due to the lacking availability of crystalline silicon in the BEOL, the waveguide material has to be replaced, resulting in the following two challenges: Most of the available dielectrics transparent at telecommunication wavelengths are non-conductive. Therefore the use of a GIG capacitor becomes mandatory. Therefore the model for the influence of doping variations should be augmented to consider correlated and uncorrelated variations in doping of the two separate sheets. Similarly, it becomes necessary to align the relative position of the idle doping level. Moreover, with the exception of amorphous silicon, the available dielectrics have drastically lower refractive index compared to silicon, resulting in a decreased index contrast and therefore confinement. The lower confinement bears two problems. It leads to an increased overlap with pads and access graphene, lowering $d$. To compensate the overlap and pad separation have to be increased, penalizing the bandwidth. Furthermore, the lowered contrast reduces the ability to guide light, which should increase the mode volume and lower $Q/V$. The exact influence on the required graphene area should be modeled in further detail going forward. The larger resonator dimensions will also pose additional computational cost on the inverse design algorithm, which can be partially compensated by less demanding mesh requirements.

Miscellaneous¶

\input{sections/Methodology/util}