Extract the electric field and angle by orders

[zhixwang]

Hi authors, can the tool extract the electric field by each order? For example, the xy-field or xz-field of order (0,1), (1,1),...

I also want to extract the angle of each order within a particular layer. There is a function in rcwa.py called diffraction_angle(), however I didn't find any introduction or explanation for this function in the manual.

I want to do some follow-up computation by orders, using both the field and the angle of it.

Thanks!

[kch3782]

Hi.

Currently, the tool does not provide such a function. Although S-parameter function does not plot the field, it has similar role for unit amplitude plane wave incidence. I think it is possible to plot using S-parameter and k-vector according to diffraction angle. When the order is [1,0], the transmission S-parameter is S, and the k-vector normalized for k0 is [kx, ky, kz], then the fields can be calculated in the form of S pol_amplitude exp(1.j k0 (kxx+kyy+kz*z)). Here, it will be necessary to adjust the amplitude according to the polarization. Another method is to fill S-matrix to zero except the order you want. Then calculate field. (order_indices could be found using an internal function sim._matching_indices(orders), and in the case of y-pol, the total number of orders must be added to order_indices.)

sim.solve_global_smatrix()
temp_S = sim.S[0]
sim.S[0] = torch.zeros_like(temp_S)
sim.S[0][order_indices,ref_order_indices] = temp_S[order_indices,ref_order_indices]
# source settings and field plot

However, since the sim.S is not based on p-pol, s-pol notation, and if x-pol and y-pol are not parallel to the p-pol and s-pol, it will not be an accurate plot.

Also, the diffraction_angle function can be used as follows. If you want to calculate the angles for the diffraction order [1,0] the output layer, You can use it with inc, azi = sim.diffraction_angle([1,0],layer='output',unit='radian') inc is the inclination angle and azi is the azimuthal angle.

diffraction_angle(self,orders,*,layer='output',unit='radian')
    Parameters
    - orders: selected diffraction orders (Recommended shape: Nx2)
    - layer: selected layer ('i', 'in', 'input' / 'o', 'out', 'output')
    - unit: unit of the output angles ('r', 'rad', 'radian' / 'd', 'deg', 'degree')
    Return
    - inclination angle (torch.Tensor), azimuthal angle (torch.Tensor)

Thank you.

[zhixwang]

Dear Author,

Many thanks for your kind and detailed reply. I am trying to follow up according to your ideas.

Please allow me to explain my simulation situation more clearly. Let's say I have 3 internal layers (layer0, layer1, layer2). Now I want to compute the E-field of a certain order, let's say order [1, 2], in layer 2. Physically as long as I can obtain the complex amplitude of the field at this order E_tilda(x,y)[1,2], and the wavevector of this order K_[1,2], I should be able to derive it as:

E_x[1,2] = E_tilda,x[1,2] exp(1jK_x[1,2]x) exp(1jK_y[1,2]y) exp(1jK_z[1,2]z) (time coherent term). E_y[1,2] = E_tilda,y[1,2] exp(1jK_x[1,2]x) exp(1jK_y[1,2]y) exp(1jK_z[1,2]z) (time coherent term).

If I understand your answer correctly, it is possible to directly get E_tilda(x,y)[1,2] and K_[1,2]. Is it correct if I compute it as follows?

K_x[1,2], K_y[1,2] = K0 sin(inc) cos(azi), K0 sin(inc) sin(azi)

, where inc, azi in layer2 come from sim.diffraction_angle([1,2],layer=2,unit='radian')

I didn't fully understand how to get E_tilda,x[1,2] and E_tilda,y[1,2] from your answers. Could you please share more about it in this regard?

Many thanks!

[kch3782]

Hi.

First, the diffraction order refers to each order when the electromagnetic field of the layer is Fourier-expanded. In addition, eigenmode refers to the electromagnetic field mode with a specific eigenvalue kz inside each layer. This eigenmode field is composed of EM fields with several diffraction orders. The eigenmodes of inhomogeneous layers have several diffraction orders, but the eigenmodes consist of single diffraction order in homogenous and isotropic layers. Therefore, several eigenmodes with different eigenvalues (kz) are involved in a specific diffraction order in the inhomogeneous layer. Thus, the field for a specific order cannot be calculated in the same way as in the previous answer for an inhomogeneous internal layer.

Is internal layer 2 in your example an inhomogenous layer? Then, is the purpose to obtain an electromagnetic field for a specific eigenmode, or to obtain a field for a specific diffraction order?

[zhixwang]

Dear Author,

Many thanks for your reply.

1. The internal layer 2 in my example is a homogenous layer.

2. The order [1, 2] in my example, is referring to the "order" in the sense of physical diffraction: [0,0] is "straight" transmission, [0,1], [1,0], [-1,0], [0.-1] share the smallest and the same diffraction angle (the same Kz) in case of a square unit-cell lattice, etc.

3. My final goal here is to extract the electromagnetic field for pairs of eigenmodes with the same Kz in a given homogeneous layer. For example, in the case of a square lattice: E_x[0,1]+E_x[0,-1]+E_x[1,0]+E_x[-1,0]. To do so, I think it's straightforward to extract the field of each order first and add them later.

I don't know how to extract the correct coefficients correctly.

I hope my explanations are clear enough for you. Thanks a lot!

[zhixwang]

Another way of explanation would be like this:

The E-field in an internal homogenous layer is a superposition of all orders of plane waves, with different complex amplitude and diffraction angles. How to extract the complex amplitude E_x_mn, E_y_mn, E_z_mn of each order of plane wave?

[kch3782]

The below codes are extracted from torcwa.rcwa.field_xy()

z_prop = 0.

if sim.source_direction == 'forward':
    C = torch.matmul(sim.C[0][layer_num],sim.E_i)
elif sim.source_direction == 'backward':
    C = torch.matmul(sim.C[1][layer_num],sim.E_i)

x_axis = x_axis.reshape([-1,1,1])
y_axis = y_axis.reshape([1,-1,1])

Kx_norm, Ky_norm = sim.Kx_norm, sim.Ky_norm
kz_norm = sim.kz_norm[layer_num]
E_eigvec = sim.E_eigvec[layer_num]
H_eigvec = sim.H_eigvec[layer_num]

Cp = torch.diag(C[:2*sim.order_N,0])
Cm = torch.diag(C[2*sim.order_N:,0])

eps_conv_inv = torch.linalg.inv(sim.eps_conv[layer_num])
mu_conv_inv = torch.linalg.inv(sim.mu_conv[layer_num])

# Phase
z_phase_p = torch.diag(torch.exp(1.j*sim.omega*kz_norm*z_prop))
z_phase_m = torch.diag(torch.exp(1.j*sim.omega*kz_norm*(sim.thickness[layer_num]-z_prop)))

# Fourier domain fields
# [diffraction order, eigenmode number]
Exy_p = torch.matmul(E_eigvec,z_phase_p)
Ex_p = Exy_p[:sim.order_N,:]
Ey_p = Exy_p[sim.order_N:,:]
Hz_p = torch.matmul(mu_conv_inv,torch.matmul(Kx_norm,Ey_p)) - torch.matmul(mu_conv_inv,torch.matmul(Ky_norm,Ex_p))
Exy_m = torch.matmul(E_eigvec,z_phase_m)
Ex_m = Exy_m[:sim.order_N,:]
Ey_m = Exy_m[sim.order_N:,:]
Hz_m = torch.matmul(mu_conv_inv,torch.matmul(Kx_norm,Ey_m)) - torch.matmul(mu_conv_inv,torch.matmul(Ky_norm,Ex_m))
Hxy_p = torch.matmul(H_eigvec,z_phase_p)
Hx_p = Hxy_p[:sim.order_N,:]
Hy_p = Hxy_p[sim.order_N:,:]
Ez_p = torch.matmul(eps_conv_inv,torch.matmul(Ky_norm,Hx_p)) - torch.matmul(eps_conv_inv,torch.matmul(Kx_norm,Hy_p))
Hxy_m = torch.matmul(-H_eigvec,z_phase_m)
Hx_m = Hxy_m[:sim.order_N,:]
Hy_m = Hxy_m[sim.order_N:,:]
Ez_m = torch.matmul(eps_conv_inv,torch.matmul(Ky_norm,Hx_m)) - torch.matmul(eps_conv_inv,torch.matmul(Kx_norm,Hy_m))

# Complex amplitude at z_prop from lower boundary
Ex_mn = torch.sum(torch.matmul(Ex_p,Cp) + torch.matmul(Ex_m,Cm),dim=1)
Ey_mn = torch.sum(torch.matmul(Ey_p,Cp) + torch.matmul(Ey_m,Cm),dim=1)
Ez_mn = torch.sum(torch.matmul(Ez_p,Cp) + torch.matmul(Ez_m,Cm),dim=1)
Hx_mn = torch.sum(torch.matmul(Hx_p,Cp) + torch.matmul(Hx_m,Cm),dim=1)
Hy_mn = torch.sum(torch.matmul(Hy_p,Cp) + torch.matmul(Hy_m,Cm),dim=1)
Hz_mn = torch.sum(torch.matmul(Hz_p,Cp) + torch.matmul(Hz_m,Cm),dim=1)

The last 6 variables are complex amplitude of the electromagnetic field at z_prop from the lower boundary in layer_num. If you want to extract the complex amplitude for [0,1] and [1,0] orders, then you can use below code.

order_indices = sim._matching_indices(torch.tensor([[0,1],[1,0]]))
Ex_mn_[order_indices] = Ex_mn[order_indices]

Or if you want to plot the field only for these orders, you can use below code with some modification.

####
Ex_mn_ = torch.zeros_like(Ex_mn)
order_indices = sim._matching_indices(torch.tensor([[0,1],[1,0]]))
Ex_mn_[order_indices] = Ex_mn[order_indices]

# Spatial domain fields
xy_phase = torch.exp(1.j * sim.omega * (sim.Kx_norm_dn*x_axis + sim.Ky_norm_dn*y_axis))
Ex = torch.sum(Ex_mn_.reshape(1,1,-1)*xy_phase,dim=2)
Ey = torch.sum(Ey_mn_.reshape(1,1,-1)*xy_phase,dim=2)
Ez = torch.sum(Ez_mn_.reshape(1,1,-1)*xy_phase,dim=2)
Hx = torch.sum(Hx_mn_.reshape(1,1,-1)*xy_phase,dim=2)
Hy = torch.sum(Hy_mn_.reshape(1,1,-1)*xy_phase,dim=2)
Hz = torch.sum(Hz_mn_.reshape(1,1,-1)*xy_phase,dim=2)

Also, there are positive mode and negative mode. (positive mode is propagate or decay in +z direction) So, if you want to see only positive mode, delete + torch.matmul(Ex_m,Cm) in the part where six variables are defined.

I hope it will be of help.

[zhixwang]

Dear author,

Many thanks for your detailed answer. It is significantly helpful.

Best wishes,