In this example, 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. The unperturbed metamaterial has been analyzed optically in the FDTD metatmaterial example. Here, the structure is analyzed electrically using the CHARGE first and optically using FDTD later.
NOTE: If you are using an older version of CHARGE (v.4.6.631 or older) then please download the simulation (.ldev) file from this archived KB page as the current files may not run with the older CHARGE.
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.
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 npdensity grid attribute in FDTD.
1. When sweeping the simulation positions of x=0 um, 3.5 um, 9.5 um, 16 um and 21.5 um, resulting data from the charge monitor will be recorded in different output .mat files, which is set in the script in "model".
2. Since no other output is needed, in the "Results" of the sweep, we set the current to output, which is arbitrary. If you do not set any result, after sweep, you will need to manually close the sweep dialogue.
3. The sweep of applied voltages is set in the Boundary Conditions table for the electrical contact named "emitter".
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 feild etc at different bias for different position.
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 npdensity 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");
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.