Inverse design of a GaP photon extractor for nitrogen-vacancy centers in diamond

Designing a light extractor is crucial for enhancing photon collection efficiency in solid-state defect qubit systems, which are valuable for quantum information and metrology applications. In high refractive index materials like diamond, total internal reflection limits the ability to collect photons emitted from qubits. By optimizing light extraction through nanophotonic structures, such as waveguides or planar dielectric structures, one can significantly increase the emission directed toward detectors, improving photoluminescence detection.

In this notebook, we use inverse design to optimize a GaP extractor structure following a similar approach discussed in Srivatsa Chakravarthi, Pengning Chao, Christian Pederson, Sean Molesky, Andrew Ivanov, Karine Hestroffer, Fariba Hatami, Alejandro W. Rodriguez, and Kai-Mei C. Fu, "Inverse-designed photon extractors for optically addressable defect qubits," Optica 7, 1805-1811 (2020) DOI: 10.1364/OPTICA.408611. The optimization aims to maximize the minimum upward flux among three simulations with the dipole source oriented in the \(x\), \(y\), and \(z\) directions. The final design achieves an order-of-magnitude enhancement in upward flux for an in-plane dipole source and over two orders of magnitude for an out-of-plane dipole source compared to the configuration without an extractor.

If you are unfamiliar with inverse design, we also recommend our intro to inverse design tutorials and our primer on automatic differentiation with tidy3d. For another example of light extractor optimization in Tidy3D, see this example.

import autograd.numpy as np
import autograd as ag
import optax
import matplotlib.pyplot as plt
import gdstk
import tidy3d as td
import tidy3d.web as web

Static Simulation Components

In this case, we are focusing on the photonic coupling to the negatively charged sharp zero-phonon line (ZPL) of nitrogen-vacancy centers at 637 nm.

lda0 = 0.637  # central wavelength
freq0 = td.C_0 / lda0  # central frequency
fwidth = freq0 / 10  # width of the source frequency range

For the extractor structure on top of diamond, a high-index material is preferred. GaP with a refractive index of 3.31 is used in this design.

# define materials
n_GaP = 3.31
GaP = td.Medium(permittivity=n_GaP**2)

n_diamond = 2.417
diamond = td.Medium(permittivity=n_diamond**2)

The GaP layer thickness is 250 nm. The design region is very compact, only 1.5 μm by 1.5 μm. The design region is discretized into 10 nm grids. To ensure the FDTD grid is commensurate with the pixel, we will use the automatic nonuniform grid with min_steps_per_wvl = lda0 / (pixel_size * n_GaP).

The minimal feature size we are aiming for is 50 nm, compatible with electron beam lithography fabrication.

h = 0.25  # thickness of the extractor layer
l = 1.5  # size of the design region
pixel_size = 0.01  # pixel size in the design region
min_steps_per_wvl = lda0 / (pixel_size * n_GaP)  # automatic grid size parameter

min_feature = 0.05  # minimal feature size

inf_eff = 1e3  # effective infinity
buffer = 0.5  # buffer spacing
run_time = 4e-13  # simulation run time

Now we define static components in the optimization. This includes the diamond substrate structure, the dipole sources, and the monitor that records the upward flux. The dipole source is located 100 nm below the top diamond surface.

In this case, we will place a FieldMonitor to record the upward radiation flux instead of directly using a FluxMonitor because differentiation with respect to flux monitor data is not supported at the moment.

# define the diamond substrate
substrate = td.Structure(
    geometry=td.Box.from_bounds(rmin=(-inf_eff, -inf_eff, -inf_eff), rmax=(inf_eff, inf_eff, 0)),
    medium=diamond,
)

# define the dipole sources oriented in the x, y, and z directions
dipole_z = -0.1  # dipole depth
pd_x = td.PointDipole(
    center=(0, 0, dipole_z),
    source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
    polarization="Ex",
)

pd_y = pd_x.updated_copy(polarization="Ey")
pd_z = pd_x.updated_copy(polarization="Ez")

# define the field monitor above the extractor to measure the upward flux
monitor_z = h + 0.4  # monitor height
monitor_size = 1.5  # monitor size

field_monitor = td.FieldMonitor(
    center=(0, 0, monitor_z),
    size=(monitor_size, monitor_size, 0),
    freqs=[freq0],
    name="field",
)

Design Region

The design region is a pixellated array of permittivity values between \(\varepsilon=1\) and \(\varepsilon=\varepsilon_{GaP}\). To define it, we use a density array whose element values are between 0 to 1, apply a conic smoothing filter, a tanh projection function, and finally scale it linearly to \(\varepsilon=1\) and \(\varepsilon=\varepsilon_{GaP}\).

from tidy3d.plugins.autograd import rescale, make_filter_and_project

# define the conic filter and tanh projection
filter_project_fn = make_filter_and_project(radius=min_feature, dl=pixel_size)


def get_density(params: np.ndarray, beta: float) -> np.ndarray:
    """Get the density of the material in the design region as a function of optimization parameters."""
    params = (params + params.T) / 2  # impose symmetry along the diagonal
    return filter_project_fn(params[..., None], beta=beta)


# define the design region bounding box
design_region = td.Box(
    center=(0, 0, h / 2),
    size=(l, l, h),
)


def get_design_region(params: np.ndarray, beta: float) -> td.Structure:
    """Get design region structure as a function of optimization parameters."""
    density = get_density(params, beta=beta)
    eps_data = rescale(density, 1, n_GaP**2)
    return td.Structure.from_permittivity_array(eps_data=eps_data, geometry=design_region)

Next we create another function that takes the design parameters and create three simulations with the dipole polarized in the \(x\), \(y\), and \(z\) directions.

def make_sims(params: np.ndarray, beta: float) -> dict[str, td.Simulation]:
    """Get simulation as a function of optimization parameters."""

    # set up the design region structure
    design_region = get_design_region(params, beta=beta)

    # create a simulation with x-oriented dipole
    sim_x = td.Simulation(
        size=(l + buffer * 2, l + buffer * 2, monitor_z + buffer * 2),
        run_time=run_time,
        structures=[substrate, design_region],
        sources=[pd_x],
        monitors=[field_monitor],
        grid_spec=td.GridSpec.auto(
            min_steps_per_wvl=min_steps_per_wvl,
        ),
    )

    # create three simulations with different dipole orientations
    sims = {
        "x": sim_x,
        "y": sim_x.updated_copy(sources=[pd_y]),
        "z": sim_x.updated_copy(sources=[pd_z]),
    }

    return sims

To ensure the simulation setup is correct, we can create and visualize an initial simulation. The initial condition we are going to use is a uniform design region with permittivity \(\varepsilon=\frac{1+\varepsilon_{GaP}}{2}\)

# number of pixels in one direction
n_pixels = int(l / pixel_size)

# initial design parameters
params0 = np.zeros((n_pixels, n_pixels)) + 0.5

# simulations with the initial design parameters
sims0 = make_sims(params0, beta=50.0)

# viasualize a simulation
sims0["x"].plot(y=0)
plt.show()

Baseline Simulations

As a baseline, we will first simulate the scenario with no extractor. This result will be compared to that with the optimized extractor later on.

# create simulations without the extractor structure
sims_ref = {axis: sims0[axis].updated_copy(structures=[substrate]) for axis in ["x", "y", "z"]}

# run the simulations
ref_results = web.run_async(simulations=sims_ref, path_dir="data")

17:22:53 CET Started working on Batch containing 3 tasks.

17:22:57 CET Maximum FlexCredit cost: 0.075 for the whole batch.

             Use 'Batch.real_cost()' to get the billed FlexCredit cost after the
             Batch has completed.

17:23:09 CET Batch complete.

To facilitate post-processing, we will develop a function to extract the upward radiation fluxes across all simulations. When applied to the baseline simulation data, this function reveals that the flux for the dipole in the vertical direction is an order of magnitude smaller than in the horizontal directions. This occurs because most emission from a vertically oriented dipole undergoes total internal reflection at the diamond-air interface, significantly limiting the extraction efficiency.

def extract_fluxes(results: td.web.BatchData) -> list:
    """Postprocessing function to extract the flux data and print them."""
    fluxes = [results[axis]["field"].flux.data[0] for axis in ["x", "y", "z"]]

    for axis, flux in zip(["x", "y", "z"], fluxes):
        print(f"E{axis} dipole flux is {flux:.3f}.")

    return fluxes


flux_ref = extract_fluxes(ref_results)

Ex dipole flux is 74.819.
Ey dipole flux is 74.819.
Ez dipole flux is 4.566.

Objective Function

Finally, before we start the optimization, we define the objective function as the minimum upward flux among the three simulations. Note that we don’t apply a simply minimum since this discontinuous function will break the gradient tracking. Instead, we use a soft minimum function given by

\[-\tau [ln(e^{-f_x/\tau}) + ln(e^{-f_y/\tau}) + ln(e^{-f_z/\tau})],\]

where \(f_i\) denotes the upward radiation flux in the \(i=[x, y, z]\) direction, and \(\tau\) is the parameter that controls the “softness” of the soft minimum function. The smaller \(\tau\) we use, the function is closer to the true minimum.

from tidy3d.plugins.autograd import smooth_min


def min_upward_flux(params: np.ndarray, beta: float) -> float:
    """Objective function for the inverse design"""
    sims = make_sims(params, beta=beta)
    batch_results = web.run_async(simulations=sims, path_dir="data", verbose=False)

    def flux_data(dim: str) -> np.ndarray:
        return np.mean(batch_results[dim]["field"].flux.data)

    fluxes = np.array([flux_data("x"), flux_data("y"), flux_data("z")])

    return smooth_min(fluxes)

# function to calculate the objective function value and its gradient with respect to the design parameters.
val_grad_fn = ag.value_and_grad(min_upward_flux)

Optimization Loop

Now we are ready to start the optimization loop. We will use the adam optimizer. In each iteration, the beta value for the tanh projection is increased on an exponential schedule to help gradually binarize the design, slowly initially an then more strongly as the optimization progresses. In each iteration, we output the objective function value as well as plot the design region density.

# hyperparameters
num_steps = 70
learning_rate = 0.1

# initialize adam optimizer with starting parameters
params = np.copy(params0)
optimizer = optax.adam(learning_rate=learning_rate)
opt_state = optimizer.init(params)

# store history
objective_history = []
params_history = [params]

# gradually increase the binarization strength
beta_min = 1
beta_max = 150
betas = np.logspace(np.log10(beta_min), np.log10(beta_max), num_steps)

for i, beta in enumerate(betas):
    # compute gradient and current objective function value
    density = get_density(params, beta)

    plt.subplots(figsize=(2, 2))
    plt.imshow(density, cmap="gray_r", vmin=0, vmax=1, origin="lower")
    plt.axis("off")
    plt.show()

    # compute objective function value and gradient
    value, gradient = val_grad_fn(params, beta=beta)

    # outputs
    print(f"step {i + 1}:")
    print(f" beta = {beta:.3e}")
    print(f" objective = {value:.3e}")

    # compute and apply updates to the optimizer based on gradient (-1 sign to maximize obj_fn)
    updates, opt_state = optimizer.update(-gradient, opt_state, params)
    params[:] = optax.apply_updates(params, updates)
    np.clip(params, 0, 1, out=params)

    # save history
    objective_history.append(value)
    params_history.append(params)

step 1:
 beta = 1.000e+00
 objective = 9.460e+00

step 2:
 beta = 1.075e+00
 objective = 5.272e+01

step 3:
 beta = 1.156e+00
 objective = 1.297e+02

step 4:
 beta = 1.243e+00
 objective = 2.179e+02

step 5:
 beta = 1.337e+00
 objective = 3.515e+02

step 6:
 beta = 1.438e+00
 objective = 4.430e+02

step 7:
 beta = 1.546e+00
 objective = 5.097e+02

step 8:
 beta = 1.663e+00
 objective = 5.225e+02

step 9:
 beta = 1.788e+00
 objective = 5.359e+02

step 10:
 beta = 1.922e+00
 objective = 5.694e+02

step 11:
 beta = 2.067e+00
 objective = 6.072e+02

step 12:
 beta = 2.223e+00
 objective = 6.193e+02

step 13:
 beta = 2.390e+00
 objective = 6.284e+02

step 14:
 beta = 2.570e+00
 objective = 6.451e+02

step 15:
 beta = 2.764e+00
 objective = 6.752e+02

step 16:
 beta = 2.972e+00
 objective = 7.063e+02

step 17:
 beta = 3.196e+00
 objective = 7.301e+02

step 18:
 beta = 3.437e+00
 objective = 7.082e+02

step 19:
 beta = 3.696e+00
 objective = 7.326e+02

step 20:
 beta = 3.974e+00
 objective = 7.371e+02

step 21:
 beta = 4.273e+00
 objective = 7.546e+02

step 22:
 beta = 4.595e+00
 objective = 7.599e+02

step 23:
 beta = 4.941e+00
 objective = 7.502e+02

step 24:
 beta = 5.313e+00
 objective = 7.399e+02

step 25:
 beta = 5.713e+00
 objective = 7.666e+02

step 26:
 beta = 6.144e+00
 objective = 7.781e+02

step 27:
 beta = 6.607e+00
 objective = 7.807e+02

step 28:
 beta = 7.104e+00
 objective = 7.468e+02

step 29:
 beta = 7.639e+00
 objective = 7.335e+02

step 30:
 beta = 8.215e+00
 objective = 7.851e+02

step 31:
 beta = 8.833e+00
 objective = 7.906e+02

step 32:
 beta = 9.499e+00
 objective = 8.208e+02

step 33:
 beta = 1.021e+01
 objective = 8.077e+02

step 34:
 beta = 1.098e+01
 objective = 8.075e+02

step 35:
 beta = 1.181e+01
 objective = 8.076e+02

step 36:
 beta = 1.270e+01
 objective = 8.229e+02

step 37:
 beta = 1.366e+01
 objective = 8.217e+02

step 38:
 beta = 1.469e+01
 objective = 8.405e+02

step 39:
 beta = 1.579e+01
 objective = 8.079e+02

step 40:
 beta = 1.698e+01
 objective = 8.306e+02

step 41:
 beta = 1.826e+01
 objective = 8.132e+02

step 42:
 beta = 1.964e+01
 objective = 8.455e+02

step 43:
 beta = 2.111e+01
 objective = 8.425e+02

step 44:
 beta = 2.270e+01
 objective = 8.529e+02

step 45:
 beta = 2.441e+01
 objective = 8.518e+02

step 46:
 beta = 2.625e+01
 objective = 8.644e+02

step 47:
 beta = 2.823e+01
 objective = 8.727e+02

step 48:
 beta = 3.036e+01
 objective = 8.749e+02

step 49:
 beta = 3.264e+01
 objective = 8.819e+02

step 50:
 beta = 3.510e+01
 objective = 8.784e+02

step 51:
 beta = 3.775e+01
 objective = 8.608e+02

step 52:
 beta = 4.059e+01
 objective = 8.508e+02

step 53:
 beta = 4.365e+01
 objective = 8.695e+02

19:07:35 CET WARNING: No connection: Retrying for 180 seconds.                  

step 54:
 beta = 4.693e+01
 objective = 8.774e+02

step 55:
 beta = 5.047e+01
 objective = 8.747e+02

step 56:
 beta = 5.427e+01
 objective = 8.818e+02

step 57:
 beta = 5.836e+01
 objective = 8.921e+02

step 58:
 beta = 6.275e+01
 objective = 8.901e+02

step 59:
 beta = 6.748e+01
 objective = 8.911e+02

step 60:
 beta = 7.256e+01
 objective = 8.963e+02

step 61:
 beta = 7.803e+01
 objective = 8.988e+02

step 62:
 beta = 8.391e+01
 objective = 9.024e+02

step 63:
 beta = 9.023e+01
 objective = 9.032e+02

step 64:
 beta = 9.702e+01
 objective = 8.989e+02

step 65:
 beta = 1.043e+02
 objective = 8.938e+02

step 66:
 beta = 1.122e+02
 objective = 8.996e+02

step 67:
 beta = 1.206e+02
 objective = 8.932e+02

step 68:
 beta = 1.297e+02
 objective = 9.012e+02

step 69:
 beta = 1.395e+02
 objective = 9.060e+02

step 70:
 beta = 1.500e+02
 objective = 9.123e+02

Optimization Result Analysis

After the optimization, we can plot the objective function as the function of iterations. A steady increase of the objective function observed.

plt.plot(objective_history, c="red", linewidth=2)
plt.xlabel("iterations")
plt.ylabel("objective function")
plt.show()

The optimum design is not necessarily the last iteration as the objective function can fluctuate a bit. We will select the optimum design that gives the highest objective function value, re-simulate it, and calculate the flux values for different dipole orientations.

# optimum design iteration index
opt_index = np.argmax(objective_history)

# optimum design parameters
params_final = params_history[opt_index]

# resimulate the optimized design
sims_opt = make_sims(params_final, beta=beta)
opt_results = web.run_async(simulations=sims_opt, path_dir="data")

# extract the fluxes
flux_opt = extract_fluxes(opt_results)

19:30:21 CET Started working on Batch containing 3 tasks.

19:30:25 CET Maximum FlexCredit cost: 0.156 for the whole batch.

             Use 'Batch.real_cost()' to get the billed FlexCredit cost after the
             Batch has completed.

19:30:48 CET Batch complete.

Ex dipole flux is 912.776.
Ey dipole flux is 912.776.
Ez dipole flux is 968.528.

It’s also helpful to calculate the relative flux enhancement compared to the baseline simulation without the extractor. The result indicates that the final design achieves an order-of-magnitude enhancement in upward flux for an in-plane dipole source and over two orders of magnitude for an out-of-plane dipole source, similar to what the authors of the referenced paper achieved.

for axis, flux_o, flux_r in zip(["x", "y", "z"], flux_opt, flux_ref):
    print(f"Relative E{axis} dipole flux enhancement is {(flux_o / flux_r):.3f}.")

Relative Ex dipole flux enhancement is 12.200.
Relative Ey dipole flux enhancement is 12.200.
Relative Ez dipole flux enhancement is 212.138.

It is important to note that the optimized design is likely not fully binarized, even with the application of a tanh projection using a large beta value. For fabrication, however, a fully binarized design is necessary. Here, we manually enforce binarization and test if that negatively affects the design performance.

In this case, the results are only marginally affected.

# fully binarize the permittivity
params_binarized = np.round(params_final)

# create simulations
sims_binarized = make_sims(params_binarized, beta=beta)

# run simulations
binarized_results = web.run_async(simulations=sims_binarized, path_dir="data")

# extract results
flux_binarized = extract_fluxes(binarized_results)

19:30:56 CET Started working on Batch containing 3 tasks.

19:31:01 CET Maximum FlexCredit cost: 0.156 for the whole batch.

             Use 'Batch.real_cost()' to get the billed FlexCredit cost after the
             Batch has completed.

19:31:15 CET Batch complete.

Ex dipole flux is 746.306.
Ey dipole flux is 746.306.
Ez dipole flux is 1044.504.

We would need to export a GDS file for the optimized design for fabrication. The GDS file generation will discretize the design into polygons. Depending on the algorithm, the discretization can potentially impact the design performance as well so it’s good to test it.

First we export a GDS file.

# gds file export
gds_path = "./misc/inv_des_diamond_light_extractor.gds"
sims_binarized["x"].to_gds_file(
    fname=gds_path,
    z=h / 2,
    permittivity_threshold=(n_GaP**2 + 1) / 2,
    frequency=freq0,
)

Then we will import the GDS file using gdstk.

# load the gds file we just created
lib_loaded = gdstk.read_gds(gds_path)

# create a cell dictionary with all the cells in the file
all_cells = {c.name: c for c in lib_loaded.cells}
print("Cell names: " + ", ".join(all_cells.keys()))

Cell names: MAIN

cell_loaded = all_cells["MAIN"]

print(cell_loaded)

Cell 'MAIN' with 24 polygons, 0 flexpaths, 0 robustpaths, 0 references, and 0 labels

We can visualize the layout.

fig, ax = plt.subplots()

for polygon in cell_loaded.polygons:
    x, y = zip(*polygon.points)
    ax.fill(x, y, edgecolor="black", fill=True, c="blue")

ax.set_aspect("equal")
plt.xlabel("x (μm)")
plt.ylabel("y (μm)")
plt.show()

Create the extractor structure from the gds geometry. Here we don’t apply any sidewall angles but if the fabrication is known to produce a certain sidewall angle, it can be tested here.

# create geometry from the gds layout
extractor_geo = td.Geometry.from_gds(
    cell_loaded, gds_layer=0, gds_dtype=0, axis=2, slab_bounds=(0, h)
)

# create the extractor structure
extractor = td.Structure(geometry=extractor_geo, medium=GaP)

Simulate the imported extractor and check its performance. Compared to the design before GDS export, we do see a small decrease in the performance but it’s quite minimal.

# create simulations with the loaded gds geometry
sims_loaded = {
    axis: sims_binarized[axis].updated_copy(structures=[substrate, extractor])
    for axis in ["x", "y", "z"]
}

# run the simulations
loaded_results = web.run_async(simulations=sims_loaded, path_dir="data")

# extract results
flux_loaded = extract_fluxes(loaded_results)

for axis, flux_l, flux_r in zip(["x", "y", "z"], flux_loaded, flux_ref):
    print(f"Relative E{axis} dipole flux enhancement is {(flux_l / flux_r):.3f}.")

19:31:23 CET Started working on Batch containing 3 tasks.

19:31:28 CET Maximum FlexCredit cost: 0.156 for the whole batch.

             Use 'Batch.real_cost()' to get the billed FlexCredit cost after the
             Batch has completed.

19:31:49 CET Batch complete.

Ex dipole flux is 766.001.
Ey dipole flux is 766.470.
Ez dipole flux is 971.992.
Relative Ex dipole flux enhancement is 10.238.
Relative Ey dipole flux enhancement is 10.244.
Relative Ez dipole flux enhancement is 212.897.

Lastly we can plot the simulation in 3D to better visualize the optimized extractor geometry.

sims_loaded["x"].plot_3d()

Final Remarks

In this notebook, we focus exclusively on a single frequency. In the referenced paper, optimization is also conducted at a single frequency, while the extractor’s performance is evaluated across a broad wavelength range. This can be achieved by specifying frequency sampling points in the FieldMonitor.

The enhancement of upward radiation flux in the optimized design can arise from two mechanisms: (1) the structure redirects the dipole radiation toward the upward direction, and (2) the total radiation from the dipole is increased due to the Purcell effect. Although we do not examine the Purcell factor in this notebook, it can be analyzed as demonstrated in the tutorial provided here. Based on the findings of the paper, the Purcell effect is not very significant in our case since the dipole is deeply buried (100 nm) under the surface.

In this example, we only used the filter to enforce the fabrication constraint. For better fabrication constraint compliance, one can also consider adding a fabrication penalty term in the objective function such as the one demonstrated in this example. Before fabrication, we should also closely examine the small features in the design to ensure fabricability.