In this example, we investigate the effect of photothermal heating in plasmonic nanostructures based on the published work in the referenced article [1]. Using Lumerical's DGTD solver for optical simulation and HEAT solver for thermal simulation, we look at two types of thermoplasmonic antennae; dipole and diabolo. In this example we consider arrays of both these antennae under varying optical intensity and compare their performance.
Requirements
Lumerical products R2018a or newer
Simulation Setup
|
Optical simulation setup (partitioned volume mode) |
Optical (DGTD)
Based on the experimental setup of Ref. [1], in this part of the example we consider two arrays of antennae, one using diabolo antennae and the other one using dipole antennae, both made out of gold. Open the diabolo_array_full.ldev project file. The array of diabolo antennae is sitting on top of a thick Sapphire (Al2O3) substrate. The dimensions of the antennae are taken from Ref. [1]. The array has a pitch of 340 nm in both X and Y axes [1]. The DGTD simulation region covers one unit cell centered on one diabolo antenna and has a width equal the period of the array in both X and Y direction. The length (along Z) of the simulation region is kept at 0.4 micron (in addition to another 0.6 micron for the two "shell" regions for the top and bottom PML layers) . A plane wave source is used to illuminate the antennae array at a wavelength of 1064 nm.
An absorption monitor is used to measure the absorbed optical power in the antenna. The monitor automatically normalizes the absorbed power with respect to the source power (CW normalization) so the 'power density' and the 'total power' is reported for an input optical power of 1 W.
Thermal (HEAT)
The thermal simulation is also performed for one unit cell centered on one diabolo antenna and so the x and y span of the HEAT simulation region is identical to the DGTD one. The main difference is that in the thermal simulation we include the entire thickness of the substrate in order to model the propagation of the heat through the substrate. The thickness of the substrate is set to 250 micron which is a typical value for commercially available sapphire wafers. An Import heat source is used to import the optical absorption data saved from the DGTD simulation. Two temperature monitors are placed surrounding the antenna to record the temperature profile. A sweep (Pin) is set up in the project file to apply a scaling factor on the heat source to perform multiple simulations under varying input powers.
The fixed temperature thermal boundary condition is applied at the bottom of the simulation region to set it at room temperature (300 K). Two convection boundary conditions are used to model the heat loss due to convection from gold to air and from sapphire to air. The model for convection is chosen to be constant with a convective heat transfer coefficient h = 10 W/m2K.
The dipole_array_full.ldev project file uses a similar setup using dipole antennae where each antenna is 215 nm long, 50 nm wide, and 50 nm thick [1].
Results and Discussion
Optical (DGTD)
Optical Resonance
Open the diabolo_array_full.ldev project file. By default the HEAT solver is disabled so the DGTD solver is the only active solver. Open the property editor for the plane wave source (source) and change the frequency range to start from 0.9 micron and stop at 1.2 micron. Open the property editor of the frequency domain EM field monitor (T) and enable the "use source limits" option. Also set the number of "frequency points" to 40. Do the same for the absorption monitor (Pabs). Run the simulation and visualize the "flux" result from the "T" monitor. View the transmission (T_front) data. The plot will show that the transmission has a minima around 1064 nm which is consistent with the findings of Ref. [1]. Alternatively, right click on the "Pabs" monitor and plot "total" to visualize the total absorbed power which also shows a maxima around 1064 nm indicating that the dipole antenna has a resonance at a wavelength of 1064 nm. The same simulation can be done for the dipole antenna using the dipole_array_full.ldev project file.
Transmission past the diabolo antenna as a function of wavelength
Power absorbed in the diabolo antenna as a function of wavelength
Optical absorption
Next, edit the property of the two monitors again (in the diabolo_array_full.ldev project file) back to the original setup (disable the "use source limits" option, set both the start and stop wavelengths at 1.064 micron to record the data at a single frequency only, and set the number of "frequency points" to 1). Run the simulation again and once the simulation is run, right click on the "Pabs" monitor and visualize the absorbed power density profile (density). Note that the plot below shows the absorbed power in log scale (can be set from chart settings) to clearly show the spatial distribution of the absorbed power. The max and min limit of the color bar has also been modified.
Absorbed power at a certain depth (Z value) inside the diabolo antenna
Use the following script commands to read the data from the monitor and save it in a .mat file to be used in the HEAT simulation.
density = getresult("DGTD::Pabs","density");
matlabsave("pabs_diabolo_array.mat",density);
Following the same approach, simulate and save the absorbed power for the dipole antenna array using the dipole_array_full.ldev project file.
Thermal (HEAT)
Open the diabolo_array_full.ldev project file. Enable the HEAT solver object and disable the DGTD solver. The heat source (heat) already has the optical absorption data loaded but you can also load the pabs_diabolo_array.mat file onto the heat source and select the Pabs attribute to import the data saved by FDTD simulation. In the DGTD simulation the absorbed power was calculated for an input power of 1 mW applied over a single antenna. From the experimental setup of Ref. [1], the incident optical beam had a FWHM of 19.5 micron and covered 2583 antennae. This means that for a total input optical power of 1 mW the power applied to a single antenna (for the unit cell) is 1/2583 mW or 3.87147e-07 W. Therefore we should use a scaling factor of 3.87147e-07 in the import heat source to run the thermal simulation for the case where the total input optical power is 1 mW.
Set the scaling factor in the heat source (in Data tab of the object properties window) to 20*3.87147e-07 (or 7.74294e-06) corresponding to an optical input of 20 mW and run the simulation. Once the simulation is run, the results will be loaded in the HEAT solver region and in the temperature monitors. Right click on the first monitor (monitor_1) and visualize temperature to see the temperature profile.
Temperature profile at the antenna and the top surface of the substrate
The diabolo_array_full.ldev project also has a parameter sweep object (Pin) created in the optimizations and sweep window to vary the input optical power from 5 mW to 125 mW in 7 steps and run the simulations (note that the value of the scale factor is P x 3.87147e-07 where P is the input power in mW as discussed earlier). The sweep object reads the temperature profile from the 2D temperature monitor (monitor_2) and saves them in a dataset (T). Since the finite element mesh of the HEAT solver is random in nature, saving the temperature distribution from all seven simulations (at different power) on the same dataset requires locking the mesh so that the mesh is identical for all the simulation files created by the parameter sweep. For this reason the mesh of the HEAT solver needs to be calculated and locked using the option in the "Simulation" menu option before running the sweep. Use the script file plot_Pin_sweep_result.lsf to lock the mesh, run the sweep, and calculate the temperature of the entire chip as a function of input optical power using the sweep result.
In order to be able to calculate the temperature of the entire chip we assume that the 2583 antennae in the center of the array that gets illuminated by the incident light are at an uniform temperature given by the peak temperature in monitor_2 and the surrounding region where the antennae are not illuminated is at a constant temperature equal to the room temperature (300 K). the temperature distribution is thus assumed to have the top-hat profile shown below (left). The total surface area is considered to be 85 micron x 85 micron as per Ref. [1]. By comparing the ratio of the illuminated and unilluminated area, the increment in average surface temperature is found to be equal to 0.0413*(Tsim - 300). The script calculates this rise in average temperature and plots it as a function of input optical power. This averaging is done so that the simulation result for a single unit cell can be compared against the measured temperature of the entire experimental setup. The figure below (right) shows the temperature rise for both the diabolo and dipole arrays as a function of input power. The plot is in agreement with Fig. 4(b) of Ref. [1].
Temperature profile for the entire array of antennae
Temperature rise versus input optical power for both antenna arrays
Related publications
[1] Z. J. Coppens, W. Li, D. G. Walker, and J. G. Valentine, "Probing and Controlling Photothermal Heat Generation in Plasmonic Nanostructures," Nano Letters, vol. 13, pp. 1023-28, 2013.