Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Normalization of intensity integrated on all orders #19

Open
GiuseppeELio opened this issue Mar 16, 2023 · 4 comments
Open

Normalization of intensity integrated on all orders #19

GiuseppeELio opened this issue Mar 16, 2023 · 4 comments

Comments

@GiuseppeELio
Copy link

Dear authors, thanks for releasing this very nice software. I tried to launch a quick simulation to understand how the code handles normalization when integrating over all orders, but I'm getting extremely large results in transmission. I've tried increasing the number of harmonics, playing with ref_order and the evanescent parameter, or also increasing the precision to sim_dtype = torch.complex128 and geo_dtype = torch.float64, but none of these seem to help. Do you have any suggestions?
Also, when the simulations wavelength matches the structure period, I see that the torcwa returns zero diffraction efficiencies.
Please, see the script below.
Thank you for your help!

import numpy as np
import torch
import matplotlib.pyplot as plt
from tqdm import tqdm

import torcwa

torch.backends.cuda.matmul.allow_tf32 = False
sim_dtype = torch.complex64
geo_dtype = torch.float32
device = torch.device('cuda')

lam_min = 400.
lam_max = 600.
lam_step = 2.
DL = int((lam_max-lam_min)/lam_step)
lam_arr = torch.linspace(lam_min, lam_max, DL, dtype=geo_dtype, device=device)

nx, ny, edge_sharpness = 256, 256, 128
L = [400., 400.] # the period coincides by chance with the first wavelength ... 
geom = torcwa.geometry(Lx=L[0], Ly=L[1], nx=nx, ny=ny, edge_sharpness=edge_sharpness, dtype=geo_dtype, device=device)

# define layer
R = 120.
layer_geom = geom.circle(R=R, Cx=L[0]/2.,Cy=L[1]/2.)
layer_thickness = 400.

substrate_eps = 1.5**2

order_N = 10
order = [order_N, order_N]

# define array with all orders
o = np.arange(-order_N, order_N+1)
all_orders = np.stack(np.meshgrid(o,o),-1).reshape(-1,2)

T_sum = torch.zeros_like(lam_arr)
R_sum = torch.zeros_like(lam_arr)

for lam_idx in tqdm(range(len(lam_arr))):
  
   sim = torcwa.rcwa(freq=1.0/lam_arr[lam_idx], order=order, L=L, dtype=sim_dtype, device=device)
   
   sim.add_input_layer(eps=1.0)
   sim.add_output_layer(eps=substrate_eps)
   sim.set_incident_angle(inc_ang=0.0, azi_ang=0.0) # normal incidence

   layer_eps = layer_geom*substrate_eps + (1.-layer_geom)
   sim.add_layer(thickness=layer_thickness, eps=layer_eps)
   
   sim.solve_global_smatrix()

   rxx = sim.S_parameters(orders=all_orders, direction='f', port='r', polarization='xx')
   txx = sim.S_parameters(orders=all_orders, direction='f', port='t', polarization='xx')
   ryy = sim.S_parameters(orders=all_orders, direction='f', port='r', polarization='yy')
   tyy = sim.S_parameters(orders=all_orders, direction='f', port='t', polarization='yy')
   
   T_sum[lam_idx] = (torch.sum(torch.abs(txx)**2) + torch.sum(torch.abs(tyy)**2))/2
   R_sum[lam_idx] = (torch.sum(torch.abs(rxx)**2) + torch.sum(torch.abs(ryy)**2))/2

eps_out, _ = sim.return_layer(layer_num=0, nx=geom.nx, ny=geom.ny)

fig, axarr = plt.subplots(1,2, figsize=(16,4))
eps_im = axarr[0].imshow(eps_out.cpu().real, extent=[-L[0]/2, L[0]/2, -L[1]/2, L[1]/2])
fig.colorbar(eps_im, ax=axarr[0])
axarr[0].set(title=f'$\epsilon$ of layer #0 (FO: {order_N})',
             xlabel='x [nm]', ylabel='y [nm]')
axarr[1].plot(lam_arr.cpu(), R_sum.cpu(), label='R')  
axarr[1].plot(lam_arr.cpu(), T_sum.cpu(), label='T')
axarr[1].set(title=f'spectrum (FO: {order_N})',
             xlabel='$\lambda$ [nm]',
             ylabel='normalized(?) intensity')
axarr[1].legend()
plt.show()

Which gives
image

@kch3782
Copy link
Owner

kch3782 commented Mar 19, 2023

Hi.

Firstly, there is an error in logic for judging the evanescent field when calculating the S-parameters of the transmission part. Thank you for identifying this issue. We will address it in version 0.1.4.2.

Even after updating the code, you may still get a T-value greater than 1 for wavelengths of around 430 nm or less due to differences in notation. Please refer to the figure below:
image

In the torcwa.rcwa.S_parameter() function, the polarization argument supports both xy-pol notation and ps-pol notation. The left side of the figure represents the E-field direction used for p-pol/s-pol notation, while the right side represents the E-field direction used for x-pol/y-pol notation. The function returns the ratio of the magnitudes of Ep and Es in the case of the left side, and returns the normalized value with power if the power_norm option is used. In ps-pol notation, k-vector, Ep, and Es are all orthogonal, so power beyond 1 does not occur.
On the other hand, in the case of the right side, the function returns the ratio of the magnitudes of Ex and Ey, and the power_norm option returns the normalized value. For power normalization of xy-pol notation, the Ez component is utilized for normalization. The value of Ez is determined by the k-vector, Ex, and Ey, and if both Ex and Ey affect Ez, the value of the S-parameter may exceed 1. In this simulation, there are diffraction orders of [1,1], [1,-1], [-1,1], and [-1,-1] below wavelength of around 430 nm for example. Although the magnitude of each power-normalized ratio of Ex and Ey may exceed 1, the contributions of Ez due to Ex and Ey are offset by the phase difference between the two components of Ex and Ey. As a result, the overall size of the actual Ez decreases, leading to energy conservation. This does not occur with ps-pol notation, so it is recommended to use p-pol and s-pol notation in the case of diffraction orders with both kx and ky component or oblique incidence.

You can use polarization='pp' or 'ss' instead of polarization='xx' in the current code.

Additionally, when the grating period and wavelength are the same, a divergence phenomenon occurs due to the Wood anomaly, caused by singular matrix operations by the mode with kz = 0. In this case, it is recommended to simulate with a small deviation in the wavelength.

Thank you and I hope this would be helpful.

@GiuseppeELio
Copy link
Author

Thank you for your feedback.
Switching to the 'pp' and 'ss' notation avoids transmittance value above one. However, we still find that energy conservation is not satisfied over the previous wavelength range.
image

For instance, at a wavelength of 460 nm, there is still a 10% of energy missing which we cannot explain.
We tried to replicate these results using another Torch based code called nannos, which also using the traditional RCWA implementation, checking for convergence over multiple Fourier orders, but it doesn't seem the problem here.
image

Do you have any suggestions? Let me know if you need the nannos script for the comparison.

@kch3782
Copy link
Owner

kch3782 commented Mar 24, 2023

Hi.

I think the diffraction efficiencies of the converted polarizations are not included in the T_sum and R_sum equations. pp and ss components are considered in the formula, but there are some polarization conversion like s-pol to p-pol (ps) and p-pol to s-pol (sp) for diffraction orders other than [0,0] order. This components should be included in the diffraction efficiency to check energy conservation.
Then, the diffraction efficiency can be calculated like the code below. There are three cases. (input p-pol, input s-pol, and the average of the two cases)

rpp = sim.S_parameters(orders=all_orders, direction='f', port='r', polarization='pp')
tpp = sim.S_parameters(orders=all_orders, direction='f', port='t', polarization='pp')
rss = sim.S_parameters(orders=all_orders, direction='f', port='r', polarization='ss')
tss = sim.S_parameters(orders=all_orders, direction='f', port='t', polarization='ss')
rps = sim.S_parameters(orders=all_orders, direction='f', port='r', polarization='ps')
tps = sim.S_parameters(orders=all_orders, direction='f', port='t', polarization='ps')
rsp = sim.S_parameters(orders=all_orders, direction='f', port='r', polarization='sp')
tsp = sim.S_parameters(orders=all_orders, direction='f', port='t', polarization='sp')

# Diffraction efficiency of input p-polarization
T_sum[lam_idx] = (torch.sum(torch.abs(tpp)**2) + torch.sum(torch.abs(tsp)**2))
R_sum[lam_idx] = (torch.sum(torch.abs(rpp)**2) + torch.sum(torch.abs(rsp)**2))
# or Diffraction efficiency of input s-polarization
T_sum[lam_idx] = (torch.sum(torch.abs(tss)**2) + torch.sum(torch.abs(tps)**2))
R_sum[lam_idx] = (torch.sum(torch.abs(rss)**2) + torch.sum(torch.abs(rps)**2))
# or Average of diffraction efficiency
T_sum[lam_idx] = ((torch.sum(torch.abs(tpp)**2) + torch.sum(torch.abs(tsp)**2)) + (torch.sum(torch.abs(tss)**2) + torch.sum(torch.abs(tps)**2)))/2
R_sum[lam_idx] = ((torch.sum(torch.abs(rpp)**2) + torch.sum(torch.abs(rsp)**2)) + (torch.sum(torch.abs(rss)**2) + torch.sum(torch.abs(rps)**2)))/2

When calculated in this way, the value of R+T is 1, resulting in energy conservation.
image

Thank you.

@GiuseppeELio
Copy link
Author

Thank you for your advice. Now the right syntax is clear to me.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants