How to define fully anisotropic materials in Tidy3D FDTD¶
Tidy3D's capabilities include modeling of non-dispersive fully anisotropic materials. In this tutorial we explain how to set up and run a simulation containing such materials. Specifically, we will consider scattering of a plane wave from fully anisotropic dielectric spheres and compare results to an equivalent setup containing only diagonally anisotropic materials.
Simulation setup¶
# standard python imports
import numpy as np
import xarray as xr
import matplotlib.pylab as plt
# tidy3D import
import tidy3d as td
from tidy3d.constants import C_0
import tidy3d.web as web
We start with defining the simulation domain as a cube with side length of 1.6 um, and selecting 11 sample frequencies in the vicinity of the wavelength of interest $\lambda = 1$ um.
sim_size = (1.6, 1.6, 1.6)
wvl_um = 1
freq0 = td.C_0 / wvl_um
fwidth = freq0 / 10
num_freqs = 11
freqs = np.linspace(freq0 - fwidth / 2, freq0 + fwidth / 2, num_freqs)
find0 = num_freqs // 2
run_time = 10 / fwidth
Next we define main permittivity and conductivity components of medium we would like to simulate.
permittivity_main = (15, 10, 5)
conductivity_main = (0.01, 0.02, 0.03)
We will create 3 different anisotropic media to use and compare in 3 different simulations. First, a diagonally anisotropic medium that will be used in a reference simulation is created using the familiar AnisotropicMedium class.
# diagonally anisotropic material
medium_diag = td.AnisotropicMedium(
xx=td.Medium(permittivity=permittivity_main[0], conductivity=conductivity_main[0]),
yy=td.Medium(permittivity=permittivity_main[1], conductivity=conductivity_main[1]),
zz=td.Medium(permittivity=permittivity_main[2], conductivity=conductivity_main[2]),
)
Next, we would like to create two fully anisotropic media obtained as rotations of the above diagonally anisotropic medium around x and z axes, respectively. Fully anisotropic materials are created using class FullyAnisotropicMedium, which requires full permittivity and conductivity tensor during initialization. Tidy3D provides a convenience class RotationAroundAxis that can be used for transformation of material property tensors. Using it we can define a fully anisotropic medium in two ways. First, one can use this class directly to perform rotation of material property tensors, and then use the resulting tensors to initialize a FullyAnisotropicMedium as following.
# Fully anisotropic media obtained as a rotation around x axis by angle phi
phi = np.pi / 3
# create diagonal tensor
permittivity_phi = np.diag(permittivity_main)
conductivity_phi = np.diag(conductivity_main)
# create rotation around axis
rotation_around_x = td.RotationAroundAxis(axis=(1, 0, 0), angle=phi)
# transform material property tensors
permittivity_phi = rotation_around_x.rotate_tensor(permittivity_phi)
conductivity_phi = rotation_around_x.rotate_tensor(conductivity_phi)
# define a fully anisotropic medium
medium_phi = td.FullyAnisotropicMedium(
permittivity=permittivity_phi, conductivity=conductivity_phi
)
The alternative way is to use a shortcut class method FullyAnisotropicMedium.from_diagonal() that performs the necessary tensor transformations internally. Note that this class method follows the signature of AnisotropicMedium class, that is, it expects main components of material properties defined as Medium's.
# Fully anisotropic media obtained as a rotation around z axis by angle theta
theta = np.pi / 7
# create rotation around axis
rotation_around_z = td.RotationAroundAxis(axis=(0, 0, 1), angle=theta)
# define a fully anisotropic medium as a transformation of diagonally anisotropic medium
medium_theta = td.FullyAnisotropicMedium.from_diagonal(
xx=td.Medium(permittivity=permittivity_main[0], conductivity=conductivity_main[0]),
yy=td.Medium(permittivity=permittivity_main[1], conductivity=conductivity_main[1]),
zz=td.Medium(permittivity=permittivity_main[2], conductivity=conductivity_main[2]),
rotation=rotation_around_z,
)
One can confirm that the three created media have identical main components by visualizing their complex refractive indexes using the built-in function .plot(). Note that in case of fully anisotropic materials the direction of each main component is provided in figure legend as they are no longer aligned with Cartesian directions.
# Visualize material properties
fig, ax = plt.subplots(1, 3, figsize=(14, 4))
medium_diag.plot(freqs, ax=ax[0])
medium_phi.plot(freqs, ax=ax[1])
medium_theta.plot(freqs, ax=ax[2])
plt.show()
Create three spheres of the same size but made of different anisotropic materials.
sphere_diag = td.Structure(geometry=td.Sphere(center=(0, 0, 0), radius=0.3), medium=medium_diag)
sphere_phi = sphere_diag.updated_copy(medium=medium_phi)
sphere_theta = sphere_diag.updated_copy(medium=medium_theta)
To set up three equivalent simulations we use TFSF sources injecting plane waves whose propagation fronts are rotated in the same way as the material properties of the three anisotropic media.
source_diag = td.TFSF(
center=(0, 0, 0),
size=(1.1, 1.1, 1.1),
source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
injection_axis=0,
direction="+",
name="tfsf",
angle_phi=0,
angle_theta=0,
pol_angle=0,
)
source_phi = source_diag.updated_copy(angle_phi=phi)
source_theta = source_diag.updated_copy(angle_theta=theta)
For a quantitative comparison of results of the three equivalent simulation we monitor total scattered flux from the sphere and field distribution in two orthogonal planes.
mnt_flux = td.FluxMonitor(center=(0, 0, 0), size=(1.3, 1.3, 1.3), freqs=freqs, name="flux")
mnt_field_yz = td.FieldMonitor(center=(0, 0, 0), size=(0, td.inf, td.inf), freqs=freqs, name="field_yz")
mnt_field_xy = td.FieldMonitor(center=(0, 0, 0), size=(td.inf, td.inf, 0), freqs=freqs, name="field_xy")
PML boundary conditions are placed on each side of the simulation domain and a simple uniform grid is used for grid discretization.
boundary_spec = td.BoundarySpec.all_sides(boundary=td.PML())
grid_spec = td.GridSpec.uniform(dl=0.005)
Combining all the components defined above we create three equivalent simulations which are rotations of each other.
sim_diag = td.Simulation(
center=(0, 0, 0),
size=sim_size,
grid_spec=grid_spec,
structures=[sphere_diag],
sources=[source_diag],
monitors=[mnt_flux, mnt_field_yz, mnt_field_xy],
run_time=run_time,
boundary_spec=boundary_spec,
)
sim_phi = sim_diag.updated_copy(structures=[sphere_phi], sources=[source_phi])
sim_theta = sim_diag.updated_copy(structures=[sphere_theta], sources=[source_theta])
Visualize simulation setups for confirmation.
fig, ax = plt.subplots(1, 3, figsize=(10, 3))
sim_diag.plot(z=0, ax=ax[0])
sim_phi.plot(z=0, ax=ax[1])
sim_theta.plot(z=0, ax=ax[2])
plt.tight_layout()
plt.show()
Results¶
Submit and run simulations on the server.
sim_data_diag = web.run(simulation=sim_diag, task_name="fully_anisotropic_diag", path="data/simulation_data_diag.hdf5", verbose=True)
sim_data_phi = web.run(simulation=sim_phi, task_name="fully_anisotropic_phi", path="data/simulation_data_phi.hdf5", verbose=True)
sim_data_theta = web.run(simulation=sim_theta, task_name="fully_anisotropic_theta", path="data/simulation_data_theta.hdf5", verbose=True)
View task using web UI at webapi.py:188 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-78d231c2-61bf-4a03-b57b-7f61b6e14d1 8v1'.
Output()
Output()
Output()
Output()
Output()
View task using web UI at webapi.py:188 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-23947891-340a-4ea7-8540-86162eeb55d 9v1'.
Output()
Output()
Output()
Output()
Output()
View task using web UI at webapi.py:188 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-169281b2-d8df-4764-a7eb-1e815dd6b28 cv1'.
Output()
Output()
Output()
Output()
Output()
To confirm that the three simulation produces close results we first consider the total scattered flux from the sphere. The slight difference in obtained results is attributed to the fact that currently Tidy3D does not support subpixel averaging for fully anisotropic materials.
# retrive flux values from flux monitor
flux_diag = sim_data_diag["flux"].flux
flux_phi = sim_data_phi["flux"].flux
# scale power by injection angle
flux_theta = sim_data_theta["flux"].flux * np.cos(theta)
# visualize
flux_diag.plot()
flux_phi.plot()
flux_theta.plot()
plt.legend(["diag", "phi", "theta"])
plt.show()
Now we look at the angular field intensity distributions around the sphere at several distances from the sphere's center.
# field sample locations
sample_radius = [0.1, 0.2, 0.4, 0.5]
sample_angle = np.linspace(0, 2 * np.pi, 100)
When comparing the reference (diagonally anisotropic) simulation with the simulation obtained by rotation around $x$ axis we use the fact that field distributions must coincide in the $yz$ plane when appropriately rotated.
# get field intensity distribution at central frequency
int_diag = sim_data_diag.get_intensity("field_yz").isel(x=0, f=find0)
int_full = sim_data_phi.get_intensity("field_yz").isel(x=0, f=find0)
fig, ax = plt.subplots(len(sample_radius), 1, constrained_layout=True, figsize=(8, 8))
for i, r in enumerate(sample_radius):
# sample points in reference simulation
y_diag = xr.DataArray(r * np.cos(sample_angle), dims="u")
z_diag = xr.DataArray(r * np.sin(sample_angle), dims="u")
# rotated sample points in fully anisotropic simulation
y_full = xr.DataArray(r * np.cos(sample_angle + phi), dims="u")
z_full = xr.DataArray(r * np.sin(sample_angle + phi), dims="u")
# interpolate at sample points
int_sampled_diag = int_diag.interp(y=y_diag, z=z_diag)
int_sampled_full = int_full.interp(y=y_full, z=z_full)
# plot comparisons
ax[i].plot(sample_angle, int_sampled_diag.data)
ax[i].plot(sample_angle, int_sampled_full.data)
ax[i].set_xlabel("Angle [rad]")
ax[i].legend(["Reference", "Fully Anisotropic"])
ax[i].set_title(f"Electric field intensity at r = {r}")
plt.show()
When comparing the reference (diagonally anisotropic) simulation with the simulation obtained by rotation around $z$ axis we use the fact that field distributions must coincide in the $xy$ plane when appropriately rotated.
# get field intensity distribution at central frequency
int_diag = sim_data_diag.get_intensity("field_xy").isel(z=0, f=find0)
# scale power by the injection angle
int_full = sim_data_theta.get_intensity("field_xy").isel(z=0, f=find0) * np.cos(theta)
fig, ax = plt.subplots(len(sample_radius), 1, constrained_layout=True, figsize=(8, 8))
for i, r in enumerate(sample_radius):
# sample points in reference simulation
x_diag = xr.DataArray(r * np.cos(sample_angle), dims="u")
y_diag = xr.DataArray(r * np.sin(sample_angle), dims="u")
# rotated sample points in fully anisotropic simulation
x_full = xr.DataArray(r * np.cos(sample_angle + theta), dims="u")
y_full = xr.DataArray(r * np.sin(sample_angle + theta), dims="u")
# interpolate at sample points
int_sampled_diag = int_diag.interp(x=x_diag, y=y_diag)
int_sampled_full = int_full.interp(x=x_full, y=y_full)
# plot comparisons
ax[i].plot(sample_angle, int_sampled_diag.data)
ax[i].plot(sample_angle, int_sampled_full.data)
ax[i].set_xlabel("Angle [rad]")
ax[i].legend(["Reference", "Fully Anisotropic"])
ax[i].set_title(f"Electric field intensity at r = {r}")
plt.show()
Again, we see that equivalent simulations produce very close results, and the small differences are attributed to the lack of subpixel averaging for fully anisotropic materials.