Learning FDTD through Examples / Graphene metamaterial absorber...

Graphene metamaterial absorber¶

Graphene, a single atomic sheet of carbon atoms arranged in a hexagonal lattice, has shown great promise for functional optical and optoelectronic devices. Compared to devices made of conventional materials, graphene-based devices have a unique tunability advantage since graphene's conductivity can be drastically modulated by electrostatic gating.

Due to its atomic thickness, it is usually very difficult to model graphene in an FDTD simulation without using a large number of grid points. Fortunately, Tidy3D natively supports a surface conductivity model such that thin material layers can be accurately simulated even with grids much larger than the actual layer thickness. More specifically, you can model graphene directly using Tidy3D's Material Library. The graphene's conductivity is described by the well established Kubo formula, which includes the contributions from the intraband and interband electronic transitions. To define graphene, a few parameters, namely the scattering rate, chemical potential, and temperature, are required as user inputs.

This notebook demonstrates how to model a graphene fishnet metamaterial absorber in the THz frequency range. By forming a Fabry-Perot resonator between the graphene layer and a metal mirror, high absorption is achieved across a bandwidth of ~ 2 THz. The design is adapted from the seminal work Andrei Andryieuski and Andrei V. Lavrinenko, "Graphene metamaterials based tunable terahertz absorber: effective surface conductivity approach," Opt. Express 21, 9144-9155 (2013).

[1]:

import numpy as np
import matplotlib.pyplot as plt

import tidy3d as td
import tidy3d.web as web

The overall design of the absorber involves two layers of graphene sheets embedded in TOPAS polymer. The bottom of the device is a metallic mirror, which forms a Fabry-Perot resonator with graphene. This resonator greatly amplifies the absorption in graphene and achieves perfect absorption when the resonance condition is met. Two designs are explored in this notebook. The first one is a simple uniform graphene layer. When the chemical potential of graphene is tuned to 0.5 eV, perfect absorption is achieved but only within a narrow frequency band. Then, we test a fishnet structure, which shows a much wider absorption band.

Uniform Graphene Sheet Absorber¶

We start with a simpler case where the absorber is made of two uniform graphene layers embedded in TOPAS polymer above a ground plate.

First, define the basic simulation parameters.

[2]:

THz = 1e12  # 1 THz = 1e12 Hz

freq0 = 2.5 * THz  # central frequency
freqs = np.linspace(0.1, 5, 300) * THz  # frequency range of interest

lda0 = td.C_0 / freq0  # central wavelength

The TOPAS polymer has a refractive index of 1.53 in the THz range with little absorption.

[3]:

n_topas = 1.53  # refractive index of the polymer
topas = td.Medium(permittivity=n_topas**2)

To define graphene conductivity, we need to specify a few parameters. Here, we consider a relaxation time of $\tau$=0.1 ps, which translates to a scattering rate $\Gamma=\frac{1}{2\tau}=\frac{6.58e-16}{2 \times 1e-13}=$0.0033 eV, where 6.58e-16 is $\hbar$ in the unit of eV*s. The temperature is at room temperature. Since there are two layers of graphene, the scaling parameter is set to 2. Note that this only works for graphene layers separated by a small distance. For real bilayer graphene or few-layer graphene in general, the conductivity is more complex as interlayer coupling and stacking order play an important role.

[4]:

gamma = 0.0033  # scattering rate
temp = 300  # temperature
scaling = 2  # number of layers.

Furthermore, we consider a unit cell size of 15 $\mu m$. In the uniform graphene sheet absorber, the distance between the graphene and the mirror is 35.9 $\mu m$.

[5]:

a = 15  # unit cell size
h = 35.9  # distance between graphene and the ground plate
offset = lda0 / 2  # distance between the flux monitor and the graphene

The periodic boundary condition is applied in both the $x$ and $y$ directions. In the $z$ direction, PML is applied in the plus direction and PEC is applied in the minus direction. The PEC mimics the metal mirror. For the grids, we use automatic nonuniform grids in all directions. The same BoundarySpec and GridSpec will be used through the notebook.

[6]:

# define a BoundarySpec
boundary_spec = td.BoundarySpec(
    x=td.Boundary.periodic(),
    y=td.Boundary.periodic(),
    z=td.Boundary(minus=td.PECBoundary(), plus=td.PML()),
)

# define a GridSpec
grid_spec = td.GridSpec.auto(min_steps_per_wvl=200, wavelength=lda0)

# simulation run time
run_time = 1e-11

Next we define a function that returns a simulation at a given graphene chemical potential. This function will later be called repeatedly to construct a simulation batch to investigate how the chemical potential affects the absorption spectrum.

[7]:

def make_sim_uniform(mu_c):

    # define graphene
    graphene = td.material_library["graphene"](
        gamma=gamma, mu_c=mu_c, temp=temp, scaling=scaling
    ).medium

    # define graphene structure
    graphene_layer = td.Structure(
        geometry=td.Box(center=(0, 0, h), size=(td.inf, td.inf, 0)), medium=graphene
    )

    # define a plane wave source
    plane_wave = td.PlaneWave(
        center=(0, 0, h + 0.1 * offset),
        size=(td.inf, td.inf, 0),
        source_time=td.GaussianPulse(freq0=freq0, fwidth=freq0 / 2),
        direction="-",
    )

    # define a flux monitor to measure reflection
    flux_monitor = td.FluxMonitor(
        center=(0, 0, h + offset),
        size=(td.inf, td.inf, 0),
        freqs=freqs,
        name="R",
    )

    # simulation domain size in z
    Lz = h + 1.1 * offset

    # define simulation
    sim = td.Simulation(
        center=(0, 0, Lz / 2),
        size=(a, a, Lz),
        grid_spec=grid_spec,
        structures=[graphene_layer],
        sources=[plane_wave],
        monitors=[flux_monitor],
        run_time=run_time,
        medium=topas,
        boundary_spec=boundary_spec,
        shutoff=1e-7,
        symmetry=(-1, 1, 0),  # symmetry is used to reduce the computational load
    )

    return sim

Specifically, we study $\mu_c=$0, 0.1, 0.2, and 0.5 eV.

[8]:

mu_cs = [0, 0.1, 0.2, 0.5]  # values of mu_c to be simulated

# define a simulation batch
sims = {f"mu_c={mu_c:.2f}": make_sim_uniform(mu_c) for mu_c in mu_cs}

Submit the simulation batch to the server.

[9]:

batch = web.Batch(simulations=sims, folder_name="default")
batch_results = batch.run(path_dir="data")

After the batch of jobs is complete, we want to plot the absorption spectra. Since the bottom ground plane is a perfect mirror, absorption is simply $A=1-R$, where $R$ is the reflection. A function plot_absorption is defined to do this plotting.

[10]:

def plot_absorption(batch_results):
    for i, mu_c in enumerate(mu_cs):
        sim_data = batch_results[f"mu_c={mu_c:.2f}"]
        A = 1 - sim_data["R"].flux
        plt.plot(freqs / THz, A, label=f"{mu_c:.2f} eV")
        plt.xlim(0, 5)
        plt.ylim(0, 1)
        plt.xlabel("Frequency (THz)")
        plt.ylabel("Absorption")
        plt.legend()

At 0.5 eV chemical potential, perfect absorption is observed around 2.3 THz. However, the bandwidth is rather small.

[11]:

plot_absorption(batch_results)

Graphene Fishnet Absorber¶

To enhance the absorption bandwidth, a fishnet design is explored, as discussed in the reference.

Similarly, we define a function here that returns a simulation at a given graphene chemical potential.

[12]:

# fishnet parameters
b = 2
w = 12

h = 17.6  # distance between the graphene and the ground plate
inf_eff = 1e3  # effective infinity


def make_sim_fishnet(mu_c):
    # define graphene
    graphene = td.material_library["graphene"](
        gamma=gamma, mu_c=mu_c, temp=temp, scaling=scaling
    ).medium

    # define the fishnet structure
    graphene_finishnet_center = td.Structure(
        geometry=td.Box.from_bounds(rmin=(0, 0, h), rmax=(w / 2, w / 2, h)),
        medium=graphene,
    )
    graphene_finishnet_right = td.Structure(
        geometry=td.Box.from_bounds(rmin=(w / 2, 0, h), rmax=(inf_eff, b / 2, h)),
        medium=graphene,
    )
    graphene_finishnet_top = td.Structure(
        geometry=td.Box.from_bounds(rmin=(0, w / 2, h), rmax=(b / 2, inf_eff, h)),
        medium=graphene,
    )

    # define a plane wave source
    plane_wave = td.PlaneWave(
        center=(0, 0, h + 0.1 * offset),
        size=(td.inf, td.inf, 0),
        source_time=td.GaussianPulse(freq0=freq0, fwidth=freq0 / 2),
        direction="-",
    )

    # define a flux monitor to measure reflection
    flux_monitor = td.FluxMonitor(
        center=(0, 0, h + offset),
        size=(td.inf, td.inf, 0),
        freqs=freqs,
        name="R",
    )

    # simulation domain size in z
    Lz = h + 1.1 * offset

    # define simulation
    sim = td.Simulation(
        center=(0, 0, Lz / 2),
        size=(a, a, Lz),
        grid_spec=grid_spec,
        structures=[
            graphene_finishnet_center,
            graphene_finishnet_right,
            graphene_finishnet_top,
        ],
        sources=[plane_wave],
        monitors=[flux_monitor],
        run_time=run_time,
        medium=topas,
        boundary_spec=boundary_spec,
        shutoff=1e-7,
        symmetry=(-1, 1, 0),
    )

    return sim

To make sure the structures are defined correctly, use the plot method to visualize the fishnet. Due to symmetry, only a quarter of it needs to be defined.

[13]:

sim = make_sim_fishnet(0.5)
ax = sim.plot(z=h)

A batch of simulations at different chemical potentials is defined and submitted to the server.

[14]:

mu_cs = [0, 0.1, 0.2, 0.5]  # values of mu_c to be simulated

# define a simulation batch
sims = {f"mu_c={mu_c:.2f}": make_sim_fishnet(mu_c) for mu_c in mu_cs}

# submit the batch to the server
batch = web.Batch(simulations=sims, folder_name="default")
batch_results = batch.run(path_dir="data")

At 0.5 eV, the absorber shows an absorption bandwidth (defined as absorption above 0.9) of ~2 THz, much better than the uniform graphene absorber case.

[15]:

plot_absorption(batch_results)

Additional Notes about Graphene¶

As briefly discussed in the introduction, graphene's conductivity includes the contributions from both the intraband and interband electronic transitions. The intraband transitions lead to a Drude-like response that is dominant at lower frequencies. The interband transitions give rise to the prominent feature around $f=2\mu_c$.

The conductivity of graphene can be inspected by using the numerical_conductivity method. Here we show an example of graphene conductivity at $\Gamma$=0.0033 eV, $\mu_c$=0.2 eV, and $T$=300 K. Below 20 THz, the conductivity is predominantly Drude-like, which comes from the intraband transitions. The interband transition feature shows up at around 100 THz (~0.4 eV).

[16]:

gamma = 0.0033  # scattering rate
mu_c = 0.2  # chemical potential
temp = 300  # temperature

freqs = np.linspace(10, 250, 100) * THz
graphene = td.material_library["graphene"](gamma=gamma, mu_c=mu_c, temp=temp)
sigma_analytical = graphene.numerical_conductivity(freqs)
plt.plot(freqs / THz, np.real(sigma_analytical * 1e6), label="Re($\sigma$) ($\mu$S)")
plt.plot(freqs / THz, np.imag(sigma_analytical * 1e6), label="Im($\sigma$) ($\mu$S)")
plt.xlabel("frequency (THz)")
plt.title("analytically calculated surface conductivity")
plt.legend()
plt.show()

Like other dispersive materials used in FDTD, graphene's conductivity needs to be fitted, which is automatically done when you define a graphene medium. The plot generated by numerical_conductivity above is the analytical result. To inspect the fitting, one can use the plot_sigma method as shown below. The fitted conductivity coincides well with the analytical result in this case.

[17]:

graphene.medium.plot_sigma(freqs)
plt.show()

[ ]: