Negative frequencies in spectra calculation - problem with materials or sources?

[Srinivas-Ganti]

Hello all,

My question today is related to material definitions in the THz regime, following my previous question here. At THz frequencies, the Drude model is sufficient to describe metals. I used Drude parameters from this source and followed the example given in the meep materials library for Al that uses Drude parameters referenced from a techinal note from Horiba

My first question - Is there anything wrong with my definitions? I am not sure as I am getting unexpected output in the simple bent waveguide example I tried next following the tutorial given in the docs.

MATERIAL DEFINITIONS

import meep as mp

um_scale = 1
eV_Hz_scale = 2.417990504024e14
eV_THz_scale = eV_Hz_scale*1e-12
eV_um_scale = um_scale * 1.23984193
thz_um_scale = 299.792458*um_scale # 1THz = 299.792458 um

wp = 1.29e16 # plasma frequency in Hz
wt = 4.12e13 # damping frequency in Hz
Au_eps_inf = 1 # instantaneous dielectric response 
wp_eV = wp/eV_Hz_scale # plasma frequency in eV
wt_eV = wt/eV_Hz_scale # damping frequency in eV

Au_Drude_frq = 1/(eV_um_scale/wp_eV)  # frequencies in meep units (inverse microns)
Au_Drude_gam = 1/(eV_um_scale/wt_eV)

# frequencies in THz are converted to um
end_wl = thz_um_scale/0.001 # 0.0001 THz
start_wl = thz_um_scale/20 # 20 THz

Au_Drude_sig = Au_eps_inf
Au_Drude_range = mp.FreqRange(min= 1/end_wl, max= 1/start_wl)
Au_Drude_susc = [
    mp.DrudeSusceptibility(
        frequency=Au_Drude_frq, gamma=Au_Drude_gam, sigma=Au_Drude_sig
    )
]

Au_Drude_THzMaterial = mp.Medium(
    epsilon=Au_eps_inf, E_susceptibilities=Au_Drude_susc, valid_freq_range=Au_Drude_range
)
JGS1_valid_range = Au_Drude_range
JGS1_THzMaterial = mp.Medium(epsilon=1.954**2,valid_freq_range=JGS1_valid_range)

[Srinivas-Ganti]

The second question is about the spectra: I am using just the JGS1_THzMaterial definition to recreate the bent waveguide example but am getting negative wavelengths and obviously something is not right in my code perhaps but I'm not able to find my mistake... Please could anyone offer a suggestion if they know?

During the simulation I get two RuntimeWarnings repeatedly, which of course are not seen when running the tutorial code as is. When I animate the fields through the waveguide it looks fine, but not the spectra. I attach the simulation script, Ez fields animation and sim cell setup I used.

RuntimeWarning: DFT frequency -0.005000000000000001 is out of material's range of 3.3356409519815205e-06-0.0667128190396304
  warnings.warn(warn_dft_fmt.format(dftf, min_freq, max_freq), RuntimeWarning)

RuntimeWarning: Note: your sources include frequencies outside the range of validity of the material models. This is fine as long as you eventually only look at outputs (fluxes, resonant modes, etc.) at valid frequencies.
  warnings.warn(warn_src, RuntimeWarning)

SIMULATION SCRIPT

import meep as mp
from MyMaterials.AuTHzDrude import *

um_scale = 1.0 

resolution = 2 # px/ um

cell_l = cell_w = 300  # cellsize in x and y directions
cell= mp.Vector3(cell_w,cell_l,0)

dpml = 10
pml_layers = [mp.PML(dpml)]

wg_w = 12.8 # width of waveguide
wg_l = 240  # length of waveguide

pad = 46 # padding distance between waveguide and cell edge

src_xpos = -cell_l/2 + dpml          # src location
src_ypos = cell_l/2 - dpml - pad

# a straight waveguide 

geometry = [mp.Block(mp.Vector3(mp.inf , wg_w, mp.inf),
            center = mp.Vector3(-cell_l/2 + dpml + wg_l/2, cell_l/2 - dpml - pad),
            material = JGS1_THzMaterial)]
            

fcen = 1/50    # lam = 50 um
df = 0.05      # <= pulse width in frequency. <=> 20 um     
sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df),
                        component=mp.Ez,
                        center=mp.Vector3(src_xpos,src_ypos),
                        size=mp.Vector3(0,wg_w))]



sim = mp.Simulation(cell_size= cell, 
                    boundary_layers = pml_layers,
                    geometry = geometry,
                    sources = sources,
                    resolution = resolution)

nfreq = 100  # number of frequencies at which to compute flux

# reflected flux
refl_fr = mp.FluxRegion(center = mp.Vector3(-0.5*cell_l+dpml+1, cell_w/2-dpml-pad),
                        size = mp.Vector3(0,2*wg_w,0))
refl = sim.add_flux(fcen, df, nfreq, refl_fr)

# transmitted flux
tran_fr = mp.FluxRegion(center = mp.Vector3(0.5*cell_l-dpml-1, cell_w/2-dpml-pad),
                        size = mp.Vector3(0,2*wg_w,0))
tran = sim.add_flux(fcen, df, nfreq, tran_fr)

# point to monitor field decay
pt = mp.Vector3(0.5*cell_l-dpml-2,cell_w/2 - dpml - pad)

sim.run(until_after_sources=mp.stop_when_fields_decayed(500,mp.Ez,pt,1e-5))

# for normalization run, save flux fields data for reflection plane

straight_refl_data = sim.get_flux_data(refl)

# save incident power for transmission plane
straight_tran_flux = mp.get_fluxes(tran)

sim.reset_meep()
            
geometry = [mp.Block(mp.Vector3(wg_l , wg_w, mp.inf),
            center = mp.Vector3(-cell_l/2 + dpml + wg_l/2, cell_l/2 - dpml - pad),
            material = JGS1_THzMaterial),
            
            mp.Block(mp.Vector3(wg_w, wg_l, mp.inf),
            center = mp.Vector3(cell_l/2 - dpml - pad ,-cell_l/2 + dpml + wg_l/2),
            material = JGS1_THzMaterial)]

sim = mp.Simulation(cell_size=cell,
                    boundary_layers=pml_layers,
                    geometry=geometry,
                    sources=sources,
                    resolution=resolution)

sim.use_output_directory()



# reflected flux
refl = sim.add_flux(fcen, df, nfreq, refl_fr)
tran_fr = mp.FluxRegion(center=mp.Vector3(cell_l/2 - dpml - pad,-0.5*cell_w+dpml+1,0),
                        size=mp.Vector3(2*wg_w,0,0))
tran = sim.add_flux(fcen, df, nfreq, tran_fr)

# for normal run, load negated fields to subtract incident from refl. fields
sim.load_minus_flux_data(refl, straight_refl_data)

pt = mp.Vector3(0.5*cell_l-dpml -pad, -0.5*cell_w+dpml+2)


# if mp.am_really_master():
#     #visualise cell with bent waveguide, src and monitors
#     plt.figure()
#     sim.plot2D(output_plane=mp.Volume(size = mp.Vector3(cell_l,cell_w,0)),
#                eps_parameters={'resolution' : 10}, plot_monitors_flag=True,
#                labels = True)
#     plt.show()


sim.run(mp.at_beginning(mp.output_epsilon),
        mp.to_appended("ez", mp.at_every(2, mp.output_efield_z)),
        until_after_sources=mp.stop_when_fields_decayed(500, mp.Ez, pt, 1e-5))

bend_refl_flux = mp.get_fluxes(refl)
bend_tran_flux = mp.get_fluxes(tran)

flux_freqs = mp.get_flux_freqs(refl)

wl = []
Rs = []
Ts = []
for i in range(nfreq):
    wl = np.append(wl, 1/flux_freqs[i])
    Rs = np.append(Rs,-bend_refl_flux[i]/straight_tran_flux[i])
    Ts = np.append(Ts,bend_tran_flux[i]/straight_tran_flux[i])
print("<<<<<<<<<<<<<<<<<<< WAVELENGTHS")
print(wl)
if mp.am_really_master():
    plt.figure()
    plt.plot(wl,Rs,'bo-',label='reflectance')
    plt.plot(wl,Ts,'ro-',label='transmittance')
    plt.plot(wl,1-Rs-Ts,'go-',label='loss')
    plt.axis([wl[-1], wl[0], 0, 1])
    plt.xlabel("wavelength (μm)")
    plt.legend(loc="upper right")
    plt.show()

[stevengj]

Every real signal is a superposition of positive and negative frequencies, so this is correct. You really need to understand this basic fact of Fourier theory if you want to look at spectra.

This looks a bit odd to me:

fcen = 1/50    # lam = 50 um
df = 0.05 

Here, your fcen is 0.02, but your bandwidth is 0.05, so your bandwidth extends over a pretty wide range (from -0.005 to +0.045). Is that really what you want or are you confused about units here?

(By default, Meep uses purely real fields, just as in actual physical systems. However, you can force it to use complex-valued fields by setting force_complex_fields=True in the Simulation constructor. In this case, you can get signals with only positive (or only negative) frequencies.)

[Srinivas-Ganti]

Thank you for your answer and explanation about the negative frequencies. I see I went overboard in my setup ^^.

I am only interested in the positive frequencies as I am trying to replicate my THz TDS experiment. In the experiment, we have a pretty wide band pulse (about 4-7 THz depending on how long we average). I will try the sim again with a more reasonable bandwidth to avoid the negative frequencies. It looked odd to me that in the result posted we see negative values for the plotted transmittance. I don't understand what that means physically here.

The real THz pulse from the spectrometer is a complicated-looking waveform, so the Gaussian pulse is only a starting point for me, although the Custom Source feature of meep looks very interesting. I am going to try and fit the experimental data and use it as input for Meep's custom source after fitting.