Joshua Mornhinweg, Harvard John A. Paulson School of Engineering and Applied Sciences

Cavity quantum electrodynamics enters the deep-strong coupling regime when the vacuum Rabi frequency $\Omega_\mathrm{R}$ of a light-matter system becomes comparable to, or larger than, the bare transition frequency $\omega_\mathrm{c}$ of the matter excitation. Reaching $\Omega_\mathrm{R}/\omega_\mathrm{c} \gtrsim 1$ requires squeezing the electromagnetic vacuum field into a mode volume far below $(\lambda/2)^3$, because the single-photon field amplitude scales as $E_\mathrm{vac} \propto 1/\sqrt{V}$.

In the terahertz (THz) range this is often done with a planar metallic metasurface evaporated directly onto a substrate. The metasurface acts as an array of sub-wavelength LC resonators that concentrate the THz vacuum field into the near-field volume just below the metal, where it overlaps with our matter system of choice.

This notebook reproduces the bare THz resonator of

J. Mornhinweg, L. K. Diebel, M. Halbhuber, M. Prager, J. Riepl, T. Inzenhofer, D. Bougeard, R. Huber, and C. Lange, Mode-multiplexing deep-strong light-matter coupling, Nature Communications 15, 1847 (2024).

The bare resonator that we can compute here supports several well-separated modes - a fundamental LC resonance at $\nu_1 \approx 0.5$ THz and a higher-order mode at $\nu_2 \approx 2$ THz - which couple simultaneously to a comb of magnetoplasmon modes. In particular, this resonator geometry features an extremely compact unit cell. Rather than isolating individual resonators, the design uses an inverted (complementary, Babinet-type) layout in which the metal is continuous and the aperture is patterned. Regions of low current density are removed and adjacent current paths of opposite phase are merged, shrinking the unit cell from a conventional 60 x 60 μm to only 30 x 32.5 μm.

**Note: Resonances will differ slightly between versions. The presented notebook was run on Tidy3D 2.9.3

What this notebook computes¶

We model a single unit cell of the metasurface with periodic boundaries, illuminated by a normally incident, broadband THz plane wave, and extract

• the transmittance, reflectance and absorptance spectra from 0.1 to 6 THz, from which the resonator mode frequencies are read off, and
• the electric near-field distributions of the first two modes just below the metal - the quantity that determines the coupling strength to the matter system.

Model simplifications. As in the near-field analysis of the paper, we simulate the bare resonator on an undoped substrate: there are no quantum wells, no two-dimensional electron gas and no external magnetic field, so no polariton splitting appears here. GaAs is treated as a lossless, dispersionless dielectric with $\varepsilon = 12.9$, which is accurate across the THz band. The gold film is 200 nm thick and is described by a surface-impedance model at the literature conductivity of gold. Finally, the array is taken to be infinite and the illumination to be a plane wave, whereas the experiment uses a finite array and a focused THz beam.

Simulation setup¶

Frequency range¶

The resonator modes of interest lie between roughly 0.5 and 2 THz, but the metasurface is designed to couple to a magnetoplasmon comb spanning several octaves, so we record the spectra over the full 0.1-6 THz window. All lengths in Tidy3D are in micrometers and all frequencies in Hz.

Materials¶

Two materials are involved. The substrate is semi-insulating GaAs, which is essentially transparent and dispersion-free below the reststrahlen band, so a constant relative permittivity $\varepsilon = 12.9$ is an excellent approximation throughout the THz range.

The metasurface is gold, and here the choice of model matters more than it might seem. A Drude metal at THz frequencies has an enormous permittivity - $|\varepsilon| \sim 10^6$, i.e. a refractive index above 2000 - so the wavelength inside the gold is only a fraction of a micron. An automatic mesher will faithfully try to resolve that, and the grid explodes: the volumetric Drude version of this simulation needs a near-uniform 0.1 um grid and roughly $7 \times 10^6$ cells, almost all of them spent inside a metal whose interior carries no useful physics.

The physics of such a large $|\varepsilon|$ is simply that the field penetrates one skin depth and stops. That is what a surface impedance boundary condition (SIBC) encodes, and is implemented in LossyMetalMedium: the metal volume is replaced by a surface impedance $Z_s = \sqrt{i \mu_0 \omega / \sigma}$ applied at its boundary, so the mesh never has to enter the gold.

Three things to note about the arguments:

• conductivity is in S/um, not S/m. Literature gold has $\sigma = 4.1 \times 10^7$ S/m, hence 41.0.
• frequency_range is required - the impedance is fitted with a pole-residue model over that band, and construction raises a validation error without it.
• thickness switches on a one-dimensional transmission-line correction for a film that is not many skin depths thick. At 200 nm and $\delta \approx 100$ nm it is a small effect (it changes $Z_s$ by about 3% at 0.57 THz and under 1% at 2.15 THz), but it is the physically correct model for a film this thin.

SIBC is accurate when the skin depth is much smaller than the feature size. Here $\delta \approx 100$ nm against 2.5 um slots - a factor of 25 - so we are comfortably inside its range of validity.

Unit-cell geometry¶

The metasurface is an inverted structure: a continuous gold film perforated by a butterfly-shaped aperture. It is easiest to construct the aperture first and then subtract it from an unbroken film, which is exactly how the design is conceived - the complementary counterpart of four LC resonators sharing their inductive loops.

The aperture has the full symmetry of the unit cell about both mirror planes, which can later also be used to reduce the simulation load.

Source and monitors¶

The metasurface is excited from above by a broadband plane wave at normal incidence. The capacitive gap runs along $x$, so the mode is driven by a field polarized along $y$; we set pol_angle to $90^\circ$ accordingly.

Along $z$ the domain is one central wavelength high, with the GaAs filling the lower half and extending through the PML so that the substrate is effectively semi-infinite and no back-surface echo appears in the spectrum.

Four monitors record the response:

• T - a flux monitor inside the GaAs, below the metasurface, giving the transmittance,
• R - a flux monitor above the source, giving the reflectance (a plane-wave source radiates in one direction only, so everything crossing this plane is reflected light),
• field_xy - the in-plane electric field 0.5 um below the metal, i.e. at the example depth of the matter system, at the two resonance frequencies,
• field_xz - a vertical cut through the cell showing how far the near field reaches into the substrate.

Boundaries, mesh and symmetry¶

The array is periodic in $x$ and $y$ and absorbing in $z$. Because the excitation is a normally incident plane wave, simple periodic boundaries are sufficient - Bloch boundaries would only be needed at oblique incidence.

With the metal reduced to a surface impedance, the mesh is driven entirely by the dielectric problem. We define dl_min as a guard rail against runaway refinement. Here we have two places where resolution is required:

• In plane, the narrow features are the 2.5 μm slots and bars. min_steps_per_wvl = 25 plus a corner refinement gives roughly 0.5 μm cells.
• In $z$, the near field below the metal is what couples to our matter excitation. It decays over a few microns. Left to itself AutoGrid puts ~2.6 μm cells at the surface, which would smear that decay into a single cell and make the near-field maps meaningless. A MeshOverrideStructure pins the step to dz_near through the top z_near microns of the substrate.

Finally, with $\mathbf{E} \parallel \hat{y}$ the field has a definite parity about both mirror planes: $E_y$ is even under $x \to -x$ and odd under $y \to -y$, giving symmetry=(1, -1, 0). This cuts the computational domain - and therefore the cost - by a factor of four.

Sanity checks before submitting¶

Running the simulation¶

Results¶

Transmittance, reflectance and absorptance¶

Tidy3D normalizes flux monitors by the power spectrum of the source, so the recorded fluxes are directly the transmittance $T$ and reflectance $R$; the absorptance follows as $A = 1 - R - T$ and, since GaAs is lossless here, it measures ohmic dissipation in the gold alone.

For an aperture (Babinet-complementary) metasurface the resonances appear as transmission peaks rather than dips: on resonance the slot radiates efficiently into the substrate, while off resonance the continuous metal film reflects almost everything.

We also plot $\sqrt{T}$, the field transmission amplitude, since that is the quantity measured directly by THz time-domain spectroscopy, where the transmitted electric field $E(t)$ is sampled electro-optically.

The two lowest peaks are the modes that carry the coupling in the experiment; compare them against the measured $\nu_1 = 0.52$ THz and $\nu_2 = 1.95$ THz. Some offset is expected from the idealizations listed at the top - nominal rather than fabricated dimensions.

Near-field distributions of the two modes¶

The coupling strength is set not by the transmittance but by how tightly the vacuum field is confined where the matter system is located. The map below shows $|\mathbf{E}|$ in a plane 0.5 um below the gold at each of the two resonance frequencies.

What to look for: an LC resonance concentrates $|\mathbf{E}|$ across the central capacitive gap, with the surrounding loop carrying the circulating current that provides the inductance, so $\nu_1$ should appear as a single lobe straddling the gap. The higher-order mode $\nu_2$ should carry additional nodes along the arms and therefore vary on a shorter length scale.

A vertical cut shows how deep the field reaches into the substrate.

The confinement depth read off this cut is the design constraint that matters for the experiment. The mode is bound to the metal surface and decays into the GaAs over a scale set by the slot dimensions - tens of micrometers at most - rather than by the free-space wavelength of several hundred micrometers. That mismatch is exactly the sub-wavelength confinement that makes $\Omega_\mathrm{R}/\omega_\mathrm{c} > 1$ attainable.