Question about non-uniform permittivities

[pwflanigan]

I'm having some trouble with the example (https://meep.readthedocs.io/en/latest/Python_Tutorials/Mode_Decomposition/#diffraction-spectrum-of-liquid-crystal-polarization-gratings and https://github.com/NanoComp/meep/blob/master/python/examples/polarization_grating.py) used to show how to implement non-homogeneous (position-dependent) permittivities.

1. I tried to visualize this custom-made material (defined with the function lc_mat) using both plot2d and get_epsilon, but in both cases the material was completely uniform (aside from the air interface were subpixel averaging was happening). Is there a way to visualize how epsilon_diag and epsilon_offdiag are implemented?
2. I'm assuming that "p" in the code is rho = (x,y,0); is that defined relative to the center of the simulation space or the center of the geometric object?
3. Could you please explain a little bit in general about how plot2d deals with non-scalar permittivities?

Thanks for the help.

[oskooi]

Any offdiagonal elements of the epsilon tensor are currently ignored by plot2D.

The epsilon data used by plot2D to generate the image is obtained using get_epsilon_grid:

meep/python/visualization.py

Lines 600 to 602 in bba8d5b

	eps_data = np.rot90(
	np.real(sim.get_epsilon_grid(xtics, ytics, ztics, eps_parameters["frequency"]))
	)

As described in the documentation, get_epsilon_grid "computes the trace of the ε(f) tensor at frequency f."

[pwflanigan]

Thanks. Could you please explain a little more about the trace thing? I was playing with some examples but I can't seem to understand how to get it to work -- for example, shouldn't these two custom-made materials have the same color in plot2D, or am I missing something?

import meep as mp
import matplotlib.pyplot as plt
import numpy as np

outer_material = mp.Medium(epsilon_diag=mp.Vector3(1,2,9))
inner_material = mp.Medium(epsilon_diag=mp.Vector3(4,4,4)) # mp.Medium(epsilon=4)

trace_outer = np.trace(outer_material.epsilon(freq=mp.inf))
trace_inner = np.trace(inner_material.epsilon(freq=mp.inf))

geometry = []
geometry.append(mp.Cylinder(radius=2, material=outer_material))
geometry.append(mp.Cylinder(radius=1, material=inner_material))

sim = mp.Simulation(resolution=50,
                    cell_size=mp.Vector3(4,4,0),
                    geometry=geometry,
                    sources=[])

plt.figure(dpi=150)
sim.plot2D()

Also, are there other examples besides the nematic grating of how to implement the position-dependent permittivity?

[stevengj]

What Meep actually stores is $\varepsilon^{-1}$, not $\varepsilon$, and if I remember correctly, what it plots is $\frac{3}{\text{trace}[\varepsilon^{-1}]}$, i.e. the harmonic mean of the eigenvalues of $\varepsilon$.

Computed here:

meep/src/array_slice.cpp

Lines 385 to 408 in bba8d5b

	else if (cS[i] == Dielectric) {
	complex<double> tr(0.0, 0.0);
	for (int k = 0; k < data->ninveps; ++k) {
	tr += (fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx, frequency) +
	fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx + ieos[2 * k], frequency) +
	fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx + ieos[1 + 2 * k], frequency) +
	fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx + ieos[2 * k] + ieos[1 + 2 * k],
	frequency));
	if (abs(tr) == 0.0) tr += 4.0; // default inveps == 1
	}
	fields[i] = IVEC_LOOP_WEIGHT(s0i, s1i, e0i, e1i, 1.0) * (4.0 * data->ninveps) / tr;
	}
	else if (cS[i] == Permeability) {
	complex<double> tr(0.0, 0.0);
	for (int k = 0; k < data->ninvmu; ++k) {
	tr += (fc->s->get_chi1inv_at_pt(imcs[k], imds[k], idx, frequency) +
	fc->s->get_chi1inv_at_pt(imcs[k], imds[k], idx + imos[2 * k], frequency) +
	fc->s->get_chi1inv_at_pt(imcs[k], imds[k], idx + imos[1 + 2 * k], frequency) +
	fc->s->get_chi1inv_at_pt(imcs[k], imds[k], idx + imos[2 * k] + imos[1 + 2 * k],
	frequency));
	if (abs(tr) == 0.0) tr += 4.0; // default invmu == 1
	}
	fields[i] = IVEC_LOOP_WEIGHT(s0i, s1i, e0i, e1i, 1.0) * (4.0 * data->ninvmu) / tr;
	}

[pwflanigan]

OK, that's good to know. But I'm still confused about the non-homogeneous permittivity -- in the nematic grating example, I took out parts not related to that issue (to save time) and included a plot2d line, but I don't see any position-dependence in the grating region. Am I not understanding this example? Is there another example that shows how to implement the position-dependence in all entries of the permittivity matrix?

import math
import matplotlib.pyplot as plt
import numpy as np
import meep as mp

resolution = 50  # pixels/μm
dpml = 1.0  # PML thickness
dsub = 1.0  # substrate thickness
dpad = 1.0  # padding thickness
k_point = mp.Vector3(0, 0, 0)
pml_layers = [mp.PML(thickness=dpml, direction=mp.X)]
n_0 = 1.55
delta_n = 0.159
epsilon_diag = mp.Matrix(
    mp.Vector3(n_0**2, 0, 0),
    mp.Vector3(0, n_0**2, 0),
    mp.Vector3(0, 0, (n_0 + delta_n) ** 2))
wvl = 0.54  # center wavelength
fcen = 1 / wvl  # center frequency

def pol_grating(d, ph, gp, nmode):
    sx = dpml + dsub + d + d + dpad + dpml
    sy = gp

    cell_size = mp.Vector3(sx, sy, 0)

    def phi(p):
        xx = p.x - (-0.5 * sx + dpml + dsub)
        if (xx >= 0) and (xx <= d):
            return math.pi * p.y / gp + ph * xx / d
        else:
            return math.pi * p.y / gp - ph * xx / d + 2 * ph

    def lc_mat(p):
        Rx = mp.Matrix(
            mp.Vector3(1, 0, 0),
            mp.Vector3(0, math.cos(phi(p)), math.sin(phi(p))),
            mp.Vector3(0, -math.sin(phi(p)), math.cos(phi(p))))
        lc_epsilon = Rx * epsilon_diag * Rx.transpose()
        lc_epsilon_diag = mp.Vector3(lc_epsilon[0].x, lc_epsilon[1].y, lc_epsilon[2].z)
        lc_epsilon_offdiag = mp.Vector3(lc_epsilon[1].x, lc_epsilon[2].x, lc_epsilon[2].y)
        return mp.Medium(epsilon_diag=lc_epsilon_diag, epsilon_offdiag=lc_epsilon_offdiag)

    geometry = [
        mp.Block(
            center=mp.Vector3(-0.5 * sx + 0.5 * (dpml + dsub)),
            size=mp.Vector3(dpml + dsub, mp.inf, mp.inf),
            material=mp.Medium(index=n_0)),
        mp.Block(
            center=mp.Vector3(-0.5 * sx + dpml + dsub + d),
            size=mp.Vector3(2 * d, mp.inf, mp.inf),
            material=lc_mat)]

    src_pt = mp.Vector3(-0.5 * sx + dpml + 0.3 * dsub, 0, 0)
    sources = [
        mp.Source(
            mp.GaussianSource(fcen, fwidth=0.05 * fcen),
            component=mp.Ez,
            center=src_pt,
            size=mp.Vector3(0, sy, 0)),
        mp.Source(
            mp.GaussianSource(fcen, fwidth=0.05 * fcen),
            component=mp.Ey,
            center=src_pt,
            size=mp.Vector3(0, sy, 0))]

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

    plt.figure(dpi=150)
    sim.plot2D()
    
    return 0

ph_uniaxial = 0  # chiral layer twist angle for uniaxial grating
ph_twisted = 70  # chiral layer twist angle for bilayer grating
gp = 6.5  # grating period
nmode = 5  # number of mode coefficients to compute
dd = 1.00

pol_grating(0.5 * dd, math.radians(ph_uniaxial), gp, nmode)

[pwflanigan]

After looking into this a little more, it's possible that the reason plot2D isn't showing a clear position-dependence in epsilon is because it just so happens that because of the way epsilon is plotted (the 3/Tr[1/epsilon] calculation) that any variation gets obscured, even though the tensors really are changing from point-to-point. I'm not 100% sure if that's it, but does that seem plausible?

Here's a more direct question: I wanted to confirm this using sim.get_epsilon(), but I got a AttributeError: 'NoneType' object has no attribute 'total_volume' message, does anyone know how to get around that? Specifying a volume with vol= didn't work.