1. Introduction

Graphene photonics has sparked an enormous interest in the optical community because of graphene’s exceptional optical properties and large light-matter interaction [1]. This is especially the case in integrated photonics due to graphene’s CMOS compatibility, but still, its potential commercial use is not quite extensive [2]. We attribute this to a limited comprehension of the fundamental physical processes occurring in various device configurations, which prevents designing optimal structures comparable to the commercially viable ones of today [3].

We focus on the combination of graphene with silicon photonics, i.e. graphene coated silicon waveguides and integrated structures based on them [4–7]. Graphene based silicon photonic devices are interesting since their fundamental principle of operation is the graphene’s absorption of light which travels parallel to it’s sheet [8, 9]. Even though this has been the main way of increasing graphene’s light matter interaction in many proposed devices, the fundamental absorption properties in this configuration are still to be fully characterized, theoretically or experimentally. In this work, we explore the graphene absorption dependency of the silicon waveguide design, i.e. waveguide’s thickness and width, using a self developed, modified 2D Finite Difference Method (FDM) and fundamental laws of electromagnetism.

In the first part of the paper, we introduce the details of the numerical method, focusing on the novelties to the algorithm we’ve introduced to simulate the graphene layer. The key idea is to use the graphene’s dynamic conductivity and fundamental boundary conditions to electromagnetically simulate absorption. In the latter part of the paper, we introduce our numerical results, i.e. the graphene induced absorption dependency of the waveguide design. We observe, for the first time, absorption peaks in the TM mode curves, as well as the shifts in the dominantly absorbed mode which can be influenced by choosing appropriate waveguide thickness and width. We substantiate our numerical predictions through the comparison with the experimental results from other research groups, which fit perfectly into our model [17, 18, 4].

We also provide qualitative explanations of the numerical, and experimental results, through the comparison of the absorption curves with in-graphene plane optical electric fields, for both modes, and show how our results can be used to optimize the performance of the state of the art devices proposed in this configuration, like optical modulators, photodetectors and polarizers.

2. Graphene modified FDM method for numerical absorption calculation

In this chapter we introduce the numerical method we’ve developed and implemented to calculate the graphene induced absorption. It is fundamentally a 2D FDM method, with modifications which take graphene into consideration. Since we are focusing on silicon photonic integrated circuits, the method is aimed at absorption at the standard λ = 1550nm optical communications wavelength, but it can be extended to other wavelengths as well.

The basic structure of interest, a graphene coated silicon wire waveguide, is schematically presented in Fig. 1(a), with the 2D cross section and the corresponding FDM ”small-grid” example shown in Fig. 1(b) (the actual grid we used in simulations consists of 1000x1000 points, including points bellow the waveguide region which correspond to the silica cladding).

 

Fig. 1 a) Graphene coated silicon wire waveguide structure. b) Cross section of the waveguide, with the depiction of an FDM grid including the waveguide and graphene (small example).

Download Full Size | PDF

Initially, we start with the standard field representations in rectangular electromagnetic media, with axis orientation as in Fig. 1 [10] (vector A corresponds to either the electric E or magnetic H vectors):

The main goal of 2D light propagation simulations is to obtain the propagation direction wavevector β as well as complete mode profiles E(x, y) and H(x, y). By taking graphene into account, β becomes complex, and from its imaginary part we can estimate the total graphene induced absorption [11]. The standard way of obtaining these values is by inserting the estimated field profiles (Eq. (1)) into the Maxwell’s curl equations, and solving the obtained system of differential equations numerically [10].

Most widely used numerical methods for solving this problem are the FEM (Finite Element Method) or, in the general case, FDTD (Finite Difference Time Domain), but both of them require approximating the graphene as an additional thin layer with an effective refractive index [6]. Since graphene is one-atom thick, we believe it is the FDM method which is more suitable, as graphene can be incorporated through its fundamental optical properties and boundary conditions, as will be shown in this chapter.

In the FDM method, the transverse plane of the waveguide is represented in a rectangular grid of equidistant points (Fig. 1(b)), where each point corresponds to a value of the structure’s relative permitivity, and it’s field value. In this numerical method, fields are expressed as unknown algebraic vectors with elements corresponding to points in the grid, and differentials as matrix operators multiplying the respected vectors, effectively introducing the field value differences of neighboring points in the grid divided by the grid resolution - the finite differences.

For additional clarity, and to introduce the way in which we uniquely modify the FDM method to include graphene into consideration, we focus on one of the magnetic field curl equations, which contains fields parallel to the graphene boundary:

When we apply the FDM method on the Eq. (2) to obtain it’s algebraic form, it becomes:

In Eq. (3), ε_r is a diagonal matrix, which elements correspond to the relative permitivity of the points in the grid - the structure presented in Fig. (1). Additionally, h_x and h_y are unknown algebraic magnetic field vectors and e_z is the unknown algebraic z–direction electric field vector. D_x and D_y are differential matrix operators which induce the differences of neigh-boring field values in the grid, normalized by the grid resolution. When all the curl differential equations are represented in this way, we obtain a system linear algebra equations which can be simplified into a regular eigenvalue problem where the eigenvalue is the desired β and the eigenvectors are the field vectors. These types of problems can easily be solved using standard numerical libraries available in most programming languages, but before actual implementation most of the equations have to be further normalized for consistency [12].

Graphene is fully described electromagnetically through its surface dynamic conductivity, which contains parts corresponding to inter, and intra-band carrier transitions [13, 14]. In the case of our target wavelength of λ = 1550nm, when the photon energy is 0.8eV, the dominant transitions are interband, and conductivity is given as [13]:

Additionally, as shown in devices employing this type of structure [4, 5, 6], the doping of graphene arising from the contact with Si at room temperature is negligible (|E_F| < ħω/2), and the dynamic conductivity is approximately σ ≈ e²/4ħ which is the value we use in simulations.

In our numerical approach we utilize the Eq. (4) to modify the standard FDM method and take graphene into account. We accomplish this by employing the magnetic field boundary condition, which states that the difference of the tangential magnetic fields on the opposite sides of the boundary has to be proportional to the tangential electric field at the boundary, with the factor of proportionality being the surface conductivity:

When graphene lies at the boundary, this surface conductivity directly corresponds to graphene’s surface dynamic conductivity given in Eq. (4), and is an ideal way of taking graphene into consideration.

However, numerical methods such as the standard FDM, FEM, and others, don’t take boundary conditions into account explicitly, as they are intrinsically satisfied within the grid of choice. Since we want to explicitly take the boundary condition into account, we return to the fundamental properties of the FDM method and use them to our advantage: Eq. (3) is a linear algebra equation, so it actually presents a system of equations itself, one equation for every point in the grid. A closer look at the equation which deals with fields around the graphene boundary (as per the ”small-grid” example shown in Fig. 1(b)) gives us:

At his point, we introduce our main modification to the algorithm: We can extend this equation with an ”imaginary” grid point, lying in between the two points in the grid, which corresponds to the boundary:

It can be observed that this extension makes no difference to the overall equation, as if there were nothing at the boundary, the two additional terms corresponding to the boundary point would cancel each other out (as per the boundary condition). However, when there is conductivity at the boundary, as is the case in the presence of graphene, the additional term becomes:

The validity of this approach can further be justified through the following reasoning: We can assume that we add ”imaginary” points to all the equations of Eq. (3) to preserve the consistency of the method. However, they cancel each other out in every instance, except for the one at the graphene boundary.

The final step is to correlate the imaginary point with the one in the existing grid, and we can do this with the following approximation:

Approximation in Eq. (9) is justified by the fact that the z field is the tangential electric field at the boundary, and the difference in it’s value between the point in the regular grid, and the ”imaginary” point which is only a half of the grid resolution away, is negligible. Finally we obtain the equation which takes graphene into consideration:

This process is repeated for every equation corresponding to the points on the opposite sides of the graphene boundary, and Eq. (3) is modified accordingly. Using the modified equation we derive the eigenvalue problem, through which solving we obtain the complex propagation direction wavevector β (eigenvalue), and the mode profiles (eigenvectors). Since the power of the propagating light is proportional to the absolute value of the field squared, the imaginary part of β corresponds to the attenuation, and we can estimate the graphene absorption (P ∼ exp(−2Im(β)z) [11]). We also use the quasi-vectorial approximation, in which we divide the general problem into two, one corresponding to the TE modes, and one to the TM modes.

The previous steps are the main reason we used the FDM method. It provided us with the perfect framework to easily implement the boundary conditions and estimate graphene’s absorption through fundamental EM modeling. This approach is unique, to the best of our knowledge, although the final result is similar to some of the approaches implemented by other groups [15, 16]. We believe this simple way of characterizing graphene could give better insights into the light matter interaction of graphene coated silicon waeguides as it is much more focused on fundamental electromagnetic theory and basic descriptions of graphene.

3. Numerical absorption results

The method described in the previous chapter was implemented in the Python programming language. We used standard Python linear algebra libraries ”NumPy” and ”SciPy” to solve the eigenvalue problem and obtain the complex propagation wavevectors and the mode profiles. The grid consisted of 1000x1000 points, equally spaced, with the grid resolution Δx = Δy = 2nm. Since we are not adding additional layers to simulate graphene, and rather take it into account through the boundary conditions, we believe this grid resolution is sufficient.

Using our method we estimated the graphene induced absorption with respect to different waveguide thicknesses and widths, and the main results are presented in Figs. 2 and 3. All dimensions were chosen to correspond to standard ones in silicon photonic chips, with an emphasis on single mode operation (for smaller waveguide widths).

 

Fig. 2 Graphene induced absorption dependency of the waveguide thickness d, for different waveguide widths w. Solid lines – Original numerical results; Triangle points – Experimentally reported results [17, 18, 4].

Download Full Size | PDF

 

Fig. 3 Graphene induced absorption dependency of the waveguide width w for different waveguide thicknesses d.

Download Full Size | PDF

For the first time, we theoretically predict a peak in the TM mode absorption curve (Fig. 2) at thicknesses d ≈ 240nm, as well as the shift in the dominantly absorbed mode at around d = 200nm waveguide thickness. Similar results can be observed in Fig. 3 with shifts in the waveguide width, for smaller thickness values. We believe these novel results could be crucial in optimizing devices which require large modulation depths, like modulators or photodetectors, based on the absorption peaks, or even creating new ones based on the design-dependent mode absorption dominance, like optical filters or polarizers.

In Fig. 2, we also included the case of the slab waveguide (w → ∞). This result was not obtained using our modified FDM method, rather we obtained it semi-analytically by calculating, and solving, the dispersion relations of the slab waveguide based structure, calculated as:

We have also compared our numerical results with previously reported experimental results (Fig. 2), and observed a very high order of correspondence. This is especially interesting in the case of the TM mode absorption curve, where the existence of the peak can clearly be anticipated from experiments, as can the reversal of the dominantly absorbed mode, both of which we are first to theoretically predict, to the best of our knowledge. We provide an extensive qualitative explanation of our findings in the following chapter.

4. Discussion

The novel numerical results from the previous chapter, which have not been reported previously by using standard methods or simulation tools, can be theoretically explained through the fundamental property of graphene, which states that it interacts with the electric fields parallel to it’s sheet, i.e. fields which lie in the graphene plane [4]. In the case of the TE mode, the field parallel to graphene is the standard, and the only dominant, electric field corresponding to the x direction in Fig. 1. However, in the TM mode case the only electric field which lies in the graphene plane is the propagation direction, longitudinal, z–electric field. The electric field usually used to describe the TM mode is the y direction electric field, but the y–field is perpendicular to the graphene sheet and doesn’t influence the graphene’s light-matter interaction.

The fact that graphene interacts with the longitudinal electric field is largely overlooked for the TM mode in literature, but could prove to be a way for explaining the remaining unknowns in experimental results. The mode profiles of the electric fields in the graphene plane, for both modes, are plotted in Fig. 4(a).

 

Fig. 4 a) Profiles of the electric fields parallel to the graphene sheet in the case of TE and TM modes. b) Dependency of the electric fields parallel to the graphene sheet of the waveguide thickness; Inset: Point used in calculations for Fig. 4(b) marked in red.

Download Full Size | PDF

To further characterize the interaction of graphene with the in-plane electric fields for both modes, we plotted their dependency of the waveguide thickness, as presented in Fig. 4(b). By comparing the results in Fig. 2 with Fig. 4(b), we can observe a very high correlation between the intensity of the in-plane fields and the total graphene induced absorption, which further substantiates our claims.

Graphene’s interaction with the in-plain fields is intuitive for the TE mode, for which the in plane, x field, also corresponds to the mode’s power distribution inside the waveguide by the Pointing’s theorem. On the other hand, TM mode’s the longitudinal, z electric field has no qualitative interpretation, and the field which corresponds to the power distribution is actually the x direction magnetic field, which doesn’t interact with graphene. The explanation that the longitudinal electric field interacts with graphene and induces absorption is slightly counter-intuitive, but it fully describes the results presented in Figs. 2 and 3 and can be qualitatively explained by the fact that if we were to characterize absorption (which is a quantum-mechanical property of materials) quantum mechanically, and through the Fermi’s Golden rule, it is this field which would be used as the perturbation to the kinetic operator for absorption calculation.

5. Conclusion

In this paper, we presented a novel way to simulate absorption in graphene coated silicon wire waveguides, and we observed reversals of the dominantly absorbed mode with waveguide design variation, as well as absorption peaks in the TM absorption curves for the first time, to the best of our knowledge. We sincerely hope our results will spark further interest in the fields of graphene enhanced silicon photonics and optoelectronics.