Defining gyrotropic materials#

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

ε = ε_{u} + i \frac{σ_{u}}{ω} + i ε_{0} G,

where $ε_{u}$ and $σ_{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,

\begin{array}{r} G = (\begin{array}{c} 0 & g_{z} & - g_{y} \\ - g_{z} & 0 & g_{x} \\ g_{y} & - g_{x} & 0 \end{array}) \end{array}

for a gyration vector $g = (\begin{matrix} g_{x} & g_{y} & g_{z} \end{matrix})$ . Denoting the electromagnetic frequency of interest as $ω_{0}$ , we can approximate the complex permittivity of material in the vicinity of $ω_{0}$ as

ε \approx ε_{u} + i \frac{1}{ω} (σ_{u} + ε_{0} G ω_{0}) .

Thus, the gyrotropic effect of a medium can be modeled by providing a modified conductivity tensor $σ = σ_{u} + ε_{0} G ω_{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 $ω_{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 $g = 0$ and another one with $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.

[1]:

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

# 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.

[2]:

# 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.

[3]:

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

Gyration vector along the propagation axis.

[4]:

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.

[5]:

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.

[6]:

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

[7]:

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.

[8]:

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

[13:07:16] 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.

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.

[9]:

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",
)

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

[10]:

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

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

[11]:

# 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,
)

Visualize simulation for setup confirmation.

[12]:

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

../_images/notebooks_Gyrotropic_25_0.png

Results#

Submit and run simulations on the server.

[13]:

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)

[13:07:17] Created task 'gyrotropic_reference' with task_id        webapi.py:188
           'fdve-0a3b96df-5a09-4192-82fb-ec14c13e0a98v1'.

           View task using web UI at                               webapi.py:190
           'https://tidy3d.simulation.cloud/workbench?taskId=fdve-              
           0a3b96df-5a09-4192-82fb-ec14c13e0a98v1'.

[13:07:18] status = queued                                         webapi.py:361

[13:07:26] status = preprocess                                     webapi.py:355

[13:07:30] Maximum FlexCredit cost: 0.025. 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.

[13:07:36] early shutoff detected, exiting.                        webapi.py:404

           status = postprocess                                    webapi.py:419

[13:07:40] status = success                                        webapi.py:426

[13:07:41] loading SimulationData from                             webapi.py:590
           data/simulation_data_gyrotropic_ref.hdf5

           Created task 'gyrotropic' with task_id                  webapi.py:188
           'fdve-cc50998f-f091-451f-a0c9-a4876a9c7b0av1'.

           View task using web UI at                               webapi.py:190
           'https://tidy3d.simulation.cloud/workbench?taskId=fdve-              
           cc50998f-f091-451f-a0c9-a4876a9c7b0av1'.

[13:07:42] status = queued                                         webapi.py:361

[13:07:51] status = preprocess                                     webapi.py:355

[13:07:55] Maximum FlexCredit cost: 0.025. 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.

[13:08:01] early shutoff detected, exiting.                        webapi.py:404

           status = postprocess                                    webapi.py:419

[13:08:05] status = success                                        webapi.py:426

           loading SimulationData from                             webapi.py:590
           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

ρ = \frac{π}{λ_{0}} (n_{-} - n_{+}) .

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

[14]:

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.

[15]:

# 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.

[16]:

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.

[17]:

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()

../_images/notebooks_Gyrotropic_37_0.png

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.

[18]:

# 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()

../_images/notebooks_Gyrotropic_39_0.png

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

[ ]:

Defining gyrotropic materials

Contents

Defining gyrotropic materials#

Simulation setup#

Results#