Active THz metamaterial

In this example, we'll show how to recreate the results of H.-T. Chen et al. This metamaterial exhibits a negative index of refraction in the THz range. First, we will optically analyze the unperturbed passive metamaterial using FDTD. Then, we will study the effect of bias-induced carrier density variations on the refractive index of the individual materials and hence the overall transmission of the metamaterial. In this step, the structure is analyzed electrically using CHARGE first and then optically using FDTD.

Passive Device

Simulation setup

The file negative_index_chen.fsp contains an example of the structure as described in Figure 1 from Chen's paper.  

The metamaterial structure here is composed of gold patterning on a GaAs substrate. Gold can be represented using a plasma material model (with plasma resonance (ω) and collision frequency (vc)), and this can be expressed as a simple conductive model in the low frequency limit when ω << vc. In this limit, the material can be approximated using the PEC (perfect electrical conductor) material.

The physical thickness of the gold is 200 nm, but this would require a very small mesh size for dz. However, since the thickness of the material is much less than the wavelength of about 130-1200 um, the metal can be represented using 2D sheets which do not require a fine mesh in the z-direction.

For the GaAs substrate, a simple constant refractive index (dielectric) model was used.  It may be possible to account for a free carrier model by adding a conductivity that depends on bias voltage.

Results (using dx = 1 nm, dy = 1 nm)

In the file provided, the mesh override region uses a coarser mesh which gives initial results in a shorter simulation time. More accurate results can be obtained by specifying a finer mesh step size. The following results are achieved by changing the mesh step size in the x and y directions in the mesh override region to 1 nm.

When the simulation is finished, run negative_index_chen_analysis.lsf to produce the results.

The transmission as seen in Figure 3 from Chen.  The peak and dip just above 1.4 THz in the spectrum is a Wood's anomaly due which occurs as the number of supported grating orders changes at this frequency.

The |E| field is similar to Figure 2a.  Rerunning the simulation with a smaller mesh size would produce a higher resolution image.

The vector plot of the surface current density at resonance. Similar to figure 2b of the paper. Once you have run the script to generate the vector plot, modify the plot properties as follows:

Active Device

Theory

A range of voltages are applied to the metamaterial via contacts and the corresponding carrier densities are calculated and recorded. The effect of the carrier densities on the refractive index of the GaAs layer is  then calculated.

To calculate the effect of the change in the carrier density on the refractive index, an FDTD simulation will be run. The np density grid attribute in FDTD will take the carrier density information and calculate the corresponding changes in the real and imaginary parts of refractive index of the material according to the Plasma-Drude formulation. For a more detailed description of  this grid attribute and the formula, please visit the Charge to index conversion page.

CHARGE Simulation

The metamaterial including the GaAs epilayer and the gold stripes can be set up in CHARGE. Open the active_thzmaterial.ldev file in CHARGE. The dimensions of the structure reflect the structures in the referenced paper by Chen et al.

In this example, The CHARGE solver uses a 2D simulation domain, several x-normal cross sections of the structure are simulated individually for a full picture of the carrier density profile across the entire metamaterial. A parameter sweep task has been set up under the sweeps and optimizations tab to sweep over the position of the simulation region for 5 different cross sections. The electron density is then collected for each of the cross sections, simulating voltages between 0 and 16 volts for every one of the cross sections. A charge monitor will record the n and p distributions and save them to .mat files for the corresponding cross sections, which will be imported into the np density grid attribute in FDTD.  

If the sweep task in CHARGE solver is run, the location of the simulation region object is varied from the center to the edge of the structure sweeping over 33 voltage points at each cross section. Once the sweep is run, we can open the different sweep files in the newly created folder and look at results such as the charge distribution, electrostatic field etc at different bias for different position.

FDTD Simulation

In FDTD, to symmetry boundaries are used and so only the first quarter of the simulation region is included, as shown in THz device. We only create the np density objects in that quarter as well. Carrier densities are used to calculate the corresponding refractive index of GaAs according to the Plasma-Drude model explained of Charge to index conversion. Note that the x span has to be specified according to the cross section of the device. This has already been done in the provided .fsp file so you can skip this.

To calculate the overall device transmission as a function of the applied voltage, a parameter sweep is used in FDTD under the sweeps and optimizations task. The file is set to only sweep voltage of 0 V, 8 V, and 16 V when the emitter voltage has been changed to positive in "model". Users with multi-processor computers or access to concurrent computing resources (e.g. a compute cluster) can take advantage of these extra processors by configuring the resources in FDTD to utilize all available cores during the parameter sweep.

Run the sweep task in the .fsp file. The sweep will run simulations for all three bias values at a range of frequencies from 0.25 THz to 2.5 THz. Use the following lines of script to plot the results once the sweep is done.  The plot is similar with the curve in Figure 3a of the referenced paper.

f = pinch(getsweepdata("sweep","f"));
T = pinch(getsweepdata("sweep","T"));
V = pinch(getsweepdata("sweep","V"));
plot(f(1:100,1)*1e-12, T(1:100,1),T(1:100,2),T(1:100,3), "Frequency (THz)", "Transmission");
legend("V=0 volts", "V=8 volts", "V=16 volts");

Discussion

Compared to the reference paper, one may find there are some differences from the above result. For example  at higher frequency, the transmission is higher than the paper. This is because, we approximate the structure as symmetric not only in x direction, but also in y direction. If removing the symmetry in y, the transmission at higher frequency will be much lower. However, since the paper did not give complete information on the refractive index of GaAs, the doping concentration and doping profile/region, and we simplified the simulation of n-GaAS and SI-GaAs as GaAS with a constant refractive index of 3.5 in FDTD and set electric properties in CHARGE, we simplified the structure as cross sections, and doping causes loss even without bias voltage, the differences are reasonable.

Related publications

H.-T. Chen et al., "Active terahertz metamaterial devices", Nature 444, 597-600 (2006)

See also

Bulk metamaterials

Metamaterial S parameter extraction  

Negative index metamaterial in GHz range using wire pairs