TIDY3D
LEARNING CENTER

Lectures
Tutorials
Applications

Example Library / Mediums / Defining spatially-varying dielectric structures

Defining spatially-varying dielectric structures

Sign up now to simulate .ipynb

When a regular Tidy3D Medium or a dispersive medium such as Lorentz is assigned to a Structure, the refractive index is homogeneous within the geometry. In some cases, such as gradient-index optics component and inhomogeneous doped structure simulations, a spatially varying refractive index distribution is desirable. In principle, this can be achieved by manually dividing the Structure into smaller sub-components and assigning a refractive index value to each component according to the spatial distribution. However, this process can be tedious and error-prone. Fortunately, Tidy3D supports both non-dispersive and dispersive media that have a spatially varying refractive index, and thus can help you achieve this result very conveniently.

In this tutorial, first we illustrate how to define structures made of custom non-dispersive medium called CustomMedium that allow for a customized refractive index spatial profile. An example is an ideal gradient index lens that has flat surfaces but a parabolic distribution of refractive index. The simulation result allows us to examine the focusing capability of the lens.

Next, we illustrate how to define structures made of custom dispersive medium. As an example, we construct a lens made of CustomLorentz that has the ideal refractive index distribution at the frequency of interest.

Finally, we explain how to construct a diagonally anisotropic medium in which each component can be independently defined as any type of custom dispersive or non-dispersive medium.

Schematic of the graded index lens

If you are new to the finite-difference time-domain (FDTD) method, we highly recommend going through our FDTD101 tutorials.

import matplotlib.pyplot as plt
import numpy as np
import tidy3d as td
import xarray as xr
from tidy3d import SpatialDataArray, web

The lens is designed to work at 1 \(\mu m\) wavelength.

lda0 = 1  # central wavelength

freq0 = td.C_0 / lda0  # central frequency

Custom Non-dispersive Medium

With CustomMedium, one can customize the spatial profile of non-dispersive refractive index inside a structure. The spatial profile is defined with SpatialDataArray where we provide scalar data such as refractive index and labeled coordinates. Below, let’s see how to define a flat lens whose index of refraction follows \(n(r)=n_0(1-Ar^2)\), where \(r\) is the radial distance to the \(z\)-axis. In this particular example, we study a lens of 20 \(\mu m\) radius and 10 \(\mu m\) thickness.

First, we define the spatial grids. The refractive index is invariant along the \(z\)-axis, and varies in the \(x\)-\(y\) plane. For the uniform \(z\)-axis, we only need to supply one grid point since CustomMedium will automatically generate uniform profiles along the axis where a single grid point is supplied. In the \(x\) and \(y\) dimensions, we set 100 grid points to fully resolve the refractive index variation.

Nx, Ny, Nz = 100, 100, 1  # number of grid points along each dimension

r = 20  # radius of the lens, unit: micron
t = 10  # thickness of the lens, unit: micron

# The coordinate for the refractive index data that includes x, y, z, and frequency
# Note: when only one coordinate is supplied along an axis, it means the medium is uniform along this axis.
X = np.linspace(-r, r, Nx)  # x grid
Y = np.linspace(-r, r, Ny)  # y grid
Z = [0]  # z grid

Next, we define a 3-dimensional array that stores the refractive index.

# define coordinate array
x_mesh, y_mesh, z_mesh = np.meshgrid(X, Y, Z, indexing="ij")
r_mesh = np.sqrt(x_mesh**2 + y_mesh**2)  # radial distance

# index of refraction array
# assign the refractive index value to the array according to the desired profile
n_data = np.ones((Nx, Ny, Nz))
n0 = 2
A = 1e-3
n_data[r_mesh <= r] = n0 * (1 - A * r_mesh[r_mesh <= r] ** 2)

Finally, we convert the numpy array to a SpatialDataArray that labels the coordinate.

# convert to dataset array
n_dataset = SpatialDataArray(n_data, coords=dict(x=X, y=Y, z=Z))

Defining CustomMedium in 3 Ways

Here, we will illustrate defining custom medium in three different ways:

  1. Using classmethod td.CustomMedium.from_nk when refractive index and extinction coefficients are readily available.
  2. Using classmethod td.CustomMedium.from_eps_raw to supply permittivity data, which can be complex-valued for lossy medium.
  3. Define permittivity and conductivity directly in td.CustomMedium.

(One can also define permittivity for each component separately via PermittivityDataset, but it is to be deprecated in the future. For custom medium that is anisotropic, it is encouraged to define through CustomAnisotropicMedium)

In any of the three ways, you can optionally define how the permittivity is interpolated to Yee-grid with interp_method that can take the value of “linear” or “nearest”. The default value is “nearest”. If the custom medium is applied to a geometry larger than the custom medium’s grid range, extrapolation is automatically applied for Yee grids outside the supplied coordinate region. When the extrapolated value is smaller (greater) than the minimal (maximal) of the supplied data, the extrapolated value will take the minimal (maximal) of the supplied data.

# # Three equivalent ways of defining custom medium for the lens

# define custom medium with n/k data
mat_custom1 = td.CustomMedium.from_nk(n_dataset, interp_method="nearest")

# define custom medium with permittivity data
eps_dataset = n_dataset**2
mat_custom2 = td.CustomMedium.from_eps_raw(eps_dataset, interp_method="nearest")

# define permittivity directly in the class
mat_custom3 = td.CustomMedium(permittivity=eps_dataset, interp_method="nearest")

Note that when the medium is lossy so that k is non-zero and permittivity is complex-valued, in Approach 3, one can directly define the conductivity in the class:

mat_custom3 = td.CustomMedium(permittivity=eps_dataset, conductivity=conductivity_dataset, interp_method="nearest")

In the other two approaches, the frequency value at each the complex-valued permittivity or n/k is defined is needed to evaluate the conductivity. There are two ways to supply the frequency information:

Simulation Setup

Define Lens Structure

# define the lens structure as a box
lens = td.Structure(
    geometry=td.Box(center=(0, 0, t / 2), size=(td.inf, td.inf, t)), medium=mat_custom1
)

Define a Source and Monitor

# define a plane wave source
plane_wave = td.PlaneWave(
    source_time=td.GaussianPulse(freq0=freq0, fwidth=freq0 / 20),
    size=(td.inf, td.inf, 0),
    center=(0, 0, -lda0 / 2),
    direction="+",
    pol_angle=0,
)


# define a field monitor in the xz plane at y=0
monitor_field_xz = td.FieldMonitor(
    center=[0, 0, 0], size=[td.inf, 0, td.inf], freqs=[freq0], name="field_xz"
)

Define a Simulation

# simulation domain size
Lx, Ly, Lz = 2 * r, 2 * r, 5 * t
sim_size = (Lx, Ly, Lz)

run_time = 2e-12  # simulation run time

# define simulation
sim = td.Simulation(
    center=(0, 0, Lz / 2 - lda0),
    size=sim_size,
    grid_spec=td.GridSpec.auto(min_steps_per_wvl=10, wavelength=lda0),
    structures=[lens],
    sources=[plane_wave],
    monitors=[monitor_field_xz],
    run_time=run_time,
    boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),  # pml is applied in all boundaries
    symmetry=(
        -1,
        1,
        0,
    ),  # symmetry is used such that only a quarter of the structure needs to be modeled.
)

Visualize the Simulation and Gradient Index Distribution

We can use the plot_eps method of Simulation to visualize the simulation setup as well as the permittivity distribution.

First, plot the \(yz\) plane at \(x\)=0.

sim.plot_eps(x=0)
plt.show()

Similarly, plot the \(xy\) plane at \(z\)=0. The spatially varying permittivity is clearly observed.

sim.plot_eps(z=t / 2)
plt.show()

Submit Simulation Job

Submit the simulation job to the server.

sim_data = web.run(
    sim,
    task_name="gradient_index_lens",
    path="data/simulation.hdf5",
    verbose=True,
)
13:17:52 UTC Created task 'gradient_index_lens' with resource_id                
             'fdve-fb7e33ef-8566-4181-beb8-ccb345fd8616' and task_type 'FDTD'.  

             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-fb7e33ef-856
             6-4181-beb8-ccb345fd8616'.                                         

             Task folder: 'default'.                                            







13:18:25 UTC Estimated FlexCredit cost: 0.674. This assumes the FDTD solver runs
             for the full simulation time; if early shutoff is reached, the     
             billed cost can be lower. Use 'web.real_cost(task_id)' to get the  
             billed FlexCredit cost after a simulation run.                     







13:18:26 UTC status = queued                                                    







             To cancel the simulation, use 'web.abort(task_id)' or              
             '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.                                                             












13:18:38 UTC status = preprocess                                                







13:18:40 UTC status = queued                                                    







13:20:28 UTC status = preprocess                                                













13:20:37 UTC starting up solver                                                 







             running solver                                                     












13:21:30 UTC early shutoff detected at 0%, exiting.                             













             status = queued                                                    












13:21:49 UTC status = preprocess                                                







13:21:56 UTC status = running                                                   







13:24:53 UTC status = postprocess                                               







13:24:57 UTC status = success                                                   













13:24:59 UTC View simulation result at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-fb7e33ef-856
             6-4181-beb8-ccb345fd8616'.                                         


















13:25:01 UTC Loading results from data/simulation.hdf5                          











Result Visualization


After the simulation is complete, we can inspect the focusing capability of the gradient-index lens by plotting the field distributions. First, plot \(E_x\) in the \(xz\) plane at \(y=0\).






sim_data.plot_field("field_xz", "Ex", vmin=-15, vmax=15)
plt.show()





















The focus is better visualized by plotting the field intensity. A strong focus about 17 \(\mu m\) from the front surface of the lens is observed.






sim_data.plot_field("field_xz", "E", "abs^2", vmin=0, vmax=300)
plt.show()



























Custom Dispersive Medium


Tidy3d supports a set of dispersive media: PoleResidue, Lorentz, Sellmeier, Drude, and Debye. Each of them has been generalized to support spatially varying profile. The class for spatially varying medium has a prefix “Custom”. For example, we can define a spatially varying Lorentz model with CustomLorentz class.


The usage of the “Custom” dispersive medium class is very similar to that of the regular medium. When defining a custom dispersive medium object, we simply replace the input fields that are scalar in the regular medium with SpatialDataArray, which defines the spatial profile with labeled coordinates. Take Lorentz model as an example, there are two input fields: a scalar eps_inf that defines the relative permittivity at infinite frequency, and a nested tuple of scalar coeffs that defines the oscillator properties including oscillator strength, frequency, and damping rate. In CustomLorentz, the scalar eps_inf turns into SpatialDataArray, and the nested tuple of scalar coeffs turns into a nested tuple of SpatialDataArray. An example is illustrated below:






# spatially varying Lorentz model that has the same permittivity profile as
# the custom non-dispersive medium at `freq0`.
eps_inf_dataset = xr.ones_like(eps_dataset)  # uniform eps_inf
f0_dataset = xr.ones_like(eps_inf_dataset) * 2 * freq0  # uniform oscillator frequency as well
gamma_dataset = xr.zeros_like(eps_inf_dataset)  # lossless oscillator
delep_dataset = (eps_dataset - 1) * 3 / 4  # non-uniform oscillator strength
mat_lorentz = td.CustomLorentz(
    eps_inf=eps_inf_dataset, coeffs=((delep_dataset, f0_dataset, gamma_dataset),)
)







Note that all input fields must be defined over the same spatial grids. In this example, the data of eps_inf and all terms in coeffs are defined over the same spatial grids.




Setup Simulation with Spatially Varying Lorentz Medium






lens_lorentz = lens.copy(update={"medium": mat_lorentz})
sim_lorentz = sim.copy(update={"structures": [lens_lorentz]})











Visualize the Simulation and Gradient Index Distribution


Let’s visualize the simulation on the same plane as before. To view the permittivity distribution at freq0, we set the kwarg freq=freq0. The permittivity distribution at freq0 looks the same as the previous custom non-dispersive medium example.






sim_lorentz.plot_eps(x=0, freq=freq0)
plt.show()

























sim_lorentz.plot_eps(z=t / 2, freq=freq0)
plt.show()

























Submit Simulation Job


Submit the simulation job to the server.






sim_data_lorentz = web.run(
    sim_lorentz,
    task_name="gradient_index_lens_lorentz",
    path="data/simulation_lorentz.hdf5",
    verbose=True,
)







13:25:09 UTC Created task 'gradient_index_lens_lorentz' with resource_id        
             'fdve-95b7981b-7251-48e6-bc56-8ce8d67cbba7' and task_type 'FDTD'.  







             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-95b7981b-725
             1-48e6-bc56-8ce8d67cbba7'.                                         







             Task folder: 'default'.                                            


















13:25:41 UTC Estimated FlexCredit cost: 0.877. This assumes the FDTD solver runs
             for the full simulation time; if early shutoff is reached, the     
             billed cost can be lower. Use 'web.real_cost(task_id)' to get the  
             billed FlexCredit cost after a simulation run.                     







13:25:42 UTC status = queued                                                    







             To cancel the simulation, use 'web.abort(task_id)' or              
             '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.                                                             












13:25:52 UTC status = preprocess                                                













13:25:58 UTC starting up solver                                                 







13:25:59 UTC running solver                                                     












13:26:30 UTC early shutoff detected at 0%, exiting.                             













             status = queued                                                    












13:27:30 UTC status = preprocess                                                







13:27:35 UTC status = running                                                   







13:31:06 UTC status = postprocess                                               







13:31:09 UTC status = success                                                   













13:31:11 UTC View simulation result at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-95b7981b-725
             1-48e6-bc56-8ce8d67cbba7'.                                         


















13:31:12 UTC Loading results from data/simulation_lorentz.hdf5                  











Result Visualization


After the simulation is complete, we inspect the focusing capability of the gradient-index lens by plotting the field distributions on the same plane as before. Again, at freq0, the result is consistent with the previous simulation.






sim_data_lorentz.plot_field("field_xz", "Ex", vmin=-15, vmax=15)
plt.show()

























sim_data_lorentz.plot_field("field_xz", "E", "abs^2", vmin=0, vmax=300)
plt.show()



























Custom Anisotropic Medium


With CustomAnisotropicMedium, one can define a diagonally anisotropic medium in which each component is spatially varying. The usage of CustomAnisotropicMedium is similar to AnisotropicMedium, except that each of its components is a spatially varying medium.


In the following, we illustrate how to define a spatially varying anisotropic medium. We define a medium that is dispersive in the zz-component, but non-dispersive in the xx and yy components. At freq0, the medium behaves isotropically, and is identical to the previous custom non-dispersive medium.




Setup Simulation with Spatially Varying Lorentz Medium






mat_anisotropic = td.CustomAnisotropicMedium(xx=mat_custom1, yy=mat_custom2, zz=mat_lorentz)











lens_anisotropic = lens.copy(update={"medium": mat_anisotropic})
sim_anisotropic = sim.copy(update={"structures": [lens_anisotropic]})











Visualize the Simulation and Gradient Index Distribution


Let’s visualize the simulation on the same plane as before. To view the permittivity distribution at freq0, we set the kwarg freq=freq0. The permittivity distribution at freq0 looks the same as the previous custom non-dispersive medium example.






sim_anisotropic.plot_eps(x=0, freq=freq0)
plt.show()

























sim_anisotropic.plot_eps(z=t / 2, freq=freq0)
plt.show()

























Submit Simulation Job


Submit the simulation job to the server.






sim_data_anisotropic = web.run(
    sim_anisotropic,
    task_name="gradient_index_lens_anisotropic",
    path="data/simulation_anisotropic.hdf5",
    verbose=True,
)







13:31:19 UTC Created task 'gradient_index_lens_anisotropic' with resource_id    
             'fdve-e73d69fa-31d1-488f-a231-7281d3c2e9a4' and task_type 'FDTD'.  







             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-e73d69fa-31d
             1-488f-a231-7281d3c2e9a4'.                                         







             Task folder: 'default'.                                            


















13:31:49 UTC Estimated FlexCredit cost: 0.877. This assumes the FDTD solver runs
             for the full simulation time; if early shutoff is reached, the     
             billed cost can be lower. Use 'web.real_cost(task_id)' to get the  
             billed FlexCredit cost after a simulation run.                     







13:31:50 UTC status = queued                                                    







             To cancel the simulation, use 'web.abort(task_id)' or              
             '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.                                                             












13:32:18 UTC status = preprocess                                                













13:32:41 UTC starting up solver                                                 







             running solver                                                     












13:33:40 UTC early shutoff detected at 0%, exiting.                             













             status = queued                                                    












13:33:55 UTC status = preprocess                                                







13:33:57 UTC status = running                                                   







13:36:17 UTC status = postprocess                                               







13:36:22 UTC status = success                                                   













13:36:24 UTC View simulation result at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-e73d69fa-31d
             1-488f-a231-7281d3c2e9a4'.                                         


















13:36:25 UTC Loading results from data/simulation_anisotropic.hdf5              











Result Visualization


After the simulation is complete, we inspect the focusing capability of the gradient-index lens by plotting the field distributions on the same plane as before. Again, at freq0, the result is consistent with the previous simulation.






sim_data_anisotropic.plot_field("field_xz", "Ex", vmin=-15, vmax=15)
plt.show()

























sim_data_anisotropic.plot_field("field_xz", "E", "abs^2", vmin=0, vmax=300)
plt.show()



























Notes:


    

  • By default, subpixel averaging is off on the surface and inside the structure made of custom medium. To apply subpixel on the interface of the structure, including exterior boundary and intersection interfaces with other structures, please set the field subpixel=True in custom medium. Here is an example on setting this option in the custom dispersive medium:
    • 







mat_lorentz = td.CustomLorentz(
    eps_inf=eps_inf_dataset,
    coeffs=((delep_dataset, f0_dataset, gamma_dataset),),
    subpixel=True,
)













          

          

            Example library
          

          

          
          

          

            
            
              
              Previous Tutorial
            
            
            
            
              Next Tutorial