TIDY3D
LEARNING CENTER

Lectures
Tutorials
Applications
Complete this course to earn 50 FlexCredits Credits are added to your account when all course items are checked in green.

Tidy3D Python Interface / Mediums / Defining gyrotropic materials

Defining gyrotropic materials

Sign up now to simulate .ipynb

Here we demonstrate the usage of Tidy3D simulation software for modeling of gyrotropic materials, that is, materials with Hermitian permittivity tensors. Currently, this is achieved by using fully anisotropic dieletric medium option with anisymmetric conductivity tensor. Specifically, let us consider modeling a material with total complex permittivity tensor \[\varepsilon = \varepsilon_u + i \frac{\sigma_u}{\omega} + i \varepsilon_0 G,\] where \(\varepsilon_u\) and \(\sigma_u\) are real symmetric matrices describing unperturbed permittivity and conductivity of the material, and \(G\) is a real antisymmetric matrix representing the magneto-optic effect. That is, \[G = \begin{pmatrix} 0 & g_z & -g_y \\ -g_z & 0 & g_x \\ g_y & -g_x & 0\end{pmatrix}\] for a gyration vector \(\boldsymbol{g} = \begin{pmatrix} g_x & g_y & g_z \end{pmatrix}\). Denoting the electromagnetic frequency of interest as \(\omega_0\), we can approximate the complex permittivity of material in the vicinity of \(\omega_0\) as \[\varepsilon \approx \varepsilon_u + i \frac{1}{\omega} \left( \sigma_u + \varepsilon_0 G \omega_0 \right).\] Thus, the gyrotropic effect of a medium can be modeled by providing a modified conductivity tensor \(\sigma = \sigma_u + \varepsilon_0 G \omega_0\) that contains both symmetric and antisymmetric parts. Note that because of this approximation, the simulation results will be most accurate in the vicinity of the frequency of interest \(\omega_0\).

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.

Simulation setup

Let us now demonstrate simulation of plane wave propagation inside a gyrotropic material. Specifically, we will consider two simulations, one with \(\boldsymbol{g} = 0\) and another one with \(\boldsymbol{g} \neq 0\). In the first, reference, simulation the injected plane wave is expected to preserve its polarization, while in the second case the polarization of the injected pulse will undergo rotation due to the gyrotropic properties of the background medium. The rotation degree can be predicted by analytical expressions. In section Results we apply it to the results of the reference simulation and compare with the results of the gyrotropic simulation.

# standard python imports
import matplotlib.pylab as plt
import numpy as np

# tidy3D import
import tidy3d as td
import tidy3d.web as web

We define an elongated simulation domain assuming wave propagation along \(z\) direction and select sample frequencies in the vicinity of the wavelength of interest.

# simulation domain
sim_length = 10
sim_size = [0.1, 0.1, sim_length]

# central frequency
wvl_um = 1
freq0 = td.C_0 / wvl_um
fwidth = freq0 / 10
freqs = np.linspace(freq0 - fwidth, freq0 + fwidth, 5)

run_time = 10 / fwidth

Properties of the unperturbed material.

eps0 = 2
permittivity_unperturbed = np.diag((eps0, eps0, eps0))
conductivity_unperturbed = np.zeros((3, 3))

Gyration vector along the propagation axis.

g = (0, 0, 0.1)
G = np.array([[0, g[2], g[1]], [-g[2], 0, g[0]], [-g[1], -g[0], 0]])

Following the procedure outlined in the beginning of this tutorial we modify the material conductivity such that it accounts for material gyrotropic properties.

conductivity = conductivity_unperturbed + G * (2 * np.pi * freq0) * td.EPSILON_0

Now we define two materials: a simple isotropic medium to use in a reference simulation and a fully anisotropic medium for modeling the magneto-optic effect.

medium_unperturbed = td.Medium(permittivity=eps0)
medium_gyrotropic = td.FullyAnisotropicMedium(
    permittivity=permittivity_unperturbed, conductivity=conductivity
)

Given the simulation geometry we place PML boundary conditions along the propagation axis and periodic boundary conditions in the orthogonal directions

boundary_spec = td.BoundarySpec(
    x=td.Boundary.periodic(),
    y=td.Boundary.periodic(),
    z=td.Boundary.pml(),
)

We will track field distribution along the propagation axis.

mnt_xz = td.FieldMonitor(
    center=(0, 0, 0),
    size=(0, 0, td.inf),
    freqs=freqs,
    name="freq_mnt_xz",
)

Since Tidy3D does not currently support injection of unidirectional plane waves into anisotropic material, we use a uniform current source to initiate a planar (bidirectional) propagating pulse.

source = td.UniformCurrentSource(
    center=(0, 0, -sim_length / 2.5),
    size=(td.inf, td.inf, 0),
    source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
    polarization="Ex",
)
04:53:54 UTC WARNING: The default value of 'current_amplitude_definition' for   
             'UniformCurrentSource' will change from 'density' to 'total' in a  
             future release. To avoid this warning and ensure consistent        
             behavior, please explicitly set                                    
             current_amplitude_definition='density' or                          
             current_amplitude_definition='total' when creating the source.

Use automatic grid discretization based on the requested number of steps per source wavelength.

grid_spec = td.GridSpec.auto(min_steps_per_wvl=30)

Combining defined above components we create two simulation, with and without gyrotropic effect.

# Propagation in an unperturbed medium
sim_unperturbed = td.Simulation(
    size=sim_size,
    grid_spec=grid_spec,
    medium=medium_unperturbed,
    sources=[source],
    monitors=[mnt_xz],
    run_time=run_time,
    boundary_spec=boundary_spec,
)

# Propagation in a gyrotropic medium
sim_gyrotropic = td.Simulation(
    size=sim_size,
    grid_spec=grid_spec,
    medium=medium_gyrotropic,
    sources=[source],
    monitors=[mnt_xz],
    run_time=run_time,
    boundary_spec=boundary_spec,
)
04:53:55 UTC WARNING: The default value of 'current_amplitude_definition' for   
             'UniformCurrentSource' will change from 'density' to 'total' in a  
             future release. To avoid this warning and ensure consistent        
             behavior, please explicitly set                                    
             current_amplitude_definition='density' or                          
             current_amplitude_definition='total' when creating the source.     

             WARNING: The default value of 'current_amplitude_definition' for   
             'UniformCurrentSource' will change from 'density' to 'total' in a  
             future release. To avoid this warning and ensure consistent        
             behavior, please explicitly set                                    
             current_amplitude_definition='density' or                          
             current_amplitude_definition='total' when creating the source.

Visualize simulation for setup confirmation.

_, ax = plt.subplots(1, 1, figsize=(2, 3))
sim_gyrotropic.plot(y=0, ax=ax)
ax.set_aspect("auto")
plt.show()

Results

Submit and run simulations on the server.

data_unperturbed = web.run(
    simulation=sim_unperturbed,
    task_name="gyrotropic_reference",
    path="data/simulation_data_gyrotropic_ref.hdf5",
    verbose=True,
)
data_gyrotropic = web.run(
    simulation=sim_gyrotropic,
    task_name="gyrotropic",
    path="data/simulation_data_gyrotropic.hdf5",
    verbose=True,
)
             Created task 'gyrotropic_reference' with resource_id               
             'fdve-eb7daf00-9301-41bc-a71d-22d7172e30b1' and task_type 'FDTD'.  

             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-eb7daf00-930
             1-41bc-a71d-22d7172e30b1'.                                         

             Task folder: 'default'.                                            







04:53:57 UTC Estimated FlexCredit cost: 0.025. 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.                     







04:53:58 UTC status = success                                                   


















04:53:59 UTC Loading results from data/simulation_data_gyrotropic_ref.hdf5      







             Created task 'gyrotropic' with resource_id                         
             'fdve-4af292bb-ea10-46eb-bb4b-98e6daf3281b' and task_type 'FDTD'.  







             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-4af292bb-ea1
             0-46eb-bb4b-98e6daf3281b'.                                         







             Task folder: 'default'.                                            


















04:54:01 UTC Estimated FlexCredit cost: 0.025. 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.                     







04:54:02 UTC status = success                                                   


















04:54:03 UTC Loading results from data/simulation_data_gyrotropic.hdf5          







Now we compare the simulation results to analytically predicted behavior. It is known that in the gyrotropic medium propagation along gyration axis allows two circularly polarized normal modes with propagation constants given by \[n_{\pm} = \sqrt{n^2 \pm G}\] and the rotatory power (rotation angle per unit length) of the medium is given by \[\rho = \frac{\pi}{\lambda_0} \left(n_- - n_+\right).\]


First we consider central frequency for which we expect simulation results to be the most accurate.






freq_ind = len(freqs) // 2







To obtain reference field distribution we extract numerical data from the reference simulation and apply rotation of polarization as predicted by the theory.






# reference (constant polarization) wave propagation
Ex_ref = (
    data_unperturbed["freq_mnt_xz"].field_components["Ex"].isel({"f": freq_ind, "x": 0, "y": 0})
)

# predicted rotatory power at chosen wavelength
wvl = td.C_0 / freqs[freq_ind]
n_minus = np.sqrt(eps0 - g[2])
n_plus = np.sqrt(eps0 + g[2])
rho = np.pi / wvl * (n_minus - n_plus)

# rotation factor along propagation direction
coord = Ex_ref.coords["z"]
src_coord = source.center[2]
dist = coord.data - src_coord
factor_x = np.cos(rho * dist).real
factor_y = -np.sin(rho * abs(dist)).real

# rotated polarizations
Ex_theory = Ex_ref * factor_x
Ey_theory = Ex_ref * factor_y







Extract simulation data at the chosen frequency.






Ex_num = data_gyrotropic["freq_mnt_xz"].field_components["Ex"].isel({"f": freq_ind, "x": 0, "y": 0})

Ey_num = data_gyrotropic["freq_mnt_xz"].field_components["Ey"].isel({"f": freq_ind, "x": 0, "y": 0})







Plot the comparison.






fig, ax = plt.subplots(2, 1, figsize=(8, 6))
Ex_num.real.plot(ax=ax[0])
Ex_theory.real.plot(ax=ax[0], ls="--")
ax[0].set_ylabel("Re[Ex]")
ax[0].legend(["Simulation", "Theory"])

Ey_num.real.plot(ax=ax[1])
Ey_theory.real.plot(ax=ax[1], ls="--")
ax[1].set_ylabel("Re[Ey]")
ax[1].legend(["Simulation", "Theory"])

plt.tight_layout()
plt.show()





















As one can see the simulation results for magneto-optic effect matches very closely the predicted behavior. Note that this comparison is made at the central wavelength. If we repeat the same comparison away from central wavelength the accuracy of simulation results will degrade because of the approximations made in defining properties of the gyrotropic medium.






# choose frequency index (out of 5)
freq_ind = 0

# reference (constant polarization) wave propagation
Ex_ref = (
    data_unperturbed["freq_mnt_xz"].field_components["Ex"].isel({"f": freq_ind, "x": 0, "y": 0})
)

# predicted rotatory power at chosen wavelength
wvl = td.C_0 / freqs[freq_ind]
n_minus = np.sqrt(eps0 - g[2])
n_plus = np.sqrt(eps0 + g[2])
rho = np.pi / wvl * (n_minus - n_plus)

# rotation factor along propagation direction
coord = Ex_ref.coords["z"]
src_coord = source.center[2]
dist = coord.data - src_coord
factor_x = np.cos(rho * dist).real
factor_y = -np.sin(rho * abs(dist)).real

# rotated polarizations
Ex_theory = Ex_ref * factor_x
Ey_theory = Ex_ref * factor_y

# extract simulation data at the chosen frequency
Ex_num = data_gyrotropic["freq_mnt_xz"].field_components["Ex"].isel({"f": freq_ind, "x": 0, "y": 0})

Ey_num = data_gyrotropic["freq_mnt_xz"].field_components["Ey"].isel({"f": freq_ind, "x": 0, "y": 0})

# plot comparison
fig, ax = plt.subplots(2, 1, figsize=(8, 6))
Ex_num.real.plot(ax=ax[0])
Ex_theory.real.plot(ax=ax[0], ls="--")
ax[0].set_ylabel("Re[Ex]")
ax[0].legend(["Simulation", "Theory"])

Ey_num.real.plot(ax=ax[1])
Ey_theory.real.plot(ax=ax[1], ls="--")
ax[1].set_ylabel("Re[Ey]")
ax[1].legend(["Simulation", "Theory"])

plt.tight_layout()
plt.show()





















Gyrotropic materials are only a special type of the FullyAnisotropicMedium. For more general cases, see the tutorial on defining fully anisotropic materials.








          

          

            Example library
          

          

          
          

          

            
            
              
              Previous Tutorial
            
            
            
            
              Next Tutorial