Unidirectional SPP from non-Hermitian metagratings

Unidirectional SPP from non-Hermitian metagratings#

Non-Hermitian systems are open systems that can exchange energy, matter, or information with their surroundings through the engineering of gain or loss materials. While quantum mechanics traditionally requires Hamiltonians to be Hermitian to produce real eigenvalues, it’s been discovered that a system with a spatially varying complex potential energy can still have real eigenvalues if it obeys a certain symmetry known as Parity-Time (PT) symmetry.

One concept that has gained particular attention in non-Hermitian photonics is the phenomenon of β€œexceptional points.” These are points in parameter space at which the eigenvalues and the eigenvectors of a non-Hermitian system coincide. Around these exceptional points, systems can have unique and often counterintuitive properties. For instance, light in these systems can exhibit uni-directional propagation, enhanced sensitivity to external parameters, or unusual scattering behavior. Non-Hermitian photonics has potential applications in many areas, including the design of new types of lasers, sensors, and photonic devices with novel properties.

This notebook demonstrates the unidirectional surface plasmon polariton (SPP) propagation from non-Hermitian metagratings. The design is adapted from Yihao Xu et al. ,Subwavelength control of light transport at the exceptional point by non-Hermitian metagratings. Sci. Adv. 9, eadf3510(2023).

Schematic of the non-Hermitian metagrating

For more interesting contents in the realm of nanophotonics, please check out the case studies on Anderson localization, hyperbolic polaritons, and plasmonic nanoantennas. For more simulation examples, please visit our examples page. If you are new to the finite-difference time-domain (FDTD) method, we highly recommend going through our FDTD101 tutorials. FDTD simulations can diverge due to various reasons. If you run into any simulation divergence issues, please follow the steps outlined in our troubleshooting guide to resolve it.

[1]:
import numpy as np
import matplotlib.pyplot as plt

import tidy3d as td
import tidy3d.web as web
from tidy3d.plugins.dispersion import DispersionFitter

Simulation Setup#

Define the basic model parameters including the simulation wavelength/frequency and the metagrating geometric parameters.

[2]:
lda0 = 1.15 # simulation wavelength
freq0 = td.C_0/lda0 # simulation frequency

w1 = 0.083 # width of the first nanostrip in the unit cell
w2 = 0.162 # width of the second nanostrip in the unit cell
t1 = 0.08  # thickness of the first nanostrip
t2 = 0.03  # thickness of the bottom layer of the second nanostrip
t3 = 0.022 # thickness of the top layer of the second nanostrip
s = 0.35  # center-to-center distance between the nanostrips
p = 1.14  # grating period
t_au = 0.06  # gold layer thickness
N = 9  # number of unit cells
inf_eff = 1e3 # effective infinity

The metagrating utilizes four materials: a silicon dioxide substrate, a gold thin film, germanium nanostrips, and chromium nanostrips. These materials can be defined in different ways. Specifically, we use the gold medium directly from the material library. The germanium is defined using the from_nk method by specifying the real and imaginary parts of the refractive index at the simulation frequency. Silicon dioxide is considered as a lossless medium. Chromium is modeled as a dispersive medium by fitting the raw data from experimental measurements.

In this simulation, we are only concerned about one frequency, so in principle the dispersive fitting is not necessary. However, chromium has \(Re(\varepsilon)<1\) at 1150 nm. In Tidy3D, we require the nondispersive medium to have \(Re(\varepsilon)\ge 1\). Therefore, it can not be defined using the from_nk method directly and dispersive fitting is the only option.

The raw data file Cr_Sytchkova.csv can be found on the refractive index database or on our github repo.

[3]:
# define gold medium
Au = td.material_library['Au']['Olmon2012stripped']

# define Ge medium
Ge = td.Medium.from_nk(n=4.6403, k=0.2519, freq=freq0)

# define SiO2 medium
SiO2 = td.Medium(permittivity=1.45**2)

# define Cr medium
fname = "misc/Cr_Sytchkova.csv"  # read the refractive index data from a csv file
fitter = DispersionFitter.from_file(fname, delimiter=",")  # construct a fitter
Cr, rms_error = fitter.fit(num_poles=4, num_tries=10)
fitted_nk = Cr.nk_model(freq0)

[14:47:33] WARNING: Unable to fit with RMS error under                fit.py:321
           'tolerance_rms' of 0.01                                              

After the fitting, we see a warning that the fitting RMS error didn’t get lower than the threshold. We can inspect the fitting result of chromium to see if the fitting is sufficiently good. Here we do see the fitting is quite resonable.

[4]:
fitter.plot(Cr)
plt.show()
../_images/notebooks_NonHermitianMetagratings_9_0.png

Next, we define the metagrating structures. The N unit cells can be defined by using a simple for loop.

[5]:
# define the nanostrips that make up the unit cells
grating = []
for i in range(N):
    geo = td.Box(center=(-s/2+i*p-(N-1)*p/2, 0, t1/2), size=(w1,td.inf,t1))
    grating.append(td.Structure(geometry=geo, medium=Ge))

    geo = td.Box(center=(s/2+i*p-(N-1)*p/2, 0, t2/2), size=(w1,td.inf,t2))
    grating.append(td.Structure(geometry=geo, medium=Ge))

    geo = td.Box(center=(s/2+i*p-(N-1)*p/2, 0, t2+t3/2), size=(w1,td.inf,t3))
    grating.append(td.Structure(geometry=geo, medium=Cr))

# define the substrate
substrate_geo = td.Box.from_bounds(rmin=(-inf_eff,-inf_eff,-inf_eff), rmax=(inf_eff,inf_eff,-t_au))
substrate = td.Structure(geometry=substrate_geo, medium=SiO2)

# define the gold film
gold_film_geo = td.Box(center=(0,0,-t_au/2), size=(inf_eff, inf_eff, t_au))
gold_film = td.Structure(geometry=gold_film_geo, medium=Au)

For excitation, we define a GaussianBeam incident on the metagrating. The waist_radius is set to five grating periods to ensure the entire grating array is covered by the beam.

To visualize the SPP propagation, a FieldMonitor is defined. Two FluxMonitors are defined to measure the power transmissions of the SPP to the left and right.

[6]:
gaussian_beam = td.GaussianBeam(
    size=(td.inf, td.inf, 0),
    center=(0, 0, 5*lda0),
    source_time=td.GaussianPulse(freq0=freq0, fwidth=freq0/10),
    direction="-",
    waist_radius=5*p,
)

field_monitor = td.FieldMonitor(
    size=(td.inf, 0, td.inf),
    freqs=[freq0],
    interval_space=(2,1,2),
    name="field"
)

flux_monitor_right = td.FluxMonitor(
    center=(23*p,0,0),
    size=(0, td.inf, 2*lda0),
    freqs=[freq0],
    name="flux_right",
)

flux_monitor_left = td.FluxMonitor(
    center=(-23*p,0,0),
    size=(0, td.inf, 2*lda0),
    freqs=[freq0],
    name="flux_left",
)
[14:47:34] WARNING: Default value for the field monitor           monitor.py:261
           'colocate' setting has changed to 'True' in Tidy3D                   
           2.4.0. All field components will be colocated to the                 
           grid boundaries. Set to 'False' to get the raw fields                
           on the Yee grid instead.                                             

Define the Tidy3D Simulation. The simulation is 2D so periodic boundary condition is applied to the direction with zero size (the \(y\) direction in this case).

[7]:
# simulation domain size
Lx = 50*p
Ly = 0
Lz = 7*lda0

# simulation run time
run_time = 5e-13

# define simulation
sim = td.Simulation(
        center=(0, 0, 2*lda0),
        size=(Lx, Ly, Lz),
        grid_spec=td.GridSpec.auto(min_steps_per_wvl=30, wavelength=lda0),
        structures=grating + [substrate, gold_film],
        sources=[gaussian_beam],
        monitors=[field_monitor, flux_monitor_right, flux_monitor_left],
        run_time=run_time,
        boundary_spec=td.BoundarySpec(
            x=td.Boundary.pml(), y=td.Boundary.periodic(), z=td.Boundary.pml()
        ),
)

Visualize the simulation setup.

[8]:
ax = sim.plot(y=0)
ax.set_aspect('auto')
plt.show()
../_images/notebooks_NonHermitianMetagratings_17_0.png

Since the grating structures are small compared to the simulation domain, we can zoom in to see them more clearly.

[9]:
ax = sim.plot(y=0)
ax.set_xlim(-6, 6)
ax.set_ylim(-0.1, 0.1)
ax.set_aspect('auto')
plt.show()
../_images/notebooks_NonHermitianMetagratings_19_0.png

The simulation settings look correct. We are ready to submit the simulation to the server.

[10]:
sim_data = web.run(
    sim, task_name="nonhermitian_metagrating", path="data/simulation_data.hdf5"
)

[14:47:35] Created task 'nonhermitian_metagrating' with task_id    webapi.py:188
           'fdve-900aa944-73ee-4b0f-8481-25b29f57fcd7v1'.                       
[14:47:36] status = queued                                         webapi.py:361
[14:47:44] status = preprocess                                     webapi.py:355
[14:47:49] Maximum FlexCredit cost: 0.059. Use                     webapi.py:341
           'web.real_cost(task_id)' to get the billed FlexCredit                
           cost after a simulation run.                                         
           starting up solver                                      webapi.py:377
           running solver                                          webapi.py:386
           To cancel the simulation, use 'web.abort(task_id)' or   webapi.py:387
           'web.delete(task_id)' or abort/delete the task in the                
           web UI. Terminating the Python script will not stop the              
           job running on the cloud.                                            
[14:48:08] early shutoff detected, exiting.                        webapi.py:404
           status = postprocess                                    webapi.py:419
[14:48:19] status = success                                        webapi.py:426
[14:48:21] loading SimulationData from data/simulation_data.hdf5   webapi.py:590

Result Analysis#

Plot \(|H_y|^2\) to visualize the SPP propagation. We can see the SPP is only propagating to the right, i.e. unidirectional propagation.

[11]:
ax = sim_data.plot_field(field_monitor_name='field', field_name='Hy', val='abs^2', vmin=0, vmax=3e-4)
ax.set_aspect('auto')
../_images/notebooks_NonHermitianMetagratings_24_0.png

To quantify the directionality, we define a contrast parameter \(C_{exc} = \frac{I_r-I_l}{I_r+I_l}\), where \(I_r\) and \(I_l\) are the transmitted power to the right and to the left, respectively. \(C_{exc}=0\) corresponds to an isotropic propagation. For the ideal unidirectional propagation, \(C_{exc}=1\) or \(C_{exc}=-1\). In this case, we observe \(C_{exc}>0.98\).

[12]:
def cal_C_exc(sim_data):
    I_r = sim_data['flux_right'].flux
    I_l = -sim_data['flux_left'].flux
    C_exc = (I_r-I_l)/(I_r+I_l)
    return C_exc.values[0]

print(f'C_exc = {cal_C_exc(sim_data):1.4f}')
C_exc = 0.9818

Sweeping Nanostrip Spacing#

The SPP propagating can be effectively tuned by changing the center-to-center distance between the nanostrips within the grating unit cell. We would like to perform a parameter sweep of the center-to-center distance \(s\) to investigate how the unidirectional SPP propagation changes. To do so, we first define a function make_sim that returns a simulation given a \(s\) value.

[13]:
def make_sim(s):

    # define the grating structures of each unit cell
    grating = []
    for i in range(N):
        geo = td.Box(center=(-s/2+i*p-(N-1)*p/2, 0, t1/2), size=(w1,td.inf,t1))
        grating.append(td.Structure(geometry=geo, medium=Ge))

        geo = td.Box(center=(s/2+i*p-(N-1)*p/2, 0, t2/2), size=(w1,td.inf,t2))
        grating.append(td.Structure(geometry=geo, medium=Ge))

        geo = td.Box(center=(s/2+i*p-(N-1)*p/2, 0, t2+t3/2), size=(w1,td.inf,t3))
        grating.append(td.Structure(geometry=geo, medium=Cr))

    # define simulation
    sim = td.Simulation(
        center=(0, 0, 2*lda0),
        size=(Lx, Ly, Lz),
        grid_spec=td.GridSpec.auto(min_steps_per_wvl=30, wavelength=lda0),
        structures=grating + [substrate, gold_film],
        sources=[gaussian_beam],
        monitors=[field_monitor, flux_monitor_right, flux_monitor_left],
        run_time=run_time,
        boundary_spec=td.BoundarySpec(
            x=td.Boundary.pml(), y=td.Boundary.periodic(), z=td.Boundary.pml()
        ),
    )

    return sim

The range of interest for \(s\) is from 350 nm to 770 nm. A batch of 5 simulations is created.

[14]:
s_array = np.linspace(0.35,0.77,5)  # collection of s for the parameter sweep

sims = {
    f"s={s:.3f}": make_sim(s) for s in s_array
}  # construct simulations for each s from s_array

Run the simulations in the batch in parallel.

[15]:
batch = web.Batch(simulations=sims)
batch_results = batch.run(path_dir="data")
[14:48:27] Created task 's=0.350' with task_id                     webapi.py:188
           'fdve-eb6d2dd3-b811-4bd0-9e62-2b8ac960291av1'.                       
[14:48:28] Created task 's=0.455' with task_id                     webapi.py:188
           'fdve-401d585f-970b-446d-99e0-39a3edabc191v1'.                       
           Created task 's=0.560' with task_id                     webapi.py:188
           'fdve-703b1ff3-a7c3-4e8c-b560-b684d164ea16v1'.                       
[14:48:29] Created task 's=0.665' with task_id                     webapi.py:188
           'fdve-e6d533b6-9fe3-4726-b9d8-5fa7cf198908v1'.                       
           Created task 's=0.770' with task_id                     webapi.py:188
           'fdve-dc0346dc-e4ac-41ce-b2d1-acfefeb7ce65v1'.                       
[14:48:31] Started working on Batch.                            container.py:475
[14:49:08] Maximum FlexCredit cost: 0.297 for the whole batch.  container.py:479
           Use 'Batch.real_cost()' to get the billed FlexCredit                 
           cost after the Batch has completed.                                  
[14:50:09] Batch complete.                                      container.py:522

Plot \(|H_y|^2\) for each simulation. We can see that the SPP propagates to the right for smaller \(s\) and to the left for larger \(s\).

\(C_{exc}\) is also computed for each simulation.

[16]:
fig, ax = plt.subplots(len(s_array),1, tight_layout=True, figsize=(8,12))

C_exc = []

# computer |H_y|^2 and C_exc for each simulation
for i, s in enumerate(s_array):
    sim_data = batch_results[f"s={s:.3f}"]
    Hy_abs = sim_data['field'].Hy.abs
    Hy_squared = Hy_abs**2
    Hy_squared.plot(x='x',y='z',ax=ax[i], vmin=0, vmax=3e-4, cmap='hot')
    ax[i].set_title(f's={s:.3f} um')
    C_exc.append(cal_C_exc(sim_data))
[14:50:12] loading SimulationData from                             webapi.py:590
           data/fdve-eb6d2dd3-b811-4bd0-9e62-2b8ac960291av1.hdf5                
[14:50:14] loading SimulationData from                             webapi.py:590
           data/fdve-401d585f-970b-446d-99e0-39a3edabc191v1.hdf5                
[14:50:16] loading SimulationData from                             webapi.py:590
           data/fdve-703b1ff3-a7c3-4e8c-b560-b684d164ea16v1.hdf5                
[14:50:18] loading SimulationData from                             webapi.py:590
           data/fdve-e6d533b6-9fe3-4726-b9d8-5fa7cf198908v1.hdf5                
[14:50:19] loading SimulationData from                             webapi.py:590
           data/fdve-dc0346dc-e4ac-41ce-b2d1-acfefeb7ce65v1.hdf5                
../_images/notebooks_NonHermitianMetagratings_35_20.png

Plot \(C_{exc}\) as a function of \(s\). \(C_{exc}\) can be tuned from >0.98 to <-0.98 by changing \(s\). Compared to the original paper, the result here is slightly different, which we attribute to slightly different material properties used.

[17]:
plt.plot(s_array*1e3, np.array(C_exc), c='red', linewidth=3)
plt.xlabel('s (nm)')
plt.ylabel('$C_{exc}$')
plt.ylim(-1,1)
plt.show()
../_images/notebooks_NonHermitianMetagratings_37_0.png
[ ]: