COMSOL 6.4 - Graphene Metamaterial Perfect Absorber

Graphene, carbon atoms arranged in a two-dimensional hexagonal lattice, has sparked tremendous research and application interests since its experimental discovery about two decades ago. Besides being ultra-thin, this magical material exhibits a plethora of interesting properties, including high electrical and thermal conductivities, high elasticity, high mechanical strength, and so on. Among various applications, a promising field is graphene-based electro-optical devices, such as photodetectors, photodiodes, and metamaterials. An additional desirable trait of graphene is that its optical response can be actively controlled by changing its Fermi energy via electrical gating. In this model, we first demonstrate how to compute the optical conductivity of graphene using the Kubo formula. The computed conductivity is then used to model a graphene-based THz metamaterial absorber (Figure 1). Due to the atomic thickness of graphene, explicit volumetric modeling of it would be computationally expensive. We show that this can easily be avoided by using the Transition Boundary Condition (TBC) to consider graphene as a 2D surface.

Figure 1: Schematic of the graphene-based THz metamaterial absorber.

Model Definition

Both the electronic intraband transitions and interband transitions contribute to the conductivity of graphene. Using the Kubo formula, it has been shown that the intraband and interband contributions are given by

(1)

(2)

where kB is the Boltzmann constant,

is the reduced Planck constant, e is the electron charge, T is the temperature, eF is the Fermi energy, τ is the relaxation time, and ω = 2πf is the angular frequency. The function H(Ω) is given by

(3)

Finally, the total 2D sheet conductivity of graphene is given by σ = σintra + σinter. In this model, we consider T = 300 K and τ = 10−13 s. The integral in σinter can be performed using the built-in integrate() operator. Note that compared to Ref. 1, j in the above equations is replaced with −j due to the phase convention in COMSOL.

For low frequencies, in the THz regime, for example, σ is dominated by the intraband transitions. For higher optical frequencies such as mid- and near-infrared, the interband contribution is important. Although in this model we aim to demonstrate a THz metamaterial, for completeness we still include both the intraband and interband contributions in the conductivity calculation.

Next, we model the graphene-based THz metamaterial absorber, which is two graphene layers patterned into a periodic fishnet structure (Figure 1) embedded in a polymer material TOPAS. The bottom plane is a metal ground plane. The refractive index of TOPAS is 1.53 in THz domain. The geometric parameters for the patterned graphene can be found in Ref. 1. The common approach to model a periodic structure like this is to use periodic boundary condition. However, in this particular case, the unit cell has mirror symmetry. Consequently, in cases where we only consider normal incidence, Perfect Electric Conductor (PEC) and Perfect Magnetic Conductor (PMC) conditions can be used as symmetry planes such that only a quarter of the unit cell needs to be modeled, which greatly reduces the simulation time. In addition, we will not model the thickness of graphene explicitly but use the TBC. This further improves the efficiency of the simulation.

Results and Discussion

The calculated graphene sheet conductivity is shown in Figure 2. It shows the real and imaginary parts of σ for a few different Fermi energies. The intraband transition leads to a Drude-like response similar to a typical metal.

Figure 2: Calculated sheet conductivity of graphene at different Fermi energies. The solid curves are the real parts and the dashed curves are the corresponding (negative) imaginary parts.

The simulated absorption spectra at different Fermi energies are shown in Figure 3. Even though the polymer matrix has no absorption, the Fabry–Perot resonator formed by the graphene and the ground plane leads to perfect absorption around 3 THz at 0.5 eV Fermi energy.

Figure 3: Absorption spectra of the graphene-based metamaterial at different Fermi energies.

Notes About the COMSOL Implementation

COMSOL provides two convenient features that make calculating the optical conductivity of graphene much easier. First of all, for complicated equations like Equation 1 and Equation 2, usually it is very important to make sure all the quantities are used with the correct and consistent unit system. In COMSOL this is not required since the unit conversion will be done automatically. Secondly, COMSOL provides a list of built-in constants. We do not need to look up the values of physical constants such as the reduced Planck constant, the Boltzmann constant, the speed of light, and the electron charge. They can be directly referred to as hbar_const, k_B_const, c_const, and e_const.

In the current model, we use the TBC to model graphene. Alternatively, graphene can be modeled as 2D surface current or a 3D layer with an effective thickness, as discussed in the blog posts: www.comsol.com/blogs/modeling-graphene-in-high-frequency-electromagnetics and www.comsol.com/blogs/should-we-model-graphene-as-a-2d-sheet-or-thin-3d-volume.

Reference

1. A. Andryieuski and A.V. Lavrinenko, “Graphene metamaterials based tunable terahertz absorber: effective surface conductivity approach,” Opt. Express, vol. 21, pp. 9144–9155, 2013.

Wave_Optics_Module/Gratings_and_Metamaterials/graphene_metamaterial_perfect_absorber

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
T	300[K]	300 K	Temperature
Ef	0.2[eV]	3.2044E-20 J	Fermi energy
tau	1e-13[s]	1E-13 s	Relexation time
d_eff	1[nm]	1E-9 m	Effective thickness of graphene
a	15[um]	1.5E-5 m	Unit cell length
b	2[um]	2E-6 m	Geometric parameter 1
d_sub	17.6[um]	1.76E-5 m	Substrate thickness
d_domain	40[um]	4E-5 m	Simulation domain height
w	12[um]	1.2E-5 m	Geometric parameter 2
n_bg	1.53	1.53	Substrate refractive index

Definitions

In the window, expand the > node.

Right-click and choose > .

In the window for , type H in the text field.

In the text field, type H.

Locate the section. In the text field, type sinh(hbar_const*x/(k_B_const*T))/(cosh(hbar_const*x/(k_B_const*T))+cosh(Ef/(k_B_const*T))).

Locate the section. In the table, enter the following settings:


Argument	Unit
x	rad/s

Locate the section. In the table, enter the following settings:


Plot	Argument	Lower limit	Upper limit	Fixed value	Unit
√	x	0	1e16	0	rad/s

Variables 1

Right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
Omega	1[rad/s]	rad/s	Variable of integral
integral	integrate((H(Omega)-H(ewfd.omega/2))/(ewfd.omega^2-4*Omega^2),Omega,0[rad/s],1e16[rad/s])	s	Integral in the interband conductivity equation
sigma_intra	((2k_B_constTe_const^2)/(pihbar_const^2))(log(2cosh(Ef/(2k_B_constT)))*(-j/(ewfd.omega-j/tau)))	S	Intraband conductivity
sigma_inter	(e_const^2/(4hbar_const))(H(ewfd.omega/2)-(j4ewfd.omega/pi)*integral)	S	Interband conductivity
sigma	sigma_intra+sigma_inter	S	Total graphene conductivity

Here, the integrate operator is used to perform the numerical integration.

Geometry 1

In the window, expand the > node, then click .

In the window for , locate the section.

From the list, choose .

Block 1 (blk1)

In the toolbar, click

In the window for , locate the section.

In the text field, type a/2.

In the text field, type d_domain.

Work Plane 1 (wp1)

In the toolbar, click

In the window for , locate the section.

In the text field, type d_sub.

Work Plane 1 (wp1) > Plane Geometry

In the window, click .

Work Plane 1 (wp1) > Square 1 (sq1)

In the toolbar, click

In the window for , locate the section.

In the text field, type w/2.

Work Plane 1 (wp1) > Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type b/2.

In the text field, type (a-w)/2.

Locate the section. In the text field, type w/2.

Work Plane 1 (wp1) > Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type (a-w)/2.

In the text field, type b/2.

Locate the section. In the text field, type w/2.

In the toolbar, click

Due to the mirror symmetry of the unit cell, it is sufficient to model a quarter of the structure, which greatly improves the simulation efficiency.

Materials

Background

In the window, under right-click and choose .

In the window for , type Background in the text field.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Refractive index, real part	n_iso ; nii = n_iso, nij = 0	n_bg	1	Refractive index

Graphene

Right-click and choose .

In the window for , type Graphene in the text field.

Click to expand the section. In the tree, select > .

Click

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Electric conductivity	sigma_iso ; sigmaii = sigma_iso, sigmaij = 0	sigma/d_eff	S/m	Basic

The 3D conductivity is the 2D sheet conductivity divided by the effective thickness.

Locate the section. From the list, choose .

Click the

button in the toolbar.

Select Boundaries 6, 8, and 10 only.

Electromagnetic Waves, Frequency Domain (ewfd)

By default, a PEC condition is applied to all exterior boundaries. We consider a plane wave polarized in the x direction that is incident on the unit cell in the normal direction. Therefore, PEC conditions apply on the surfaces perpendicular to the polarization direction, while PMC conditions are used on the surfaces parallel to the polarization direction.

Perfect Magnetic Conductor 1

In the toolbar, click