Chiral optical structures, characterized by their absence of mirror symmetry, interact distinctively with left circularly polarized (LCP) and right circularly polarized (RCP) light, resulting in circular dichroism (CD) responses. This phenomenon is crucial in diverse research fields, such as drug synthesis, chiral sensing, and optical communication. However, natural optical materials exhibit very weak CD due to limited chiral light-matter interaction strength. In contrast, optical metasurfaces, as artificially engineered materials, can exhibit electromagnetic properties surpassing those of natural materials. Specifically, optical metasurfaces with chiral structures as unit cells can demonstrate strong chiroptical properties with high CD.

In this notebook, we demonstrate a wavelength-tunable infrared chiral metasurface, incorporating the phase-change material GST-225, with a parallelogram-shaped resonator as its unit cell. The metasurface's resonance wavelength is tuned through the phase transition of GST-225. The design is based on Haotian Tang, Liliana Stan, David A. Czaplewski, Xiaodong Yang, and Jie Gao, "Wavelength-tunable infrared chiral metasurfaces with phase-change materials," Opt. Express 31, 21118-21127 (2023) DOI:10.1364/OE.489841.

For more metamaterial and metasurface examples, please see the dielectric metasurface absorber, the graphene metamaterial absorber, the Fano metasurface, and the gradient metasurface reflector.

[1]:

import numpy as np
import matplotlib.pyplot as plt

import tidy3d as td
import tidy3d.web as web
from tidy3d.plugins.dispersion import FastDispersionFitter, AdvancedFastFitterParam

Simulation Setup¶

In this example, we investigate the wavelength range of 2-3 µm.

[2]:

lda0 = 2.5  # central wavelength
freq0 = td.C_0 / lda0  # central frequency
ldas = np.linspace(2, 3, 101)  # wavelength range
freqs = td.C_0 / ldas  # frequency range
fwidth = 0.5 * (np.max(freqs) - np.min(freqs))  # width of the source frequency range

Define Geometric Parameters¶

The physical meaning of the parameters is schematically shown in Fig. 1(b) of the paper.

[3]:

Px = 1.25  # unit cell period in the x direction
Py = 0.4  # unit cell period in the y direction
A = 0.24  # width of the Al resonator
B = 0.21  # width of the slot
alpha = 30 * np.pi / 180  # tilted angle of the slot
t_gst = 0.075  # thickness of the GST layer
t_al2o3 = 0.025  # thickness of the Al2O3 capping layer
t_al = 0.06  # thickness of the Al layer

Define Material Media¶

For aluminum, we directly use the predefined medium from the material library. Alumina is defined by a Sellmeier model.

[4]:

Al = td.material_library["Al"]["RakicLorentzDrude1998"]

coeffs = [(1.431, 0.0727**2), (0.651, 0.119**2), (5.341, 18.028**2)]
Al2O3 = td.Sellmeier(coeffs=coeffs)

The crystalline and amorphous states of GST are defined by using the FastDispersionFitter to fit wavelength-dependent refractive index profiles stored in CSV files.

[5]:

fitter = FastDispersionFitter.from_file("misc/GST_A_nk.csv", skiprows=0, delimiter=",")
advanced_param = AdvancedFastFitterParam(weights=(1, 1))
GST_A, _ = fitter.fit(max_num_poles=4, advanced_param=advanced_param, tolerance_rms=0.05)
fitter.plot(GST_A)
plt.show()

Output()

16:46:29 Eastern Standard Time WARNING: Unable to fit with weighted RMS error   
                               under 'tolerance_rms' of 0.05

No description has been provided for this image

[6]:

fitter = FastDispersionFitter.from_file("misc/GST_C_nk.csv", skiprows=0, delimiter=",")
GST_C, _ = fitter.fit(max_num_poles=4, advanced_param=advanced_param, tolerance_rms=0.05)
fitter.plot(GST_C)
plt.show()

Output()

Define Structures in the Unit Cell¶

The unit cell includes a 200 nm-thick aluminum mirror layer on a silicon substrate, followed by a 75 nm GST-225 spacer layer. Above this, there is a 25 nm-thick alumina capping layer, and the assembly is topped with a 60 nm-thick aluminum layer, patterned into parallelogram-shaped resonators. The resonators are defined as PolySlabs with calculated vertices. Since we need to simulate the metasurface with the GST in two states, we define a make_gst_layer(state) function to construct the GST layer given the state.

[7]:

inf_eff = 1e2  # effective infinity

# calculate the vertices of the resonator on the left
vertices = [
    (-inf_eff, A / 2),
    (-inf_eff, -A / 2),
    (-(B / np.cos(alpha) + A * np.tan(alpha)) / 2, -A / 2),
    (-(B / np.cos(alpha) - A * np.tan(alpha)) / 2, A / 2),
]

# define the left resonator structure
left_resonator = td.Structure(
    geometry=td.PolySlab(vertices=vertices, axis=2, slab_bounds=(0, t_al)), medium=Al
)

# calculate the vertices of the resonator on the right
vertices = [
    (inf_eff, A / 2),
    (inf_eff, -A / 2),
    ((B / np.cos(alpha) - A * np.tan(alpha)) / 2, -A / 2),
    ((B / np.cos(alpha) + A * np.tan(alpha)) / 2, A / 2),
]

# define the right resonator structure
right_resonator = td.Structure(
    geometry=td.PolySlab(vertices=vertices, axis=2, slab_bounds=(0, t_al)), medium=Al
)

# define the Al2O3 buffer layer
buffer_layer = td.Structure(
    geometry=td.Box(center=(0, 0, -t_al2o3 / 2), size=(td.inf, td.inf, t_al2o3)), medium=Al2O3
)

# define a function to construct the GST layer given the state
def make_gst_layer(state):
    if state == "amorphous":
        medium = GST_A
    elif state == "crystalline":
        medium = GST_C
    else:
        raise ValueError("state must be `amorphous` or `crystalline`")

    gst_layer = td.Structure(
        geometry=td.Box(center=(0, 0, -t_al2o3 - t_gst / 2), size=(td.inf, td.inf, t_gst)),
        medium=medium,
    )
    return gst_layer


# define the Al substrate
al_substrate = td.Structure(
    geometry=td.Box.from_bounds(
        rmin=(-inf_eff, -inf_eff, -inf_eff), rmax=(inf_eff, inf_eff, -t_al2o3 - t_gst)
    ),
    medium=Al,
)

Define Source and Monitors¶

A circularly polarized plane wave is just two perpendicularly linearly polarized plane waves with a 90 degree phase difference. Here we define a circular_polarized_plane_wave(pol) function to construct the circularly polarized plane wave given the polarization, either left or right.

Since the metasurface is fabricated on a metal mirror, the transmission will be zero. We only need to define a FluxMonitor to measure reflection ($R$). Absorption can be calculated simply by $A=1-R$. To further inspect the resonant field patterns, we define two FieldMonitors in the $xz$ and $xy$ planes. To reduce the result data size, we only record fields at the resonant wavelength when GST is in the amorphous state, which is predetermined to be 2.4 µm.

[8]:

def circular_polarized_plane_wave(pol):
    # define a plane wave polarized in the x direction
    plane_wave_x = td.PlaneWave(
        source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
        size=(td.inf, td.inf, 0),
        center=(0, 0, 0.3 * lda0),
        direction="-",
        pol_angle=0,
    )

    # determine the phase difference given the polarization
    if pol == "left":
        phase = -np.pi / 2
    elif pol == "right":
        phase = np.pi / 2
    else:
        raise ValueError("pol must be `left` or `right`")

    # define a plane wave polarized in the y direction with a phase difference
    plane_wave_y = td.PlaneWave(
        source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth, phase=phase),
        size=(td.inf, td.inf, 0),
        center=(0, 0, 0.3 * lda0),
        direction="-",
        pol_angle=np.pi / 2,
    )

    return [plane_wave_x, plane_wave_y]


# add a flux monitor to detect transmission
monitor_r = td.FluxMonitor(
    center=[0, 0, 0.5 * lda0], size=[td.inf, td.inf, 0], freqs=freqs, name="R"
)

freq_res = td.C_0 / 2.4  # resonant frequency

# add field monitors to see the field profiles at the absorption peak frequency when GST is in the amorphous state
monitor_field_xz = td.FieldMonitor(
    center=(0, 0, 0), size=(td.inf, 0, td.inf), freqs=[freq_res], name="field_xz"
)

monitor_field_xy = td.FieldMonitor(
    center=(0, 0, 0), size=(td.inf, td.inf, 0), freqs=[freq_res], name="field_xy"
)

Define Simulation¶

Now we are ready to set up the simulation. In total we need to run four simulations to cover both polarizations and GST states. To make it more convenient, we define a make_sim(pol, state) function.

[9]:

run_time = 3e-13  # simulation run time

# simulation domain box
sim_box = td.Box.from_bounds(
    rmin=(-Px / 2, -Py / 2, -t_al2o3 - t_gst - lda0 / 2), rmax=(Px / 2, Py / 2, t_al + lda0 / 2)
)

# define a function to construct a simulation given the polarization and GST state
def make_sim(pol, state):
    sim = td.Simulation(
        center=sim_box.center,
        size=sim_box.size,
        grid_spec=td.GridSpec.auto(min_steps_per_wvl=20, wavelength=lda0),
        structures=[
            left_resonator,
            right_resonator,
            buffer_layer,
            make_gst_layer(state),
            al_substrate,
        ],
        sources=circular_polarized_plane_wave(pol),
        monitors=[monitor_r, monitor_field_xz, monitor_field_xy],
        run_time=run_time,
        boundary_spec=td.BoundarySpec(
            x=td.Boundary.periodic(), y=td.Boundary.periodic(), z=td.Boundary.pml()
        ),
    )
    return sim


# make one simulation and plot it to inspect
sim_left_A = make_sim("left", "amorphous")
sim_left_A.plot_3d()

In addition, we can inspect the grid to ensure it is sufficiently fine.

[10]:

ax = sim_left_A.plot(y=0)
sim_left_A.plot_grid(y=0, ax=ax)
ax.set_ylim(-0.1, 0.1)
ax.set_ylim(-0.2, 0.2)
ax.set_aspect("auto")

Run Simulation Batch¶

To make the simulations more efficient, we define a Batch and submit it to run all four simulations in parallel.

[11]:

# define four simulations
sims = {
    "LCP_A": sim_left_A,
    "LCP_C": make_sim("left", "crystalline"),
    "RCP_A": make_sim("right", "amorphous"),
    "RCP_C": make_sim("right", "crystalline"),
}

# create a batch and run all sims in parallel
batch = web.Batch(simulations=sims)
batch_results = batch.run(path_dir="data")

16:46:32 Eastern Standard Time Created task 'LCP_A' with task_id                
                               'fdve-98d8bf11-9921-427a-8a87-20f9327406ea' and  
                               task_type 'FDTD'.

                               View task using web UI at                        
                               'https://tidy3d.simulation.cloud/workbench?taskId
                               =fdve-98d8bf11-9921-427a-8a87-20f9327406ea'.

Output()

16:46:33 Eastern Standard Time Created task 'LCP_C' with task_id                
                               'fdve-75ab7a60-84bd-47a9-ad50-fdf455195bbc' and  
                               task_type 'FDTD'.

                               View task using web UI at                        
                               'https://tidy3d.simulation.cloud/workbench?taskId
                               =fdve-75ab7a60-84bd-47a9-ad50-fdf455195bbc'.

Output()

                               Created task 'RCP_A' with task_id                
                               'fdve-7b63bc41-38b1-4661-af82-485ce0b9db93' and  
                               task_type 'FDTD'.

                               View task using web UI at                        
                               'https://tidy3d.simulation.cloud/workbench?taskId
                               =fdve-7b63bc41-38b1-4661-af82-485ce0b9db93'.

Output()

16:46:34 Eastern Standard Time Created task 'RCP_C' with task_id                
                               'fdve-90760f28-e71c-44d6-9af1-c9dd9d30b7a6' and  
                               task_type 'FDTD'.

                               View task using web UI at                        
                               'https://tidy3d.simulation.cloud/workbench?taskId
                               =fdve-90760f28-e71c-44d6-9af1-c9dd9d30b7a6'.

Output()

16:46:38 Eastern Standard Time Started working on Batch.

16:47:09 Eastern Standard Time Maximum FlexCredit cost: 1.194 for the whole     
                               batch.

                               Use 'Batch.real_cost()' to get the billed        
                               FlexCredit cost after the Batch has completed.

Output()

16:52:03 Eastern Standard Time Batch complete.

Result Analysis¶

After the simulations are complete, we can extract the absorption spectra for all simulations.

[12]:

# extract absorption from four simulations
A_LCP_A = 1 - batch_results["LCP_A"]["R"].flux / 2
A_LCP_C = 1 - batch_results["LCP_C"]["R"].flux / 2
A_RCP_A = 1 - batch_results["RCP_A"]["R"].flux / 2
A_RCP_C = 1 - batch_results["RCP_C"]["R"].flux / 2

# plot absorption spectra
plt.plot(ldas, A_LCP_A, "red", label="LCP_A")
plt.plot(ldas, A_RCP_A, "red", linestyle="--", label="RCP_A")
plt.plot(ldas, A_LCP_C, "blue", label="LCP_C")
plt.plot(ldas, A_RCP_C, "blue", linestyle="--", label="RCP_C")
plt.ylabel("Absorption")
plt.xlabel("Wavelength ($\mu m$)")
plt.legend()
plt.show()

Output()

16:52:06 Eastern Standard Time loading simulation from                          
                               data\fdve-98d8bf11-9921-427a-8a87-20f9327406ea.hd
                               f5

Output()

16:52:08 Eastern Standard Time loading simulation from                          
                               data\fdve-75ab7a60-84bd-47a9-ad50-fdf455195bbc.hd
                               f5

Output()

16:52:09 Eastern Standard Time loading simulation from                          
                               data\fdve-7b63bc41-38b1-4661-af82-485ce0b9db93.hd
                               f5

Output()

16:52:11 Eastern Standard Time loading simulation from                          
                               data\fdve-90760f28-e71c-44d6-9af1-c9dd9d30b7a6.hd
                               f5

With the absorption spectra, we can further calculate CD for both GST states. Note that compared to the publication, the results are similar but not identical, which is likely due to the slightly different refractive indices for the materials.

[13]:

# plot CD spectra
plt.plot(ldas, A_LCP_A - A_RCP_A, "red", label="A")
plt.plot(ldas, A_LCP_C - A_RCP_C, "blue", label="C")
plt.title("Circular dichroism (CD)")
plt.ylabel("CD")
plt.xlabel("Wavelength ($\mu m$)")
plt.legend()
plt.show()

Lastly, we can plot the field distributions at the resonant wavelength when GST is in the amorphous state to visualize the resonant mode. We can see that with the left circularly polarized incident, the response is much stronger, leading to a much higher absorption as seen above.

[14]:

fig, ax = plt.subplots(2, 1, tight_layout=True)
batch_results["LCP_A"].plot_field(
    field_monitor_name="field_xy", field_name="E", val="abs", vmin=0, vmax=500, ax=ax[0]
)
batch_results["RCP_A"].plot_field(
    field_monitor_name="field_xy", field_name="E", val="abs", vmin=0, vmax=500, ax=ax[1]
)
plt.show()

16:52:12 Eastern Standard Time loading simulation from                          
                               data\fdve-98d8bf11-9921-427a-8a87-20f9327406ea.hd
                               f5

                               loading simulation from                          
                               data\fdve-7b63bc41-38b1-4661-af82-485ce0b9db93.hd
                               f5

Similarly, plot the field distributions in the $xz$ plane.

[15]:

fig, ax = plt.subplots(1, 2, figsize=(8, 3), tight_layout=True)
batch_results["LCP_A"].plot_field(
    field_monitor_name="field_xz", field_name="E", val="abs", vmin=0, vmax=500, ax=ax[0]
)
ax[0].set_ylim(-0.2, 1)
batch_results["RCP_A"].plot_field(
    field_monitor_name="field_xz", field_name="E", val="abs", vmin=0, vmax=500, ax=ax[1]
)
ax[1].set_ylim(-0.2, 1)
plt.show()

16:52:13 Eastern Standard Time loading simulation from                          
                               data\fdve-98d8bf11-9921-427a-8a87-20f9327406ea.hd
                               f5

                               loading simulation from                          
                               data\fdve-7b63bc41-38b1-4661-af82-485ce0b9db93.hd
                               f5