Unexpected field enhancement trend in terajet simulation

[SebastianD9]

Hello,

I am simulating a terajet formed by a dielectric cube (3D), and I am observing a result that seems inconsistent with the literature.

According to published results on terajets, increasing the refractive index of the particle should lead to higher field enhancement:

However, in my simulations I observe the opposite trend:
For refractive index n = 1.41, the field enhancement agrees well with literature values.
For n = 2.1, the enhancement is significantly lower — even smaller than for n = 1.41.

This seems unphysical and suggests that there might be an issue in my setup.

Has anyone encountered a similar issue in Meep or terajet simulations?

Any suggestions would be greatly appreciated.

Simulation details:
Wavelength: 1.5 mm (≈ 0.2 THz)
Geometry: dielectric cube (size = 2λ)
Source: continuous wave Gaussian beam (Ex polarization)
Normalization: field intensity normalized to reference simulation (no cube)
Quantity analyzed: max(E²) behind the cube

Code:

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

lambda_mm = 1.5  #1.5mm = 0.2THz
fcen = 1.0 / lambda_mm  
n_cube = 2.1   
L_cube = 2 * lambda_mm  
resolution = 20  

dpml = 1.5 * lambda_mm  
z_before = 3.0 * lambda_mm   
z_after = 8.0 * lambda_mm    
sx_total = 10.0 * lambda_mm  
sy_total = 10.0 * lambda_mm  
sz_total = z_before + L_cube + z_after + 2 * dpml  

cell = mp.Vector3(sx_total, sy_total, sz_total) 
z_cube_center = -sz_total / 2 + dpml + z_before + L_cube / 2 
z_shadow_meep = z_cube_center + L_cube / 2  

z_src = -sz_total / 2 + dpml + 0.5 * lambda_mm  
beam_width = 3.0 * lambda_mm  

n_periods = 30  
t_end = sz_total * 2.0 + n_periods / fcen + L_cube * n_cube

def gaussian_beam(p):
    return np.exp(-((p.x / beam_width)**2 + (p.y / beam_width)**2))

def get_sources():
    return [
        mp.Source(
            src=mp.ContinuousSource(frequency=fcen, is_integrated=True),  
            component=mp.Ex,  
            center=mp.Vector3(0, 0, z_src),  
            size=mp.Vector3(sx_total, sy_total, 0),  
            amp_func=gaussian_beam  
        )
    ]

print("reference")
sim_ref = mp.Simulation(
    cell_size=cell,
    geometry=[],
    sources=get_sources(),
    boundary_layers=[mp.PML(thickness=dpml)],
    resolution=resolution,
    dimensions=3
)

dft_ref = sim_ref.add_dft_fields(
    [mp.Ex, mp.Ey], fcen, fcen, 1,
    center=mp.Vector3(),
    size=mp.Vector3(sx_total, 0, sz_total)  
)

sim_ref.run(until=t_end)

ex_ref = sim_ref.get_dft_array(dft_ref, mp.Ex, 0)
ey_ref = sim_ref.get_dft_array(dft_ref, mp.Ey, 0)
E2_ref_array = np.abs(ex_ref)**2 + np.abs(ey_ref)**2

E2_ref_val = np.max(E2_ref_array)  

sim_ref.reset_meep()

geometry = [
    mp.Block(
        size=mp.Vector3(L_cube, L_cube, L_cube),
        center=mp.Vector3(0, 0, z_cube_center),
        material=mp.Medium(index=n_cube, D_conductivity=0)  
    )
]

print("main")
sim_main = mp.Simulation(
    cell_size=cell,
    geometry=geometry,
    sources=get_sources(),
    boundary_layers=[mp.PML(thickness=dpml)],
    resolution=resolution,
    dimensions=3
)

dft_main_xz = sim_main.add_dft_fields(
    [mp.Ex, mp.Ey], fcen, fcen, 1,
    center=mp.Vector3(),
    size=mp.Vector3(sx_total, 0, sz_total)
)

sim_main.run(until=t_end)

ex_raw = sim_main.get_dft_array(dft_main_xz, mp.Ex, 0)
ey_raw = sim_main.get_dft_array(dft_main_xz, mp.Ey, 0)
E2_main = np.abs(ex_raw)**2 + np.abs(ey_raw)**2

E2_norm = E2_main / E2_ref_val  

nx, nz = E2_norm.shape  
x_coords = np.linspace(-sx_total / 2, sx_total / 2, nx) 
z_coords = np.linspace(-sz_total / 2, sz_total / 2, nz) - z_shadow_meep  

z_mask = z_coords > 0.1  
maximum_behind = E2_norm[:, z_mask]
max_idx = np.unravel_index(np.argmax(maximum_behind), maximum_behind.shape)
z_max_val = z_coords[z_mask][max_idx[1]]
x_max_val = x_coords[max_idx[0]]

maximum = np.max(maximum_behind)

#print(f"max: {np.max(E2_main):.2f}")
#print(f"Reference: {E2_ref_val:.2f}")
#print(f"Max: {maximum:.2f}")
#print(f"Zmax: {z_max_val:.2f} mm")

plt.figure(figsize=(11, 6))

im = plt.pcolormesh(z_coords, x_coords, E2_norm, cmap='jet', shading='auto', vmin=0, vmax=maximum)

ax = plt.gca()
ax.set_aspect('equal', adjustable='box')

rect_z = [z_cube_center - L_cube / 2 - z_shadow_meep,
          z_cube_center + L_cube / 2 - z_shadow_meep,
          z_cube_center + L_cube / 2 - z_shadow_meep,
          z_cube_center - L_cube / 2 - z_shadow_meep,
          z_cube_center - L_cube / 2 - z_shadow_meep]

h = L_cube / 2

rect_x = [-h, -h, h, h, -h]
plt.plot(rect_z, rect_x, color='cyan', lw=2, alpha=0.9)  
plt.fill(rect_z, rect_x, color='cyan', alpha=0.15)         

plt.scatter([z_max_val], [x_max_val], color='white', s=40, edgecolors='black', zorder=5)
plt.text(0.98, 0.98, f'z_max = {z_max_val:.2f} mm',
         transform=ax.transAxes,
         color='black', fontweight='bold', fontsize=11,
         ha='right', va='top',
         bbox=dict(facecolor='white', edgecolor='gray', boxstyle='round,pad=0.3'))

plt.xlim(-7, 8)
plt.ylim(-4, 4)
cb = plt.colorbar(im)
cb.set_label("Wzmocnienie Intensywności")
plt.xlabel('Oś propagacji Z (mm)')
plt.ylabel('Oś poprzeczna X (mm)')

plt.tight_layout()
plt.savefig('terajet_3d.png', dpi=200)