Luneburg lens is a prototypical gradient index (GRIN) optical component. A classical Luneburg lens is a spherical lens with a spatially varying refractive index profile following $n(r)=n_0\sqrt{2-(r/R)^2}$, where $r$ is the radial distance, $R$ is the radius of the lens, and $n_0$ is the refractive index of the ambient environment. Plane wave incident on a Luneburg lens will be focused to a point on the surface of the lens. Compared to a usual refractive lens, Luneburg lens is abberation-free and coma-free, which enables a wide range of applications in modern optical systems.
However, it is practically difficult to construct such a lens due to the required gradient index distribution. In the microwave regime, high-gain antennas based on a Luneburg lens design can be achieved by, for example, using concentric ceremics shells with different densities. In the optical frequencies, such an approach is generally not applicable.
In this notebook, we demonstrate the numerical simulation of a practical 3D optical Luneburg lens. The structure consists of a large number of subwavelength unit cells. Using an effective medium approach, each unit cell can be approximated by a local effective index, which can be tuned by the filling fraction of the dielectric polymer in the unit cell. By varying the filling fraction of each unit cell such that the local effective index follows $n(r)=n_0\sqrt{2-(r/R)^2}$, a Luneburg lens is constructed. This design is adapted from Zhao, Y. Y. et al. Three-dimensional Luneburg lens at optical frequencies. Laser Photonics Rev. 10, 665–672 (2016). In the simulation, a linearly polarized plane wave is launched towards the Luneburg lens. Through the visualization of the field distribution, the focusing capability of the lens can be assessed. We also compare the practical Luneburg lens design with the idealized case, which is simulated using CustomMedium. The comparison result shows the practical Luneburg lens design is very optimal.
For more lens examples such as the metalens in the visible frequency range and the spherical Fresnel lens, please visit our examples page.
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
import tidy3d as td
import tidy3d.web as web
from tidy3d import ScalarFieldDataArray
from tidy3d import PermittivityDataset
Simulation Setup¶
The Luneburg lens is designed to work in the mid-IR frequencies around 6.25 $\mu m$. Therefore, the spectrum of the source in the simulation is around this wavelength.
lda0 = 6.25 # central wavelength
ldas = np.linspace(5.25, 7.25, 10) # simulation wavelength range
freq0 = td.C_0 / lda0 # central frequency
freqs = td.C_0 / ldas # simulation frequency range
fwidth = 0.5 * (
np.max(freqs) - np.min(freqs)
) # width of the frequency gaussian distribution
The period of the unit cell is 2 $\mu m$.
a = 2 # period of the unit cell
Only two materials are involved in this model -- the dielectric polymer and air. The polymer has a refractive index of 1.52 in mid-IR.
n_d = 1.52 # refractive index of the dielectric polymer
n_0 = 1 # refractive index of air
dielectric = td.Medium(permittivity=n_d**2)
air = td.Medium(permittivity=n_0**2)
The unit cell is a simply cubic with voids. We define the width of the polymer frames to be $w$. By tuning $w$ from 0 to 0.5$a$, the fillig fraction $f$ is changed from 0 to 1. Since the Luneburg lens structure consists of a large number of unit cells with varying geometries, it is convenient to define a function called build_unit_cell that takes in $w$ and the center coordinates and returns a unit cell structure. This function can then be called systematically later to construct the whole lens.
def build_unit_cell(w, x, y, z):
unit_cell = []
unit_cell.append(
td.Structure(
geometry=td.Box(center=(x, y, z), size=(a, a, a)), medium=dielectric
)
)
unit_cell.append(
td.Structure(
geometry=td.Box(center=(x, y, z), size=(a - 2 * w, a - 2 * w, a)),
medium=air,
)
)
unit_cell.append(
td.Structure(
geometry=td.Box(center=(x, y, z), size=(a - 2 * w, a, a - 2 * w)),
medium=air,
)
)
unit_cell.append(
td.Structure(
geometry=td.Box(center=(x, y, z), size=(a, a - 2 * w, a - 2 * w)),
medium=air,
)
)
return unit_cell
In this particular design, the radius of the Luneburg lens is 20 $\mu m$, i.e. 10 unit cells. The effective index at each site is the discretized version of $n(r)=n_0\sqrt{2-(r/R)^2}$. For our unit cell, the effective index scales approximately linearly with the filling fraction. Therefore, it is easy to obtain the desirable filling fraction at each unit cell.
N = 10 # number of unit cells from 0 to R
R = N * a # radius of the Luneburg lens
r = np.linspace(a / 2, R - a / 2, N) # distance of each unit cell to the origin
n_r = np.sqrt(2 - (r / R) ** 2) # desired effective index at each unit cell
f_r = (n_r - n_0) / (n_d - n_0) # corresponding filling fraction at each unit cell
Since we wrote the build_unit_cell function with $w$ as the input argument, we need to know $w$ at each unit cell. This can be done through the relationship that $f=1-\frac{(a-2w)^2(a+4w)}{a^3}$. Here we simply use the fsolve function from the Scipy library to solve for $w$ with the given $f$.
w_r = np.zeros(N) # width of the polymer frame at each unit cell
# solve for w_r from f_r
for i, f in enumerate(f_r):
def func(w):
return 1 - (a - 2 * w) ** 2 * (a + 4 * w) / a**3 - f
solution = fsolve(func, 0.5)
w_r[i] = solution.item()
With the obtained $w$ as a function of radial distance, we are ready to construct the Luneburg lens strucutre. This can be done easily through calling the build_unit_cell function in a neasted loops over the $x$, $y$, and $z$ coordinates of each unit cell.
Thanks to the symmetries, we only need to build a quater of the Luneburg lens structure. This drastically reduces the total number of Structures as well as the number of grid points of the model.
luneburg_lens = []
for x in r:
for y in r:
for z in np.linspace(-R + a / 2, R - a / 2, 2 * N):
r_local = np.sqrt(
x**2 + y**2 + z**2
) # radial distance of the unit cell
# build an unit cell if the radial distance is smaller or equal to the lens radius
if r_local <= R:
luneburg_lens.extend(
build_unit_cell(np.interp(r_local, r, w_r), x, y, z)
)
Next, we define a PlaneWave source and two FieldMonitors, one in the $xz$ plane at $y=0$ and one in the $xy$ plane at $z=R$. The focus is supposed to be at $z=R$ so we can visualize the focal spot through the second FieldMonitor.
# define a plane wave source
plane_wave = td.PlaneWave(
source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
size=(td.inf, td.inf, 0),
center=(0, 0, -11 * a),
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 field monitor in the xy plane at z=R
monitor_field_xy = td.FieldMonitor(
center=[0, 0, R], size=[td.inf, td.inf, 0], freqs=[freq0], name="field_xy"
)
Finally, we are ready to define the simulation with the above defined structures, source, and monitors.
# simulation domain size
Lx, Ly, Lz = 25 * a, 25 * a, 30 * a
sim_size = (Lx, Ly, Lz)
run_time = 3e-12 # simulation run time
# define simulation
sim = td.Simulation(
center=(0, 0, 2 * a),
size=sim_size,
grid_spec=td.GridSpec.auto(min_steps_per_wvl=30, wavelength=lda0),
structures=luneburg_lens,
sources=[plane_wave],
monitors=[monitor_field_xz, monitor_field_xy],
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.
)
To visualize the simulation, especially the structures, we use the plot function of Simulation. More specifically, we can look at the slice through the $xz$ plane at $y=0$. From the plot, we can see the unit cells with varying filling fraction as well as the source and monitors are set up correctly.
sim.plot(y=0)
plt.show()
Submit the simulation job to the server. We name the simulation data sim_data_practical to distinguish the simulation data for the idealizd Luneburg lens later on.
sim_data_practical = web.run(
sim,
task_name="practical_luneburg_lens",
path="data/simulation_data.hdf5",
verbose=True,
)
09:56:26 Eastern Daylight Time Created task 'practical_luneburg_lens' with task_id 'fdve-28c6f122-a548-4fd4-ade7-45393936e1b9' and task_type 'FDTD'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId =fdve-28c6f122-a548-4fd4-ade7-45393936e1b9'.
Output()
09:56:36 Eastern Daylight Time 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.
Output()
09:56:45 Eastern Daylight Time status = preprocess
09:56:53 Eastern Daylight Time Maximum FlexCredit cost: 0.205. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
starting up solver
running solver
Output()
09:57:22 Eastern Daylight Time early shutoff detected at 48%, exiting.
status = postprocess
Output()
09:57:25 Eastern Daylight Time status = success
09:57:26 Eastern Daylight Time View simulation result at 'https://tidy3d.simulation.cloud/workbench?taskId =fdve-28c6f122-a548-4fd4-ade7-45393936e1b9'.
Output()
09:57:28 Eastern Daylight Time loading simulation from data/simulation_data.hdf5
Result Visualization¶
First, let's visualize the field intensity as well as $E_x$ in the $xz$ plane. A strong intensity spot is observed around $z=R$, indicating the good focus capability of the designed Luneburg lens. From the $E_x$ plot, we can see the wave front gradually converges due to the locally varying effective refractive index.
fig, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(10, 5))
# plot field intensity at the xz plane
sim_data_practical.plot_field(
field_monitor_name="field_xz", field_name="E", val="abs^2", ax=ax1, vmin=0, vmax=25
)
# plot Ex at the xz plane
sim_data_practical.plot_field(
field_monitor_name="field_xz", field_name="Ex", ax=ax2, vmin=-5, vmax=5
)
plt.show()
To further investigate the focusing, we plot several intensity profiles around the focus. The intensity and width of the focus can be clearly observed from this plot.
zs = [21, 23, 25, 27] # z coordinates of the field profile slices
# plot intensity profiles
fig, ax = plt.subplots()
for z in zs:
I = sim_data_practical.get_intensity("field_xz").sel(z=z, method="nearest")
I.plot(ax=ax, label=f"z={z} $\mu m$")
ax.legend()
ax.set_title("Field intensity")
plt.show()
Lastly, plot the intensity distribution in the $xy$ plane to visualize the focal spot shape.
sim_data_practical.plot_field(
field_monitor_name="field_xy", field_name="E", val="abs^2", vmin=0, vmax=25
)
plt.show()
Comparison to the Ideal Luneburg Lens¶
As a comparison, we simulate an ideal Luneburg lens whose local refractive index follows exactly according to $n(r)=n_0\sqrt{2-(r/R)^2}$. This can be realized in Tidy3D using the CustomMedium, which conveniently defines a spatially varying refractive index profile within one structure.
First, we define the spatial and frequency grids and define a 4-dimensional array that stores the refractive index.
Nx, Ny, Nz, Nf = 100, 100, 100, 1 # number of grid points along each dimension
X = np.linspace(-R, R, Nx) # x grid
Y = np.linspace(-R, R, Ny) # y grid
Z = np.linspace(-R, R, Nz) # z grid
freqs = [freq0] # frequency grid
# define coordinate array
x_mesh, y_mesh, z_mesh, freq_mesh = np.meshgrid(X, Y, Z, freqs, indexing="ij")
r_mesh = np.sqrt(x_mesh**2 + y_mesh**2 + z_mesh**2)
# index of refraction array
# assign the refractive index value to the array according to the desired profile
n_data = np.ones((Nx, Ny, Nz, Nf))
n_data[r_mesh <= R] = np.sqrt(2 - (r_mesh[r_mesh <= R] / R) ** 2)
The numpy array is converted to a ScalarFieldDataArray that labels the coordinate. A custome medium is then defined using the classmethod td.CustomMedium.from_nk. Finally, the lens structure is defined.
# convert to dataset array
n_dataset = ScalarFieldDataArray(n_data, coords=dict(x=X, y=Y, z=Z, f=freqs))
# define custom medium based on the dataset
mat_custom = td.CustomMedium.from_nk(n_dataset, interp_method="nearest")
# define the ideal luneburg lens structure
lens = td.Structure(geometry=td.Sphere(radius=R), medium=mat_custom)
Most simulation setup from the previous section can be readily reused. Therefore, we simply copy the previous simulation and only update the structures to contain the ideal Luneburg lens. Using the plot_eps method we can visualize the gradient permittivity distribution.
sim = sim.copy(update={"structures": [lens]})
sim.plot_eps(y=0)
plt.show()
Submit the simulation job to the server.
sim_data_ideal = web.run(
sim, task_name="ideal_luneburg_lens", path="data/simulation_data.hdf5", verbose=True
)
09:57:40 Eastern Daylight Time Created task 'ideal_luneburg_lens' with task_id 'fdve-999b5e17-b15f-455b-873c-2f12e20240c1' and task_type 'FDTD'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId =fdve-999b5e17-b15f-455b-873c-2f12e20240c1'.
Output()
09:57:44 Eastern Daylight Time 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.
Output()
09:57:50 Eastern Daylight Time status = preprocess
09:57:53 Eastern Daylight Time Maximum FlexCredit cost: 0.040. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
starting up solver
running solver
Output()
09:58:06 Eastern Daylight Time early shutoff detected at 20%, exiting.
status = postprocess
Output()
09:58:11 Eastern Daylight Time status = success
View simulation result at 'https://tidy3d.simulation.cloud/workbench?taskId =fdve-999b5e17-b15f-455b-873c-2f12e20240c1'.
Output()
09:58:13 Eastern Daylight Time loading simulation from data/simulation_data.hdf5
We perform the same postprocessing visualizations as in the previous section. The results of the designed Luneburg lens are very similar to the idealized case, confirming the validity of the design using the effective index approach.
fig, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(10, 5))
# plot field intensity at the xz plane
sim_data_ideal.plot_field(
field_monitor_name="field_xz", field_name="E", val="abs^2", ax=ax1, vmin=0, vmax=25
)
# plot Ex at the xz plane
sim_data_ideal.plot_field(
field_monitor_name="field_xz", field_name="Ex", ax=ax2, vmin=-5, vmax=5
)
plt.show()
# plot intensity profiles
fig, ax = plt.subplots()
for z in zs:
I = sim_data_ideal.get_intensity("field_xz").sel(z=z, method="nearest")
I.plot(ax=ax, label=f"z={z} $\mu m$")
ax.legend()
ax.set_title("Field intensity")
plt.show()
sim_data_ideal.plot_field(
field_monitor_name="field_xy", field_name="E", val="abs^2", vmin=0, vmax=25
)
plt.show()
Lastly, as a direct comparison, we plot the field intensity around the focus for the practical and idealized Luneburg lens together. The resutls are nearly identical, which validates the design of the practical Lunebug lens.
# compare ideal and practical Luneburg lens intensity profiles at fix z
z = 21
fig, ax = plt.subplots()
I_practical = sim_data_practical.get_intensity("field_xz").sel(z=z, method="nearest")
I_practical.plot(ax=ax, label="Practical Luneburg lens")
I_ideal = sim_data_ideal.get_intensity("field_xz").sel(z=z, method="nearest")
I_ideal.plot(linestyle="--", ax=ax, label="Ideal Luneburg lens")
ax.legend()
ax.set_title("Field intensity")
plt.show()
