Emitters and fields in the time domain, rebuilt from scratch.

The project and where each part is addressed

The task wording is from the call and the QUEST project description; the split into modules is mine.

Checks against independent results

Built for this study

• A nodal DGTD solver (Legendre–Gauss–Lobatto elements of any order, upwind flux, low-storage RK4) with auxiliary differential equations for any sum of Drude–Lorentz poles.
• Emitters as sheets of N-level density matrices on element interfaces, coupled to the field through the jump condition of the numerical flux.
• Transfer-matrix references, including the exact complex pole of an emitter in a layered structure.
• Mie theory with a complex-frequency pole search for the quasinormal modes of a dispersive sphere.
• The Lifshitz–McLachlan Casimir–Polder integral, with both power-law limits as unit tests; tests run in CI on every push.

How each step would grow in the PhD

• Solver: 3D DGTD on unstructured meshes for real silicon nanostructures, with the same ADE and emitter coupling.
• Emitters: density matrices with dephasing and incoherent decay (Lindblad terms), many-level molecules, and emitters at several positions.
• Beyond semiclassics: the semiclassical model gives the Purcell rate in the weak-excitation limit but no vacuum fluctuations; quantised treatments for single emitters.
• Resonances: normalised quasinormal modes of 3D silicon resonators, and their use as a reduced model for the emitter dynamics.
• Atom–surface forces: Casimir–Polder near nanostructured silicon, where the planar formula no longer applies.

Gold: the two-band model of Ref. [1]

Discretising the interband integral into Lorentz poles

A Gauss–Legendre rule with $N_G$ nodes turns the integral into $N_G$ Lorentz poles, i.e. $N_G$ ADEs. Few poles leave unphysical peaks below 600 nm, as in Fig. 2 of [1]. With the columns $\omega_g$ and $\gamma$ of Table 1 in arXiv v1 exchanged, the model reproduces the peaks of Fig. 3 of [1] (); as printed it does not. I have not checked the journal version.

Silicon: one Lorentz model cannot cover both regions

Fitted from 250 nm, the Lorentz tails of the UV critical points leave a spurious absorption κ = at 1278 nm, where silicon is transparent. The near-infrared fit (700–1450 nm) brings it down to , with |Δn|/n ≤ . Modules C and D use the near-infrared model and module E the broadband one, whose static limit ε(0) = is close to the known 11.7.

Method and result

• Gold fit: 49 points of Johnson & Christy, residuals relative to $|n+i\kappa|$, Monte Carlo start and least squares. $\hbar\omega_p$ = eV and $\hbar\omega_g$ = eV; χ² = , against with the published parameters.
• Discretise, then refit (Sec. 4 of [1]): with $N_G$ = 3 poles the refitted model reaches χ² = .
• ADE check: integrating the ADEs for a short pulse and Fourier transforming gives back χ(ω) to (gold, 4 ADEs) and (silicon).

Reflection and transmission: DGTD against the transfer matrix

Convergence

Error of T(λ) for a 300 nm slab, 700–1450 nm. It falls exponentially with $N$ to for a constant ε; with the silicon ADEs it levels off near , a floor I have not traced. Against $h$ the measured slopes are for $N$ = 1–4, faster than $N+1$, as expected from the superconvergent dispersion error of DG [3].

Rabi flopping (nonlinear check)

A dilute two-level sheet in vacuum, driven by a resonant pulse of area θ: the final population follows $\sin^2(\theta/2)$ to . The solver keeps the counter-rotating terms that the $\sin^2$ law drops; the deviation shrinks for a longer pulse of the same area.

Purcell effect in the silicon cavity

What the three views show

• Position: decay rate of a weakly excited two-level sheet, fitted from $\rho_{ee}(t)$, against the transfer-matrix Purcell factor (cavity Q = , peak factor relative to bulk silicon in this 1D geometry).
• Frequency: a three-level V emitter with one transition at the cavity mode and one at 1.12 ω₀, inside the mirror stop band: the first decays times faster than in bulk silicon, the second times, both from the same run.
• Beyond Purcell: the Purcell formula assumes Γ ≪ κ. For Γ/κ up to the DGTD rate exceeds it by up to ; the exact complex pole of the same linear problem, $\omega_e^2-\omega^2-i\omega A Z_{\rm eff}(\omega)=0$, agrees with the DGTD to .
• In the weak-excitation limit the semiclassical decay equals that of a classical dipole, so the Purcell factor (an LDOS ratio) is the same as in the quantum treatment; the absolute rate of a single emitter is not.

Normal-mode splitting and saturation

With a dense sheet at the cavity centre the transmission peak splits. For a weak probe the Bloch equations are linear, and the DGTD matches the transfer matrix with a Lorentz sheet to . Stronger probes saturate the emitters and the splitting closes; this regime has no linear reference.

Purcell factor at nm from the surface

At 852 nm: tangential (exact) against from the MD pole alone; radial . A radial dipole does not couple to the MD mode at all. Adding the quadrupole poles makes the tangential estimate worse, not better, and the pole terms diverge at long wavelength where the true coefficients vanish as a power of the size parameter: a pole term alone lacks the non-resonant background of each coefficient.

Quasinormal modes of the sphere

Purcell factor at 852 nm against the gap. The error of the MD pole in the tangential F − 1 is : one resonance works close to the surface and fails as F approaches 1.

From van der Waals to Casimir–Polder

Near the surface $U \to -C_3/z^3$; far from it, $-3\hbar c\,\alpha(0)/(32\pi^2\varepsilon_0 z^4)$ for an ideal mirror. Both limits are unit tests of the integral. For the ideal mirror the potential is halfway between them at nm.

Result and comparison

• The two D lines give α(0) = a.u., of the measured 401.0 a.u. [9]; the core electrons are left out.
• Ref. [7] computes the same coefficients with relativistic coupled-cluster polarisabilities and tabulated optical data. Its valence parts agree with mine within ; its totals, with the core, are larger (4.84 a.u. for the ideal mirror).
• Silicon gives of the ideal-mirror C₃. Near nanostructured silicon the planar formula no longer holds; that is where a DGTD calculation of the Green's tensor, as in [10], is needed.

Simplifications

• The solver is 1D (normal incidence, emitters as sheets); 3D nanostructures are treated with Mie theory, not with DGTD.
• Emitters are semiclassical (no vacuum fluctuations, no dephasing); the decay comes from the field the emitters radiate themselves.
• Emitter strengths are chosen for affordable run times, not for a single G-centre; rates are compared in ratios, which do not depend on that choice.
• SiO₂ has a constant index 1.447; silicon in modules C–D keeps a small spurious absorption from the Lorentz model.
• The Cs polarisability has the D lines only; the pole approximation of module D keeps no background term.

Reproduce

python -m venv .venv && .venv/bin/pip install -e . pytest
.venv/bin/python scripts/run_moduleA.py   # materials (1 min)
.venv/bin/python scripts/run_moduleB.py   # DGTD validation (1 min)
OMP_NUM_THREADS=1 .venv/bin/python scripts/run_moduleC.py   # emitters (~15 min, 10 cores)
.venv/bin/python scripts/run_moduleD.py   # Mie, QNMs
.venv/bin/python scripts/run_moduleE.py   # Casimir-Polder
.venv/bin/python scripts/plots.py && .venv/bin/python scripts/build_site.py
.venv/bin/pytest

Publication figures (PDF)