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).
If you are new to the finite-difference time-domain (FDTD) method, we highly recommend going through our FDTD101 tutorials. For more simulation examples, please visit our examples page. FDTD simulations can diverge due to various reasons. If you run into any simulation divergence issues, please follow the steps outlined in our troubleshooting guide to resolve it.
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.
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.
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.
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$.
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.
# 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.
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.
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.
batch = web.Batch(simulations=sims, folder_name="default")
batch_results = batch.run(path_dir="data")
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- e7f9c895-55f6-411b-9941-8f6e30852b7av1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- f2d0defb-9bf8-4085-9a71-23505ffd1ea7v1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 095e3ae5-0053-4177-829c-e820b2a53350v1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- d6b08b2f-efef-4cb0-864d-19c302858c67v1'.
[11:29:02] Started working on Batch. container.py:475
[11:29:30] Maximum FlexCredit cost: 0.100 for the whole batch. container.py:479 Use 'Batch.real_cost()' to get the billed FlexCredit cost after the Batch has completed.
[11:29:37] Batch complete. container.py:522
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.
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.
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.
# 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.
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.
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")
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 61224bd7-a026-4ce7-80e2-66f5162ad1f2v1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 100e69b5-312e-424f-95df-fb32aa20d691v1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 40cd4b72-e67b-421e-aeec-88aa4d845b57v1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 6e529e36-0a18-41b9-a2ce-bc13a60623aev1'.
[11:29:45] Started working on Batch. container.py:475
[11:30:06] Maximum FlexCredit cost: 0.100 for the whole batch. container.py:479 Use 'Batch.real_cost()' to get the billed FlexCredit cost after the Batch has completed.
[11:30:18] Batch complete. container.py:522
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.
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).
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.
graphene.medium.plot_sigma(freqs)
plt.show()
