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

\[\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,

\[\begin{split}G = \begin{pmatrix} 0 & g_z & -g_y \\ -g_z & 0 & g_x \\ g_y & -g_x & 0\end{pmatrix}\end{split}\]

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\).

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.

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

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)

[16:57:08] Created task 'gyrotropic_reference' with task_id                                           webapi.py:186
           'fdve-010efdb5-62aa-4da9-9fea-d51749934a32v1'.

           View task using web UI at                                                                  webapi.py:188
           'https://tidy3d.simulation.cloud/workbench?taskId=fdve-010efdb5-62aa-4da9-9fea-d51749934a3              
           2v1'.

[16:57:12] status = queued                                                                            webapi.py:321

[16:57:19] status = preprocess                                                                        webapi.py:315

[16:57:26] Maximum FlexCredit cost: 0.025. Use 'web.real_cost(task_id)' to get the billed FlexCredit  webapi.py:338
           cost after a simulation run.

           starting up solver                                                                         webapi.py:342

           running solver                                                                             webapi.py:352

[16:57:33] early shutoff detected, exiting.                                                           webapi.py:366

           status = postprocess                                                                       webapi.py:383

[16:57:39] status = success                                                                           webapi.py:390

[16:57:41] loading SimulationData from data/simulation_data_gyrotropic_ref.hdf5                       webapi.py:568

[16:57:41] Created task 'gyrotropic' with task_id 'fdve-99c887ed-1681-492e-a374-2781ca415175v1'.      webapi.py:186

           View task using web UI at                                                                  webapi.py:188
           'https://tidy3d.simulation.cloud/workbench?taskId=fdve-99c887ed-1681-492e-a374-2781ca41517              
           5v1'.

[16:57:43] status = queued                                                                            webapi.py:321

[16:57:49] status = preprocess                                                                        webapi.py:315

[16:57:56] Maximum FlexCredit cost: 0.025. Use 'web.real_cost(task_id)' to get the billed FlexCredit  webapi.py:338
           cost after a simulation run.

           starting up solver                                                                         webapi.py:342

           running solver                                                                             webapi.py:352

[16:58:02] status = postprocess                                                                       webapi.py:383

[16:58:08] status = success                                                                           webapi.py:390

[16:58:09] loading SimulationData from data/simulation_data_gyrotropic.hdf5                           webapi.py:568

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.

[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

[ ]:

Defining gyrotropic materials

Contents

Defining gyrotropic materials#

Simulation setup#

Results#