stevengj
commented
6 years ago Given the complex mode coefficients of the linearly-polarized modes, there is a simple formula to compute the coefficients of the circularly polarized modes.
Closed oskooi closed 5 years ago
stevengj
commented
6 years ago Given the complex mode coefficients of the linearly-polarized modes, there is a simple formula to compute the coefficients of the circularly polarized modes.
stevengj
commented
5 years ago Simple approach with existing code:
Do two simulations: (a) one for TM (z-polarized) and (b) one for TE (xy-polarized) input. Let's say the light is incident from the x direction, and we have mirror symmetry in y. The TM input is ODD_Z+EVEN_Y, and the TE input is EVEN_Z+ODD_Y.
In each simulation, compute the complex mode amplitudes for both TM and TE outputs (i.e. all of the mode parities). If you have mirror symmetry in y, then the y parity is preserved, but the z parity is not because z symmetry is broken by the anisotropy of the medium in this example. So, for the TM input you will calculate both ODD_Z+EVEN_Y and EVEN_Z+EVEN_Y coefficients. For the TE input you will compute both ODD_Z+ODD_Y and EVEN_Z+ODD_Y coefficients.
Circularly polarized light is TM ± iTE. Since this is a linear problem, the same relationship applies to the outputs: the mode amplitudes you would get from circularly polarized light should be (TM amplitudes) ± i(TE amplitudes). There is a subtlety, however, because the TM amplitudes are EVEN_Y and the TE amplitudes are ODD_Y: you should think of each TM amplitude as (+y)+(-y) and each TE amplitude as (+y)-(-y) if you decompose them into beams diffracting in the ±y directions. So, if you have (TM) + i(TE) circular input, then to get the corresponding +y output amplitude you should compute (TM amplitudes) + i(TE amplitudes), whereas to get the -y diffracted output you compute (TM amplitudes) - i(TE amplitudes). Vice versa for (TM) - i(TE) circular input.
(This approach actually gives you a factor of two speed advantage compared to simulations with circularly polarized inputs, because it allows you to exploit the y mirror symmetry.)
oskooi
commented
5 years ago The structure in this example breaks the mirror symmetry in y because the nematic director of the liquid crystal is rotating around x. (Fig. 1a of the reference shows the nematic director at just a single cross section where the twist angle happens to be 0° and therefore appears to be mirror symmetric in x. However, for any non-zero twist angle this symmetry is broken.) Exploiting the symmetry to obtain a factor of two speed advantage is therefore negated.
The following script computes the first 5 TM (ODD_Z+EVEN_Y) and TE (EVEN_Z+ODD_Y) mode coefficients for a TM (ODD_Z+EVEN_Y) input. As usual, a normalization run with an empty cell is used to obtain the input flux of the source. At the end of the simulation, the complex mode coefficients as well as the transmittance is displayed for each mode. For reference, the reflectance and transmittance computed using the Poynting flux are also shown.
# polarization grating from C. Oh and M.J. Escuti, Optics Letters, Vol. 33, No. 20, pp. 2287-9, 2008
# note: reference uses z as the propagation direction and y as the out-of-plane direction; this script uses x and z, respectively
import meep as mp
import math
resolution = 50 # pixels/um
wvl = 0.54 # center wavelength
fcen = 1/wvl # center frequency
df = 0.05*fcen # frequency width
dpml = 1.0 # PML thickness
dsub = 1.0 # substrate thickness
dpad = 1.0 # padding thickness
d = 1.7 # chiral layer thickness
ph = math.radians(70) # twist angle of each chiral layer
gp = 6.5 # grating period
sx = dpml+dsub+d+d+dpad+dpml
sy = gp
cell_size = mp.Vector3(sx,sy,0)
pml_layers = [mp.PML(thickness=dpml,direction=mp.X)]
# twist angle of nematic director; from equation 1b
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
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))
# return the anisotropic permittivity tensor for a uniaxial, twisted nematic liquid crystal
def lc_mat(p):
# rotation matrix for rotation around x axis
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)]
# linearly polarized planewave: TM (mp.Ez) or TE (mp.Hz)
src_cmpt = mp.Ez
src_pt = mp.Vector3(-0.5*sx+dpml+0.3*dsub,0,0)
sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df), component=mp.Ez, center=src_pt, size=mp.Vector3(0,sy,0), amplitude=1 if src_cmpt==mp.Ez else 0),
mp.Source(mp.GaussianSource(fcen,fwidth=df), component=mp.Hz, center=src_pt, size=mp.Vector3(0,sy,0), amplitude=1 if src_cmpt==mp.Hz else 0)]
k_point = mp.Vector3(0,0,0)
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
boundary_layers=pml_layers,
k_point=k_point,
sources=sources,
default_material=mp.Medium(index=n_0))
refl_pt = mp.Vector3(-0.5*sx+dpml+0.5*dsub,0,0)
refl_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,sy,0)))
sim.run(until_after_sources=100)
input_flux = mp.get_fluxes(refl_flux)
input_flux_data = sim.get_flux_data(refl_flux)
sim.reset_meep()
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
boundary_layers=pml_layers,
k_point=k_point,
sources=sources,
geometry=geometry)
refl_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,sy,0)))
sim.load_minus_flux_data(refl_flux,input_flux_data)
tran_pt = mp.Vector3(0.5*sx-dpml-0.5*dpad,0,0)
tran_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=tran_pt, size=mp.Vector3(0,sy,0)))
sim.run(until_after_sources=300)
nmode = 5 # number of mode coefficients to compute
Tsum = 0
res1 = sim.get_eigenmode_coefficients(tran_flux, range(1,nmode+1), eig_parity=mp.ODD_Z+mp.EVEN_Y)
coeffs1 = res1.alpha
res2 = sim.get_eigenmode_coefficients(tran_flux, range(1,nmode+1), eig_parity=mp.EVEN_Z+mp.ODD_Y)
coeffs2 = res2.alpha
for nm in range(nmode):
kdom = res1.kdom[nm]
angle = math.degrees(math.acos(kdom.x/fcen))
tran1 = abs(coeffs1[nm,0,0])**2/input_flux[0]
tran2 = abs(coeffs2[nm,0,0])**2/input_flux[0]
print("mode1:, {}, {:.2f}, {:.5f}{:+.5f}j, {:.8f}".format(nm,angle,coeffs1[nm,0,0].real,coeffs1[nm,0,0].imag,tran1))
print("mode2:, {}, {:.2f}, {:.5f}{:+.5f}j, {:.8f}".format(nm,angle,coeffs2[nm,0,0].real,coeffs2[nm,0,0].imag,tran2))
Tsum += tran1+tran2
print("tran:, {:.6f}".format(Tsum))
r_flux = mp.get_fluxes(refl_flux)
t_flux = mp.get_fluxes(tran_flux)
Rflux = -r_flux[0]/input_flux[0]
Tflux = t_flux[0]/input_flux[0]
print("flux:, {:.6f}, {:.6f}, {:.6f}".format(Rflux,Tflux,Rflux+Tflux))(Note: two sources with different polarizations are defined even though only one is used in the 2d simulation. This is necessary to ensure that all field components are used in the normalization run and not just those pertaining to the particular polarization of the source. Specifying extra_materials using an anisotropic permittivity was not sufficient in this regard.)
The output from this script is shown below. The values for the mode coefficients as well as the transmittance are non-zero only for the first order which is consistent with the reference. Also, the sum of the transmittance values from the mode coefficients (0.474278) is exactly half the transmittance obtained using the Poynting flux (0.947084).
mode1:, 0, 0.00, -0.07646-0.14027j, 0.00082210
mode2:, 0, 0.00, -0.04977-0.09338j, 0.00036070
mode1:, 1, 4.77, -2.34913+0.23330j, 0.17950932
mode2:, 1, 4.77, 3.00414-0.29894j, 0.29358311
mode1:, 2, 9.56, -0.00423+0.00761j, 0.00000244
mode2:, 2, 9.56, 0.00427+0.00115j, 0.00000063
mode1:, 3, 14.43, 0.00014-0.00003j, 0.00000000
mode2:, 3, 14.43, 0.00018-0.00010j, 0.00000000
mode1:, 4, 19.41, 0.00000+0.00000j, 0.00000000
mode2:, 4, 19.41, 0.00000-0.00000j, 0.00000000
tran:, 0.474278
flux:, 0.052870, 0.947084, 0.999954The results for a TE (EVEN_Z+ODD_Y) input are shown below and are also similar.
mode1:, 0, 0.00, -0.07608-0.14283j, 0.00036289
mode2:, 0, 0.00, 0.11645+0.21409j, 0.00082306
mode1:, 1, 4.77, -4.57737+0.44810j, 0.29313212
mode2:, 1, 4.77, -3.58358+0.35179j, 0.17967491
mode1:, 2, 9.56, -0.00646-0.00146j, 0.00000061
mode2:, 2, 9.56, -0.00649+0.01106j, 0.00000228
mode1:, 3, 14.43, -0.00028+0.00014j, 0.00000000
mode2:, 3, 14.43, 0.00021-0.00005j, 0.00000000
mode1:, 4, 19.41, -0.00000+0.00000j, 0.00000000
mode2:, 4, 19.41, 0.00000+0.00000j, 0.00000000
tran:, 0.473996
flux:, 0.052871, 0.947084, 0.999954Following your prescription for computing the mode coefficients of +y output for a (TM)+i(TE) circular input, the transmittance based on these mode coefficients (which involve halving the input flux in the denominator) correctly indicates that only the m=+1 diffraction order is non-zero. The following are the transmittance values for the first five orders:
0.00235
0.94126
0.00001
0.00000
0.00000Similarly, computing the mode coefficients of -y output for a (TM)-i(TE) circular input indicates that only the transmittance for the m=-1 diffraction order is non-zero:
0.00238
0.94348
0.00000
0.00000
0.00000Pending further testing, the results of the reference seem to therefore have been successfully reproduced.
oskooi
commented
5 years ago The following script which is adapted from the recent tutorial example computes the transmittance spectrum of a uniaxial homogeneous polarization grating for a circular-polarized planewave input (the tutorial involved a linear-polarized source). The calculation is based on the strategy outlined previously in this issue of doing two separate simulations using a linear-polarized planewave input with orthogonal polarizations (ODD_Z+EVEN_Y and EVEN_Z+ODD_Y) and then summing/diffing the mode coefficients appropriately.
import meep as mp
import math
import numpy as np
import argparse
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)
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
df = 0.05*fcen # frequency width
nmode = 5 # number of mode coefficients to compute
def pol_grating(d,gp,amp_ez,amp_hz):
sx = dpml+dsub+d+d+dpad+dpml
sy = gp
cell_size = mp.Vector3(sx,sy,0)
pml_layers = [mp.PML(thickness=dpml,direction=mp.X)]
def phi(p):
return math.pi*p.y/gp
# return the anisotropic permittivity tensor for a uniaxial, twisted nematic liquid crystal
def lc_mat(p):
# rotation matrix for rotation around x axis
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)]
# linear-polarized planewave pulse source
src_pt = mp.Vector3(-0.5*sx+dpml+0.3*dsub,0,0)
sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df), component=mp.Ez, center=src_pt, size=mp.Vector3(0,sy,0), amplitude=amp_ez),
mp.Source(mp.GaussianSource(fcen,fwidth=df), component=mp.Hz, center=src_pt, size=mp.Vector3(0,sy,0), amplitude=amp_hz)]
symmetries = []
if amp_ez == 1 and amp_hz == 0:
symmetries = [mp.Mirror(mp.Y,phase=+1)]
elif amp_ez == 0 and amp_hz == 1:
symmetries = [mp.Mirror(mp.Y,phase=-1)]
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
boundary_layers=pml_layers,
k_point=k_point,
sources=sources,
symmetries=symmetries,
extra_materials=[lc_mat(mp.Vector3())],
default_material=mp.Medium(index=n_0))
refl_pt = mp.Vector3(-0.5*sx+dpml+0.5*dsub,0,0)
refl_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,sy,0)))
sim.run(until_after_sources=100)
input_flux = mp.get_fluxes(refl_flux)
input_flux_data = sim.get_flux_data(refl_flux)
sim.reset_meep()
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
boundary_layers=pml_layers,
k_point=k_point,
symmetries=symmetries,
sources=sources,
geometry=geometry)
refl_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,sy,0)))
sim.load_minus_flux_data(refl_flux,input_flux_data)
tran_pt = mp.Vector3(0.5*sx-dpml-0.5*dpad,0,0)
tran_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=tran_pt, size=mp.Vector3(0,sy,0)))
sim.run(until_after_sources=300)
if amp_ez == 1 and amp_hz == 0:
res1 = sim.get_eigenmode_coefficients(tran_flux, range(1,nmode+1), eig_parity=mp.ODD_Z+mp.EVEN_Y)
res2 = sim.get_eigenmode_coefficients(tran_flux, range(1,nmode+1), eig_parity=mp.EVEN_Z+mp.EVEN_Y)
elif amp_ez == 0 and amp_hz == 1:
res1 = sim.get_eigenmode_coefficients(tran_flux, range(1,nmode+1), eig_parity=mp.ODD_Z+mp.ODD_Y)
res2 = sim.get_eigenmode_coefficients(tran_flux, range(1,nmode+1), eig_parity=mp.EVEN_Z+mp.ODD_Y)
angles1 = [math.degrees(math.acos(kdom.x/fcen)) for kdom in res1.kdom]
angles2 = [math.degrees(math.acos(kdom.x/fcen)) for kdom in res2.kdom]
for nm in range(nmode):
print("mode1:, {}, {:.2f}, {:.5f}{:+.5f}j".format(nm,angles1[nm],res1.alpha[nm,0,0].real,res1.alpha[nm,0,0].imag))
print("mode2:, {}, {:.2f}, {:.5f}{:+.5f}j".format(nm,angles2[nm],res2.alpha[nm,0,0].real,res2.alpha[nm,0,0].imag))
return input_flux[0], angles1, angles2, res1.alpha[:,0,0], res2.alpha[:,0,0];
if __name__ == '__main__':
parser = argparse.ArgumentParser()
parser.add_argument('-dd', type=float, default=1.7, help='chiral layer thickness (default: 1.7 μm)')
parser.add_argument('-gp', type=float, default=6.5, help='grating period (default: 6.5 μm)')
args = parser.parse_args()
input_flux1, angles1a, angles1b, coeffs1a, coeffs1b = pol_grating(args.dd,args.gp,1,0)
input_flux2, angles2a, angles2b, coeffs2a, coeffs2b = pol_grating(args.dd,args.gp,0,1)
tran1 = np.zeros(nmode)
tran2 = np.zeros(nmode)
tran1[0] = abs(coeffs1a[0]+1j*coeffs2b[0])**2/(input_flux1+input_flux2)
tran2[0] = abs(coeffs1a[0]-1j*coeffs2b[0])**2/(input_flux1+input_flux2)
tran1[1:] = abs(coeffs1a[1:]+coeffs1b[:-1]+1j*(coeffs2a[:-1]+coeffs2b[1:]))**2/(input_flux1+input_flux2)
tran2[1:] = abs(coeffs1a[1:]+coeffs1b[:-1]-1j*(coeffs2a[:-1]+coeffs2b[1:]))**2/(input_flux1+input_flux2)
print("tran:, 0, 0, {:.5f}".format(0.5*(tran1[0]+tran2[0])))
for m in range(1,nmode):
print("tran:, +{}, +{:.2f}, {:.5f}".format(m,angles1a[m],tran1[m]))
print("tran:, -{}, -{:.2f}, {:.5f}".format(m,angles1a[m],tran2[m]))Note: when adding the mode coefficients for non-zero orders, the arrays need to be offset from each other in order for the mode coefficients to correspond to the same angles.
The diffraction efficiency can be computed analytically: η0=cos2(πΔnd/λ), η±1=0.5*(1∓S)sin2(πΔnd/λ) where S is the ellipticity of the polarization (e.g., CCW: S=+1, CW: S=-1). We can verify the correctness of the script using two runs: (1) Δnd/λ=0.5, and (2) Δnd/λ=1.0.
Δnd/λ=0.5
tran:, 0, 0, 0.00066
tran:, +1, +4.77, 1.34654
tran:, -1, -4.77, 1.34253
tran:, +2, +9.56, 0.00003
tran:, -2, -9.56, 0.00002
tran:, +3, +14.43, 0.00002
tran:, -3, -14.43, 0.00002
tran:, +4, +19.41, 0.00002
tran:, -4, -19.41, 0.00002The results are incorrect: the transmittance for the m=+1 and -1 orders should be 0.94 and 0 (or vice versa) rather than 1.34654 and 1.34654. (the reflectance is 0.06 which is why the transmittance is not 1.0).
Δnd/λ=1.0
tran:, 0, 0, 0.92771
tran:, +1, +4.77, 0.00424
tran:, -1, -4.77, 0.00375
tran:, +2, +9.56, 0.00001
tran:, -2, -9.56, 0.00003
tran:, +3, +14.43, 0.00002
tran:, -3, -14.43, 0.00002
tran:, +4, +19.41, 0.00002
tran:, -4, -19.41, 0.00002The results are (mostly) correct: the m=0 order is approximately 0.93 and the m=±1 orders are 0 (although the results are still off slightly from their expected values).
stevengj
commented
5 years ago Two fixes:
Instead of an Hz source, use an Ey source, to make sure its phase and amplitude exactly correspond to a 90° rotation of the Ez source.
As I understand your code, the 1a and 2a coefficients correspond to the amplitude of the s polarization, while the 1b and 2b coefficients are the amplitude of the p polarization. The total power in a given diffraction order is |s|²+|p|². So, to get the total power in a diffraction order for circularly polarized input, you want |1a±i2a|² + |1b±i2b|²
oskooi
commented
5 years ago Even with the two suggested changes, the results are still incorrect for Δnd/λ=0.5 where only one of the +1 or -1 orders should be non-zero. Rather, as previously, the transmittance in both orders is non-zero and equal. For Δnd/λ=1, the results are correct (same as before).
Δnd/λ=0.5
tran:, 0, 0, 0.00064
tran:, +1, +4.77, 0.47147
tran:, -1, -4.77, 0.47147
tran:, +2, +9.56, 0.00001
tran:, -2, -9.56, 0.00001
tran:, +3, +14.43, 0.00001
tran:, -3, -14.43, 0.00001
tran:, +4, +19.41, 0.00001
tran:, -4, -19.41, 0.00001Δnd/λ=1.0
tran:, 0, 0, 0.92773
tran:, +1, +4.77, 0.00147
tran:, -1, -4.77, 0.00145
tran:, +2, +9.56, 0.00001
tran:, -2, -9.56, 0.00001
tran:, +3, +14.43, 0.00001
tran:, -3, -14.43, 0.00001
tran:, +4, +19.41, 0.00001
tran:, -4, -19.41, 0.00001import meep as mp
import math
import numpy as np
import argparse
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)
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
nmode = 5 # number of mode coefficients to compute
def pol_grating(d,gp,amp_ez):
sx = dpml+dsub+d+d+dpad+dpml
sy = gp
cell_size = mp.Vector3(sx,sy,0)
pml_layers = [mp.PML(thickness=dpml,direction=mp.X)]
def phi(p):
return math.pi*p.y/gp
# return the anisotropic permittivity tensor for a uniaxial, twisted nematic liquid crystal
def lc_mat(p):
# rotation matrix for rotation around x axis
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)]
# linear-polarized planewave pulse source
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 if amp_ez else mp.Ey, center=src_pt, size=mp.Vector3(0,sy,0))]
symmetries = []
if amp_ez:
symmetries = [mp.Mirror(mp.Y,phase=+1)]
else:
symmetries = [mp.Mirror(mp.Y,phase=-1)]
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
boundary_layers=pml_layers,
k_point=k_point,
sources=sources,
symmetries=symmetries,
default_material=mp.Medium(index=n_0))
sim.init_sim()
if amp_ez:
sim.fields.require_component(mp.Hz)
else:
sim.fields.require_component(mp.Ez)
refl_pt = mp.Vector3(-0.5*sx+dpml+0.5*dsub,0,0)
refl_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,sy,0)))
sim.run(until_after_sources=100)
input_flux = mp.get_fluxes(refl_flux)
input_flux_data = sim.get_flux_data(refl_flux)
sim.reset_meep()
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
boundary_layers=pml_layers,
k_point=k_point,
symmetries=symmetries,
sources=sources,
geometry=geometry)
refl_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,sy,0)))
sim.load_minus_flux_data(refl_flux,input_flux_data)
tran_pt = mp.Vector3(0.5*sx-dpml-0.5*dpad,0,0)
tran_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=tran_pt, size=mp.Vector3(0,sy,0)))
sim.run(until_after_sources=300)
if amp_ez:
res1 = sim.get_eigenmode_coefficients(tran_flux, range(1,nmode+1), eig_parity=mp.ODD_Z+mp.EVEN_Y)
res2 = sim.get_eigenmode_coefficients(tran_flux, range(1,nmode+1), eig_parity=mp.EVEN_Z+mp.EVEN_Y)
else:
res1 = sim.get_eigenmode_coefficients(tran_flux, range(1,nmode+1), eig_parity=mp.ODD_Z+mp.ODD_Y)
res2 = sim.get_eigenmode_coefficients(tran_flux, range(1,nmode+1), eig_parity=mp.EVEN_Z+mp.ODD_Y)
angles1 = [math.degrees(math.acos(kdom.x/fcen)) for kdom in res1.kdom]
angles2 = [math.degrees(math.acos(kdom.x/fcen)) for kdom in res2.kdom]
for nm in range(nmode):
print("mode1:, {}, {:.2f}, {:.5f}{:+.5f}j".format(nm,angles1[nm],res1.alpha[nm,0,0].real,res1.alpha[nm,0,0].imag))
print("mode2:, {}, {:.2f}, {:.5f}{:+.5f}j".format(nm,angles2[nm],res2.alpha[nm,0,0].real,res2.alpha[nm,0,0].imag))
return input_flux[0], angles1, angles2, res1.alpha[:,0,0], res2.alpha[:,0,0]
if __name__ == '__main__':
parser = argparse.ArgumentParser()
parser.add_argument('-dd', type=float, default=0.85, help='chiral layer thickness (default: 0.85 μm)')
parser.add_argument('-gp', type=float, default=6.5, help='grating period (default: 6.5 μm)')
args = parser.parse_args()
input_flux1, angles1a, angles1b, coeffs1a, coeffs1b = pol_grating(args.dd,args.gp,1)
input_flux2, angles2a, angles2b, coeffs2a, coeffs2b = pol_grating(args.dd,args.gp,0)
tran1 = np.zeros(nmode)
tran2 = np.zeros(nmode)
tran1[0] = (abs(coeffs1a[0])**2 + abs(coeffs2b[0])**2)/(input_flux1+input_flux2)
tran2[0] = (abs(coeffs1a[0])**2 + abs(coeffs2b[0])**2)/(input_flux1+input_flux2)
tran1[1:] = (abs(coeffs1a[1:]+1j*coeffs2a[:-1])**2 + abs(coeffs1b[:-1]+1j*coeffs2b[1:])**2)/(input_flux1+input_flux2)
tran2[1:] = (abs(coeffs1a[1:]-1j*coeffs2a[:-1])**2 + abs(coeffs1b[:-1]-1j*coeffs2b[1:])**2)/(input_flux1+input_flux2)
tran1 = 0.5*tran1
tran2 = 0.5*tran2
print("tran:, 0, 0, {:.5f}".format(tran1[0]+tran2[0]))
for m in range(1,nmode):
print("tran:, +{}, +{:.2f}, {:.5f}".format(m,angles1a[m],tran1[m]))
print("tran:, -{}, -{:.2f}, {:.5f}".format(m,angles1a[m],tran2[m]))Δnd/λ=0.5 Ez source
mode1:, 0, 0.00, -0.12669-0.03749j
mode2:, 0, 4.77, -1.97321+3.20535j
mode1:, 1, 4.77, 1.98577-3.20884j
mode2:, 1, 9.56, 0.00229+0.00413j
mode1:, 2, 9.56, 0.00897-0.02405j
mode2:, 2, 14.43, 0.00001-0.00012j
mode1:, 3, 14.43, 0.01153-0.01964j
mode2:, 3, 19.41, -0.00000-0.00000j
mode1:, 4, 19.41, 0.01155-0.01976j
mode2:, 4, 24.54, -0.00000+0.00000jHz source
mode1:, 0, 4.77, 1.98027-3.19884j
mode2:, 0, 0.00, 0.14403+0.00967j
mode1:, 1, 9.56, -0.00252-0.00428j
mode2:, 1, 4.77, 1.97990-3.21514j
mode1:, 2, 14.43, -0.00002+0.00012j
mode2:, 2, 9.56, 0.01046-0.02367j
mode1:, 3, 19.41, 0.00000+0.00000j
mode2:, 3, 14.43, 0.01267-0.01947j
mode1:, 4, 24.54, 0.00000-0.00000j
mode2:, 4, 19.41, 0.01268-0.01959jThe polarization grating has been correctly set up in the script. This is verified by results shown in the figure below for snapshots of Ez for 4 different cases involving a circularly-polarized input planewave: Δnd/λ=0.5 and 1.0, and input chirality of Ez+iEy and Ez-iEy. As expected, the chirality of the input planewave determines the first diffraction order (m=±1) for Δnd/λ=0.5. The angle of the first diffracted order (±4.77°) also agrees with the expected analytic result. Results are similar for Ey snapshots.
Is it possible that there is a bug in mode decomposition for the case of a 2d cell where all field components are used?
oskooi
commented
5 years ago It seems the correct results are produced by modifying the transmittance expression as follows:
original
tran1[1:] = (abs(coeffs1a[1:]+1j*coeffs2a[:-1])**2 + abs(coeffs1b[:-1]+1j*coeffs2b[1:])**2)/(input_flux1+input_flux2)
tran2[1:] = (abs(coeffs1a[1:]-1j*coeffs2a[:-1])**2 + abs(coeffs1b[:-1]-1j*coeffs2b[1:])**2)/(input_flux1+input_flux2)modified
tran1[1:] = (abs(coeffs1a[1:]+coeffs2a[:-1])**2 + abs(coeffs1b[:-1]-coeffs2b[1:])**2)/(input_flux1+input_flux2)
tran2[1:] = (abs(coeffs1a[1:]-coeffs2a[:-1])**2 + abs(coeffs1b[:-1]+coeffs2b[1:])**2)/(input_flux1+input_flux2)The main change is that the 1j prefactor multiplying the mode coefficients with ODD_Z parity for an Ey linearly-polarized input has been removed.
Δnd/λ=0.5
tran:, 0, 0, 0.00064
tran:, +1, +4.77, 0.94294
tran:, -1, -4.77, 0.00000
tran:, +2, +9.56, 0.00001
tran:, -2, -9.56, 0.00001
tran:, +3, +14.43, 0.00001
tran:, -3, -14.43, 0.00001
tran:, +4, +19.41, 0.00001
tran:, -4, -19.41, 0.00001Note that only the m=+1 order has non-zero transmittance.
Δnd/λ=1.0
tran:, 0, 0, 0.92773
tran:, +1, +4.77, 0.00001
tran:, -1, -4.77, 0.00291
tran:, +2, +9.56, 0.00001
tran:, -2, -9.56, 0.00001
tran:, +3, +14.43, 0.00001
tran:, -3, -14.43, 0.00001
tran:, +4, +19.41, 0.00001
tran:, -4, -19.41, 0.00001Note that only the m=0 order has non-zero transmittance.
As another verification, changing the amplitude of the Ey source from +1 to -1 which in effect flips the chirality of the input source results in the m=-1 having non-zero transmittance.
tran:, 0, 0, 0.00064
tran:, +1, +4.77, 0.00000
tran:, -1, -4.77, 0.94294
tran:, +2, +9.56, 0.00001
tran:, -2, -9.56, 0.00001
tran:, +3, +14.43, 0.00001
tran:, -3, -14.43, 0.00001
tran:, +4, +19.41, 0.00001
tran:, -4, -19.41, 0.00001Maybe the issue is that MPB is not choosing the phase for the Ez and Hz polarized modes in the same way? (Once we implement the set_planewave feature you will be able to choose whatever phase you want.) That is, the phase is chosen deterministically, but it is not necessarily chosen in the same way for different polarizations.
If you call get-eigenmode-field (or whatever it is called) for Ez at the center of one mode and for Hz at the center of the other mode, you will see the relative phases. In particular, look at the ratio Hz/Ez — if that isn't one, then the phases are being chosen differently. (I'm guessing from above that they differ by a factor of ±i.) Then you could correct the mode amplitude α of the Hz mode by multiplying α by Ez/Hz, I think.
oskooi
commented
5 years ago The ratio Hz/Ez is not one as shown in the figure below (bottom right inset). The imaginary part of Hz/Ez is 0 but the real part is a non-zero constant (-1.5). These results were generated using the accompanying script which computes two separate eigenmodes (Ez:ODD_Z+EVEN_Y and Hz:EVEN_Z+ODD_Y) for an arbitrary band (5) of a dielectric medium with n=1.5 at a k_point of (0,0,0). The dominant planewave for both eigenmodes is (1.445683,-0.400000,0.000000). The other three insets show the real and imaginary parts as well as the intensity profile of these two eigenmodes. Results are similar for other bands.
import meep as mp
import matplotlib.pyplot as plt
import numpy as np
resolution = 50 # pixels/μm
dpml = 1 # PML thickness
dmat = 10 # length of material
sx = dpml+dmat+dpml
sy = 10
cell_size = mp.Vector3(sx,sy,0)
pml_layers = [mp.Absorber(thickness=dpml,direction=mp.X)]
fcen = 1 # center frequency (wavelength = 1 μm)
band_num = 5
glass = mp.Medium(index=1.5)
k = mp.Vector3(0,0,0)
eig_parity = mp.ODD_Z+mp.EVEN_Y
symmetries = [mp.Mirror(mp.Y)]
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
boundary_layers=pml_layers,
k_point=k,
symmetries=symmetries,
default_material=glass)
sim.init_sim()
mon_pt = mp.Vector3(0.5*sx-dpml-0.2*dmat,0,0)
data_oddz = sim.get_eigenmode(fcen, mp.X, mp.Volume(center=mon_pt,size=mp.Vector3(0,sy,0)), band_num, k, parity=eig_parity, resolution=resolution, verbose=True)
sim.reset_meep()
eig_parity = mp.EVEN_Z+mp.ODD_Y
symmetries = [mp.Mirror(mp.Y, phase=-1)]
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
boundary_layers=pml_layers,
k_point=k,
symmetries=symmetries,
default_material=glass)
sim.init_sim()
data_evenz = sim.get_eigenmode(fcen, mp.X, mp.Volume(center=mon_pt,size=mp.Vector3(0,sy,0)), band_num, k, parity=eig_parity, resolution=resolution, verbose=True)
y = np.linspace(-0.5*sy,0.5*sy,sy*resolution)
Ez_data = np.asarray([data_oddz.amplitude(mp.Vector3(0,yc,0),mp.Ez) for yc in y])
Hz_data = np.asarray([data_evenz.amplitude(mp.Vector3(0,yc,0),mp.Hz) for yc in y])
phase = Hz_data/Ez_data
plt.figure(dpi=300)
plt.subplot(2,2,1)
plt.plot(y,np.real(Ez_data),'bo-',label=r'real (E$_z$)')
plt.plot(y,np.real(Hz_data),'ro-',label=r'real (H$_z$)')
plt.xlabel('y coordinate')
plt.ylabel(r'E$_z$ or H$_z$')
plt.legend(loc='lower right')
plt.title('real part')
plt.subplot(2,2,2)
plt.plot(y,np.imag(Ez_data),'bo-',label=r'imag (E$_z$)')
plt.plot(y,np.imag(Hz_data),'ro-',label=r'imag (H$_z$)')
plt.xlabel('y coordinate')
plt.ylabel(r'E$_z$ or H$_z$')
plt.legend(loc='lower right')
plt.title('imaginary part')
plt.subplot(2,2,3)
plt.plot(y,abs(Ez_data)**2,'bo-',label=r'E$_z$')
plt.plot(y,abs(Hz_data)**2,'ro-',label=r'H$_z$')
plt.xlabel('y coordinate')
plt.ylabel(r'|E$_z$|$^2$ or |H$_z$|$^2$')
plt.legend(loc='lower right')
plt.title('intensity')
plt.subplot(2,2,4)
plt.plot(y,np.real(phase),'bo-',label='real')
plt.plot(y,np.imag(phase),'ro-',label='imag')
plt.xlabel('y coordinate')
plt.ylabel('H$_z$/E$_z$')
plt.legend(loc='center')
plt.title('relative phase')
plt.tight_layout()
plt.show()The mode coefficients for the first band (i.e., zeroth order) of 2d vacuum for a Ez and Ey source have a phase difference of exactly π (i.e., the coefficients differ by a scalar factor of -1). This is demonstrated using the following script.
import meep as mp
dpml = 1
dair = 10
sx = dpml+dair+dpml
sy = 5
cell_size = mp.Vector3(sx,sy,0)
pml_layers = [mp.PML(thickness=dpml,direction=mp.X)]
resolution = 25
fcen = 1.0
src_pt = mp.Vector3(-0.5*sx+dpml+0.1*dair)
mon_pt = mp.Vector3(+0.5*sx-dpml-0.1*dair)
k = mp.Vector3(0,0,0)
def compute_zero_difforder(src_ez):
sources = [mp.Source(mp.GaussianSource(fcen, 0.2*fcen),
center=src_pt,
size=mp.Vector3(0,sy,0),
component=mp.Ez if src_ez else mp.Ey)]
symmetries = [mp.Mirror(direction=mp.Y, phase=+1 if src_ez else -1)]
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
boundary_layers=pml_layers,
symmetries=symmetries,
k_point=k,
sources=sources)
mode_mon = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=mon_pt, size=mp.Vector3(0,sy,0)))
sim.run(until_after_sources=100)
res = sim.get_eigenmode_coefficients(mode_mon, [1], eig_parity=mp.ODD_Z+mp.EVEN_Y if src_ez else mp.EVEN_Z+mp.ODD_Y)
coeffs = res.alpha
print("mode:, {}, {:.5f}{:+.5f}j, {:.5f}".format("ez" if src_ez else "ey",coeffs[0,0,0].real,coeffs[0,0,0].imag,abs(coeffs[0,0,0])**2))
if __name__ == '__main__':
compute_zero_difforder(True)
compute_zero_difforder(False)The output displays the mode coefficient for each source as well as its squared magnitude (i.e., power):
mode:, ez, 0.10776+0.00400j, 0.01163
mode:, ey, -0.10776-0.00400j, 0.01163Does this eiπ phase get introduced somewhere in get_eigenmode?
stevengj
commented
5 years ago Okay, good to know — everything is clear and correct, and the UI issue will be fixed by the set-planewave functionality.
The mode decomposition feature currently only works for linearly-polarized modes. It would be useful to extend this to circularly-polarized modes. An example application which would be nice to set up as a tutorial in the documentation is described in Optics Letters, Vol. 33, No. 20, pp. 2287-9, 2008 (pdf). This involves a linearly-polarized planewave incident on a 2d slab of cholesteric liquid crystal. The transmitted fields are diffracted planewaves at various angles with left or right circular polarization. The mode coefficients of these circularly-polarized modes could be used to calculate the diffraction efficiency of each order.