implementing graphene, epsilon infinity value of 5.5

[mdfeinstein]

Hello Nikolaos, I have looked at your code for graphene in a coordinate stretch and tried to implement the same in my simulations. When working with a coordinate stretch only along the wave propagation direction (normal to the graphene surface), I got what seemed like reasonable results (comparing to simulations done in Lumerical as well as experimental results), but it was a little hard to tell without increasing the resolution in the graphene plane, which led me to attempting a coordinate stretch in the graphene plane as well. When doing this, though, I found the results were highly dependent on the stretch factors in the plane. When taking graphene out of the simulation (leaving gold rods on a substrate only), the results were stable with respect to the stretch factors, so I assume it has something to do with the graphene coordinate transformation.

My main question is, how did you land on an epsilon infinity value of 5.5? it seems from me from the theory it should depend on the background relative permittivity, though I am not completely sure what to do when there are different background materials on each side (and potentially spatially dependent as well, or dispersive materials).

Also, I have included below a function I wrote for creating graphene material object. The user specifies the Fermi Energy and scattering rate, as well as optionally a set of compression factors, a relative epsilon for the graphene, and which direction is the normal. I used the constants in your code to convert from eV to meep frequency units, and to calculate the Omega_p by scaling your value of 0.64 eV by the desired eV value. However, as stated, I feel that this method is not stable with respect to the stretch factors. If you spot anything wrong, please do let me know!

def graphene_maker_3dstretch(ef, tau, resolution, compression=[1,1,1], epsilon=1, normal=mp.X):
    '''
    returns a mp.Medium with a Drude susceptibility approximating graphene's intraband response.
    ef is fermi energy in eV
    tau is scatering rate in eV
    resolution is meep resolution
    compression is 3-iterable of x,y,z compression factors. 
    '''
    ev_to_f = lambda ev : ev / 1.23984193 
    Omega_p= np.sqrt(ef/0.64)*3.225094 # eV, value for 0.34nm thickness (https://iopscience.iop.org/article/10.1088/1742-6596/129/1/012004/meta),  Graphene Fermi level=0.64eV 
    Gama = tau # 
    epsilon_inf= 5.5 # epsilon infinity
    Sigma= None
    #define sigma tensor by setting normal component to 0
    sigma_diag = [1,1,1]
    sigma_diag[normal]=0
    sigma_diag = mp.Vector3(*sigma_diag) 
    # eV (normalize to meep thickness 1/(resolution*stretch) taking into account resolution and stretch
    Omega_p=Omega_p*np.sqrt(0.34*10**-9)/np.sqrt(1/resolution*10**-6)/np.sqrt(1/compression[normal]) 
    omega = ev_to_f(Omega_p)
    gamma = ev_to_f(tau)
    drude = mp.DrudeSusceptibility(frequency=omega, gamma=gamma, sigma=Sigma, sigma_diag = sigma_diag) 
    trf = mp.Matrix(diag=mp.Vector3(*compression))
    graphene = mp.Medium( E_susceptibilities = [drude], epsilon_diag=mp.Vector3(epsilon,epsilon, epsilon) )
    graphene.transform(trf)
    return graphene

[Nikolaos-Matthaiakakis]

Hello, it has been a while since I worked on this so I don't remember all the details, and as I am very busy with other projects and research proposals I don't currently have the time to investigate sadly.

This method can be tricky, so instead I decided to brute force the graphene simulations for my work with a 2nm resolution and periodic boundaries instead of PML orthogonal to the meta-atom. This way you can avoid large simulation areas by avoiding the space requirements for the PML. Typically experiments use graphene meta-atom arrays instead of isolated graphene meta-atoms so either way this method is usually more suitable.

But if you want to continue working on this method, maybe you should consider checking how changing the stretching factors could potentially affect your effective thickness of graphene and if the conductivity should be adjusted accordingly. I am saying that without having checked your code though, maybe things are already fine and there is another issue that causes the problem.

Regarding the eps_inf value please use the relevant value for your case. If you want to compare with the lumerical model you should probably use the same exact values used there.