2d Cell with Out-of-Plane Wavevector


A 2d cell with a k_point that has a non-zero component (e.g., planewaves incident on a 2d structure from an out-of-plane direction, propagating modes of fiber waveguides with 2d cladding cross section, etc.) could be accomplished by a "1-pixel-thick" 3d simulation with complex fields (this is done by the kz_2d="3d" option). However, Meep will by default model the dependence using a modified 2d simulation with complex fields (this is done by the kz_2d="complex" option), which improves performance with no loss in accuracy (as demonstrated below). As a further optimization, Meep can also model this problem with "real" fields via kz_2d="real/imag", but the fields must be interpreted in a special way.

Mathematically, an dependence () of the fields can be treated by including an cross-product in the curls of Maxwell's equations, which couples the and polarizations. Since this term is complex (), the electromagnetic fields are complex.

However, an additional trick is possible. Since the complex term only couples and polarizations to one another, then we can choose the polarization ( fields) to be purely real while the polarization ( fields) is purely imaginary. The kz_2d="real/imag" option does precisely this, but stores the purely imaginary polarization as the "real" parts of the fields, both internally and in the output. So, if you use the kz_2d="real/imag" and you output both and , the output will be real numbers, but you should multiply by to get the actual fields. This requires some care, which is why this option is not the default. The good news is that calculations of flux, energy, forces and similar real quantities are insensitive to this implicit factor (the cancels), so the built-in calculations do the right thing with no modification.

As a demonstration of this feature, consider the example of computing the reflectance of a planar air/dielectric interface for a planewave incident from the out-of-plane direction. This is a slightly modified version of Tutorials/Basics/Angular Reflectance Spectrum of a Planar Interface. The interface normal is along the axis. An -polarized planewave is incident at an angle of 19.4° counter clockwise around the -axis with 0° along the direction (i.e., the plane of incidence is which corresponds to the -polarization). The 2d cell is in the plane (with a dimension of a single pixel). The reflectance at a wavelength of 1 μm is computed using separate simulations for kz_2d of "real/imag", "complex", and "3d". The three results are compared with the analytic Fresnel equation.

The simulation script is examples/refl-angular-kz2d.py.

import meep as mp
import math

def refl_planar(theta, kz_2d):
    resolution = 100

    dpml = 1.0
    sx = 10
    sx = 10 + 2*dpml
    cell_size = mp.Vector3(sx)
    pml_layers = [mp.PML(dpml)]

    fcen = 1.0

    # plane of incidence is XZ
    k = mp.Vector3(z=math.sin(theta_r)).scale(fcen)

    sources = [mp.Source(mp.GaussianSource(fcen,fwidth=0.2*fcen),
                         component=mp.Ey,
                         center=mp.Vector3(-0.5*sx+dpml))]

    sim = mp.Simulation(cell_size=cell_size,
                        boundary_layers=pml_layers,
                        sources=sources,
                        k_point=k,
                        kz_2d=kz_2d,
                        resolution=resolution)

    refl_fr = mp.FluxRegion(center=mp.Vector3(-0.25*sx))
    refl = sim.add_flux(fcen, 0, 1, refl_fr)

    sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ey, mp.Vector3(-0.5*sx+dpml), 1e-9))

    empty_flux = mp.get_fluxes(refl)
    empty_data = sim.get_flux_data(refl)
    sim.reset_meep()

    # add a block with n=3.5 for the air-dielectric interface
    geometry = [mp.Block(mp.Vector3(0.5*sx,mp.inf,mp.inf),
                         center=mp.Vector3(0.25*sx),
                         material=mp.Medium(index=3.5))]

    sim = mp.Simulation(cell_size=cell_size,
                        geometry=geometry,
                        boundary_layers=pml_layers,
                        sources=sources,
                        k_point=k,
                        kz_2d=kz_2d,
                        resolution=resolution)

    refl = sim.add_flux(fcen, 0, 1, refl_fr)
    sim.load_minus_flux_data(refl, empty_data)

    sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ey, mp.Vector3(-0.5*sx+dpml), 1e-9))

    refl_flux = mp.get_fluxes(refl)
    freqs = mp.get_flux_freqs(refl)

    Rmeep = -refl_flux[0]/empty_flux[0]
    return Rmeep


# rotation angle of source: CCW around Y axis, 0 degrees along +X axis
theta_r = math.radians(19.4)

Rmeep_real_imag = refl_planar(theta_r,"real/imag")
Rmeep_complex = refl_planar(theta_r,"complex")
Rmeep_3d = refl_planar(theta_r,"3d")

n1=1
n2=3.5

# compute angle of refracted planewave in medium n2
# for incident planewave in medium n1 at angle theta_in
theta_out = lambda theta_in: math.asin(n1*math.sin(theta_in)/n2)

# compute Fresnel reflectance for S-polarization in medium n2
# for incident planewave in medium n1 at angle theta_in
Rfresnel = lambda theta_in: math.fabs((n2*math.cos(theta_out(theta_in))-n1*math.cos(theta_in))/(n2*math.cos(theta_out(theta_in))+n1*math.cos(theta_in)))**2

print("refl:, {:.16f} (real/imag), {:.16f} (complex), {:.16f} (3d), {:.16f} (analytic)".format(Rmeep_real_imag,Rmeep_complex,Rmeep_3d,Rfresnel(theta_r)))

The Meep results are identical to within six decimal digits and agree well with the analytic theory with an error of less than 0.7%.

refl:, 0.3272338236967464 (real/imag), 0.3272338244372344 (complex), 0.3272330216564413 (3d), 0.3293821216165117 (analytic)