Density of states calculation

[pavlo2014]

Dear All,

I am puzzled with the output of dft_ldos. Please find below the band diagrams of a Bragg stack calculated using MPB and MEEP, DOS spectrum and the python code. I am using one of example codes supplied with MEEP here to demonstrate the problem.

The band diagram as calculated using MPB

As calculated with MEEP

The spectrum of density of states (DOS) calculated with MEEP

ABS(DOS) if you want to see it

The bands calculated with MPB and MEEP agree well. However, the DOS spectrum indicates states in the band gap from 0.22 to 0.43. There is also a peak in DOS at f=0.4. I would not expect to see any states in the band gap.

Am I calculating/plotting the DOS wrong or its is time to review my understanding the theory of DOS?

Please feel free to point me towards a suitable source of information to save your time.

best wishes,

pavlo

`import meep as mp
from meep import mpb
import matplotlib.pyplot as plt
import numpy as np

n_lo = 1.0
n_hi = 3.0

w_hi = n_lo / (n_hi + n_lo) # a quarter_wave stack

geometry_lattice = mp.Lattice(size=mp.Vector3(1)) # 1d cell
default_material = mp.Medium(index=n_lo)
geometry = mp.Cylinder(material=mp.Medium(index=n_hi), center=mp.Vector3(), axis=mp.Vector3(1),
radius=mp.inf, height=w_hi)

kx = 0.5
k_points = mp.interpolate(50, [mp.Vector3(), mp.Vector3(kx)])

resolution = 64
num_bands = 3

def main():

# 1. Band diagrams using MPB 

ms = mpb.ModeSolver(num_bands=num_bands,
                    k_points=k_points,
                    geometry_lattice=geometry_lattice,
                    geometry=[geometry],
                    resolution=resolution,
                    default_material=default_material)

ms.run_te()
freqs = ms.all_freqs.tolist()

fig = plt.figure(figsize=(20,15))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])

for row in freqs:
    for elem in row:
        plt.scatter(k_points[freqs.index(row)][0], elem.real, color='black')
    
plt.plot(k_points, k_points, color='red')

ax.set_xlabel('k-point')
ax.set_ylabel("f (c/a)")
ax.set_xlim([0, 0.5])
ax.set_ylim([0, 0.8])    

plt.grid()
plt.savefig('bandsMPB.png')
fig.show()

# 2. Band diagrams using MEEP

cell = mp.Vector3(1, 1)

s = mp.Source(src=mp.GaussianSource(0.3, fwidth=1.0), component=mp.Hz,
              center=mp.Vector3(0))
              
sim = mp.Simulation(cell_size=cell, geometry=[geometry], sources=[s], resolution=resolution)

freqsMEEP = sim.run_k_points(300, k_points)

fig = plt.figure(figsize=(20,15))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])

for row in freqsMEEP:
    for elem in row:
        plt.scatter(k_points[freqsMEEP.index(row)][0], elem.real, color='black')
    
plt.plot(k_points, k_points, color='red')

ax.set_xlabel('k-point')
ax.set_ylabel("f (c/a)")
ax.set_xlim([0, 0.5])
ax.set_ylim([0, 0.8])

plt.grid()
plt.savefig('bandsMEEP.png')
fig.show()

# 3. DOS diagrams using MEEP

sim.reset_meep()

sim.run(mp.dft_ldos(0.0, 0.8, 1000), until_after_sources=10000)


ldos = sim.ldos_data
f = np.linspace(0.0, 0.8, len(ldos))

fig = plt.figure(figsize=(20,15))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])

plt.plot(f, ldos)

ax.set_xlabel('f (c/a)')
ax.set_ylabel("DOS")

plt.grid()
plt.savefig('ldos.png')
fig.show()   

if name == 'main':
main()`

[stevengj]

cell = mp.Vector3(1, 1)

This looks wrong. You have a 2d unit cell in the Meep calculation, which will give you additional eigenvalues due to band-folding effects.

I once asked about something like this as an exam question, actually see problem 2 from my 2007 midterm and solutions thereof.

[pavlo2014]

Dear Steven,

Many thanks for your response and links to the exam questions and answers. They are useful, it is always good to learn something new. This is an interesting effect and it deserves attention.

Yes, defining the the two-dimensional cell was wrong. Somehow I merged your code and my old code for a 2D simulation and ended up with this error.

The cell is now defined as "cell = mp.Vector3(1)".

The band diagram looks correct.

But the LDOS spectrum still have the peak at f=0.4.

Is the cause of this peak known? There are no eigenvalues at these frequencies of the band digram and I would not expect any states to be at this frequency too.

best wishes,
pavlo

[stevengj]

(Note also that dft_ldos computes the local density of states (LDOS), not the DOS.)

[pavlo2014]

Thank you, noted.

[stevengj]

Something weird is going on here because the LDOS should always be ≥ 0.