unit test and tutorial for LDOS in planar cavity using cylindrical coordinates

[oskooi]

Based on a suggestion in #2076 (comment).

In cylindrical coordinates, there is a noticeable difference in the Purcell enhancement factor compared to the analytic result at the Van Hove singularity at a cavity thickness of 1.5 which does not seem to go away with resolution. Perhaps this could be due to the E_r polarized point-dipole source at r=0?

[stevengj]

Second source is redundant.

[stevengj]

Mismatch right at the van-Hove singularity is probably due to the fact that the cylindrical PML in Meep is currently a "quasi"-PML, and the errors in the PML have a more severe effect as the group velocity goes to zero, which happens right at the singularity. Solution: either increase the PML thickness, or increase the non-PML thickness, or both.

[stevengj]

(Note that what you have here is a circularly polarized source, while the theory is for linearly polarized, but that's okay. Linearly polarized is just the average of left- and right-circularly polarized, and the powers in the two cases just add, so it should be the same result.)

[oskooi]

Solution: either increase the PML thickness, or increase the non-PML thickness, or both.

At a resolution of 200 for a cavity thickness of 1.5 (the Van Hove singulariy), I tried increasing both the PML and non-PML thicknesses but the LDOS value via the PE factor is unchanged around 1.1. For comparison, the theoretical result is 1.5555. That's a fairly large error of ~30%. This could suggest there is still a bug in the Er source at r=0 but then again the PE factor at all other points (i.e., which are not Van Hove singularities) do agree with the theory.

[stevengj]

I can't think of any explanation other than convergence problems that would cause it to be wrong at one frequency but not others. Maybe try going to much thicker, e.g. by 5x or 10x.

(Probably the critical parameter is the thickness of the quasi-PML.)

You could also double-check that a circularly polarized source in 3d doesn't change the result.

[oskooi]

I think I have figured out what the issue here is which is not related to the quasi PML in cylindrical coordinates. Rather, similar to what was described in #1277, small changes in the cavity thickness (z dimension of the cell) due to discretization errors introduced by the choice of the grid resolution can produce large changes in the LDOS which is exacerbated because it is a Van Hove singularity. This is demonstrated in the plot below showing the measured Purcell enhancement factor at the Van Hove singularity (at a cavity thickness of 0.5 in units of wavelength in the cavity medium) as a function of the grid resolution. The PE factor varies by a factor of more than 2 when changing the resolution by just one unit.

In my previous results which were computed at a resolution of 200 pixels/μm, I just happened to be unlucky and chose a resolution in which the discretization errors dominated the results. We can point this out in the tutorial.

import meep as mp
import numpy as np


n = 2.4            # refractive index of surrounding medium                                                                                                  
wvl = 1.0          # wavelength (in vacuum)                                                                                                                  

fcen = 1/wvl
sources = [mp.Source(src=mp.GaussianSource(fcen,fwidth=0.2*fcen),
                     component=mp.Er,
                     center=mp.Vector3())]


def bulk_ldos(res):
    L = 5.0        # length of non-PML region in r/z                                                                                                         
    dpml = 1.0     # thickness of PML                                                                                                                        
    sr = L+dpml
    sz = L+2*dpml
    cell_size = mp.Vector3(sr,0,sz)

    pml_layers = [mp.PML(dpml)]

    sim = mp.Simulation(resolution=res,
                        cell_size=cell_size,
                        boundary_layers=pml_layers,
                        sources=sources,
                        dimensions=mp.CYLINDRICAL,
                        m=-1,
                        default_material=mp.Medium(index=n))

    sim.run(mp.dft_ldos(fcen,0,1),
            until_after_sources=mp.stop_when_fields_decayed(20,
                                                            mp.Er,
                                                            mp.Vector3(),
                                                            1e-9))

    return sim.ldos_data[0]


def cavity_ldos(res,sz):
    Lr = 25.0        # length of non-PML region in r                                                                                                         
    dpml = 5.0       # thickness of PML in r                                                                                                                 
    sr = Lr+dpml
    cell_size = mp.Vector3(sr,0,sz)

    pml_layers = [mp.PML(dpml,direction=mp.R)]

    sim = mp.Simulation(resolution=res,
                        cell_size=cell_size,
                        boundary_layers=pml_layers,
                        sources=sources,
                        dimensions=mp.CYLINDRICAL,
                        m=-1,
                        default_material=mp.Medium(index=n))

    sim.run(mp.dft_ldos(fcen,0,1),
            until_after_sources=mp.stop_when_fields_decayed(20,
                                                            mp.Er,
                                                            mp.Vector3(),
                                                            1e-7))

    return sim.ldos_data[0]

def purcell_enh_theory(ct):
    return (3*np.fix(ct+0.5)/(4*ct) +
            (4*np.power(np.fix(ct+0.5),3) -
             np.fix(ct+0.5)) /
            (16*np.power(ct,3)))


if __name__ == '__main__':
    # Van Hove singularity                                                                                                                                   
    # units of wavelength in cavity medium                                                                                                                   
    cavity_thickness = 1.50

    gap = cavity_thickness*wvl/n
    pe_theory = purcell_enh_theory(cavity_thickness)
    print("enh-theory:, {:.6f} (gap), {:.6f} (theory)".format(gap,pe_theory))

    for res in range(30,101):
        ldos_bulk = bulk_ldos(res)
        ldos_cavity = cavity_ldos(res,gap)
        pe_meep = ldos_cavity/ldos_bulk
        print("enh:, {}, {:.6f}".format(res,pe_meep))

[codecov-commenter]

Codecov Report

Merging #2082 (85121e9) into master (1265469) will increase coverage by 0.06%.
The diff coverage is n/a.

@@            Coverage Diff             @@
##           master    #2082      +/-   ##
==========================================
+ Coverage   73.11%   73.18%   +0.06%     
==========================================
  Files          17       17              
  Lines        4895     4941      +46     
==========================================
+ Hits         3579     3616      +37     
- Misses       1316     1325       +9     

Impacted Files	Coverage Δ
python/adjoint/optimization_problem.py	54.67% <0.00%> (-1.80%)	⬇️
python/simulation.py	76.93% <0.00%> (+<0.01%)	⬆️
python/adjoint/objective.py	92.96% <0.00%> (+0.89%)	⬆️

[stevengj]

Exactly at the frequency of the van Hove singularity, the LDOS is not really well defined, because there is a discontinuity at that point. If you look at frequencies slightly below or above this point you should get more sensible convergence results with resolution.