28  Metamaterialien und Doppelbrechung

1968 stellte der sowjetische Physiker Victor Veselago eine Frage, die sich auf einer Briefmarke formulieren lässt: Was wäre, wenn ein Material gleichzeitig negatives \(\varepsilon\) und negatives \(\mu\) hätte? Seine Antwort auf dem Papier war verblüffend — Wellen liefen darin rückwärts, Brechung würde auf die falsche Seite knicken, eine flache Platte würde zur Linse. Und dann passierte: nichts. 32 Jahre lang blieb die Arbeit eine Kuriosität, denn kein bekanntes Material der Welt hat beides zugleich. Erst im Jahr 2000 baute die Gruppe um David Smith in San Diego ein „Material”, das es doch kann — zusammengesetzt aus Kupferdrähten und kleinen Ring-Schwingkreisen, beides Zentimeter groß. Der Trick: Wenn die Bausteine viel kleiner sind als die Wellenlänge, fragt die Welle nicht, woraus ein Medium besteht. Sie spürt nur gemittelte Antworten — ein effektives \(\varepsilon\), ein effektives \(\mu\). Solche maßgeschneiderten Kunst-Medien heißen Metamaterialien.

Das ist die Gegenrichtung zu Kapitel 26: Dort waren die Bausteine vergleichbar mit der Wellenlänge (\(a \sim \lambda\)), und alles drehte sich um Interferenz — Bandlücken, Bragg-Spiegel. Hier sind die Bausteine viel kleiner (\(a \ll \lambda\)), und die Welle sieht ein homogenes Medium mit Eigenschaften, die kein Naturstoff liefert. In diesem Kapitel bauen wir Veselagos Medium in Meep nach und messen alle drei Vorhersagen. Danach drehen wir den Spieß um und schauen auf Materialien, die zwei verschiedene Brechungsindizes gleichzeitig haben — je nachdem, wie das Licht polarisiert ist: Doppelbrechung, vom Kalkspat-Doppelbild bis zur λ/4-Platte, die aus linearem Licht das zirkulare macht, das wir in Kapitel 4 versprochen haben.

Lernziele

Nach diesem Kapitel kannst du …

  1. … erklären, warum \(\varepsilon < 0\) allein eine Welle sperrt (Kapitel 16), \(\varepsilon < 0\) und \(\mu < 0\) sie aber wieder laufen lassen — und warum dann \(n = -\sqrt{\varepsilon\mu}\) die richtige Wahl ist (Energie vorwärts, Phase rückwärts),
  2. … das historische Smith-Experiment nachbauen: Drähte allein sperren, Ringe allein sperren, beide zusammen öffnen ein Durchlassfenster — der Fingerabdruck des negativen Index,
  3. … negative Brechung und die Veselago-Flachlinse simulieren und weißt, woran reale Metamaterialien kranken (Verluste wohnen an der Resonanz, Kapitel 15),
  4. anisotrope Materialien (\(\varepsilon\) als Tensor) in Meep ansetzen und Doppelbrechung messen: zwei Polarisationen, zwei Brechungswinkel, zwei Indizes,
  5. … das Kalkspat-Doppelbild erklären (Walk-off: der außerordentliche Strahl wandert seitlich, obwohl er senkrecht einfällt) und eine λ/4-Platte dimensionieren, die linear polarisiertes Licht zirkular macht.

28.1 Verbotene Vorzeichen?

Fangen wir mit dem an, was wir schon wissen. Die Wellengleichung aus Kapitel 4 verknüpft die Wellenzahl mit den Materialkonstanten:

\[ k^2 = \omega^2 \varepsilon_0\mu_0\,\varepsilon\mu \quad\Longrightarrow\quad n^2 = \varepsilon\mu . \]

Ist nur \(\varepsilon\) negativ, wird \(k^2 < 0\) — die Wellenzahl ist imaginär, die Welle klingt exponentiell ab, statt zu laufen. Genau das haben wir zweimal gemessen: an der Ionosphären-Klippe in Kapitel 16 (Drude-Medium unter der Plasmafrequenz) und im Rohr unter dem Cutoff in Kapitel 17. Ein Medium mit \(\varepsilon < 0\) ist ein Spiegel, kein Leiter.

Aber jetzt Veselagos Schachzug: Sind \(\varepsilon\) und \(\mu\) negativ, ist ihr Produkt wieder positiv — \(k^2 > 0\), die Welle läuft! Und hier sitzt die Falle, die sich als nächstes aufdrängt:

WarnungNaheliegende Vermutung: „n = √((−1)·(−1)) = +1 — also verhält sich das Medium wie Vakuum“

Warum sie naheliegt: \(n^2 = \varepsilon\mu = (-1)(-1) = 1\), und die Schulregel sagt: Wurzel ziehen, fertig, \(n = 1\). Dann wäre das Doppel-Minus-Medium von Vakuum nicht zu unterscheiden — Veselagos Frage wäre langweilig.

Was stattdessen stimmt: Die Wurzel hat zwei Zweige, und Maxwell wählt für uns. Schreibe die beiden Rotationsgleichungen für eine ebene Welle (Kapitel 4): \(\vec{k} \times \vec{E} = \omega \mu_0\mu \vec{H}\) und \(\vec{k} \times \vec{H} = -\omega \varepsilon_0\varepsilon \vec{E}\). Bei positiven Materialkonstanten bilden \(\vec{E}\), \(\vec{H}\), \(\vec{k}\) ein rechtshändiges Dreibein — und der Energiefluss \(\vec{S} = \vec{E} \times \vec{H}\) zeigt entlang \(\vec{k}\). Drehst du die Vorzeichen von \(\varepsilon\) und \(\mu\) um, klappen beide Gleichungen um: Das Dreibein wird linkshändig, \(\vec{S}\) zeigt jetzt entgegen \(\vec{k}\). Die Energie fließt von der Quelle weg (das muss sie), also zeigt die Phasenausbreitung zur Quelle hin. Genau das bedeutet \(n = -1\): Snellius gilt weiter, aber mit negativem Winkel — und die Phase läuft rückwärts. Veselago nannte solche Stoffe darum „linkshändige Medien”.

Eine Sorge müssen wir noch ausräumen, bevor wir bauen: Darf es ein Medium mit \(\varepsilon = \mu = -1\) überhaupt geben? Die Energiedichte aus Kapitel 18, \(u = \tfrac{1}{2}\varepsilon_0\varepsilon E^2 + \tfrac{1}{2}\mu_0\mu H^2\), wäre damit negativ — ein Feld, das Energie erzeugt, wenn man es anfacht. Das wäre tatsächlich Unsinn. Die Auflösung: Diese Energieformel gilt nur für Konstanten ohne Frequenzabhängigkeit. In einem dispersiven Medium gilt stattdessen \(u = \tfrac{1}{2}\varepsilon_0\frac{\mathrm{d}(\omega\varepsilon)}{\mathrm{d}\omega}E^2 + \dots\) — und diese Ableitungen sind bei unseren Drude- und Lorentz-Modellen positiv, auch wo \(\varepsilon\) und \(\mu\) selbst negativ sind. Merksatz: Negative Materialkonstanten gibt es nur mit Dispersion. Ein \(n = -1\) bei allen Frequenzen wäre unphysikalisch; ein \(n = -1\) bei einer Frequenz ist erlaubt. Warum Kausalität und Dispersion so eng verheiratet sind, ist das Thema von Kapitel 29 (Kramers-Kronig).

28.2 Die Zutaten: Drähte und Ringe

Woher nehmen, wenn die Natur es nicht liefert? Smiths Antwort von 2000 besteht aus zwei Bausteinen, die wir beide schon kennen:

Baustein 1 — das Drahtgitter (macht \(\varepsilon < 0\)). Ein Gitter aus dünnen Metalldrähten, Abstand viel kleiner als die Wellenlänge, wirkt auf eine Welle mit \(\vec{E}\) parallel zu den Drähten wie ein verdünntes Plasma: Die Elektronen können den Drähten entlang schwingen, aber die Induktivität der dünnen Drähte (Kapitel 22) macht sie träge. Das Ergebnis ist genau die Drude-Antwort aus Kapitel 16,

\[ \varepsilon(f) = 1 - \frac{f_p^2}{f^2 + \mathrm{i}\gamma f}, \]

nur dass die „Plasmafrequenz” \(f_p\) jetzt keine Materialkonstante ist, sondern eine Stellschraube: Drahtabstand und Drahtdicke legen sie fest. Unterhalb von \(f_p\) ist \(\varepsilon < 0\) — die Ionosphäre aus dem Baumarkt. (In Übung 28.3 baust du dieses Gitter aus echten PEC-Stäbchen und misst seine Klippe.)

Baustein 2 — der Split-Ring (macht \(\mu < 0\)). Ein aufgeschnittener Metallring ist ein Schwingkreis im Sinn von Kapitel 19: Der Ring ist die Induktivität \(L\), der Schlitz der Kondensator \(C\). Das Magnetfeld der Welle durchsetzt den Ring und treibt einen Kreisstrom — ein getriebener Resonator (Kapitel 15). Und jetzt kommt der Phasensprung an der Resonanz ins Spiel, den wir dort gemessen haben: Knapp oberhalb der Ringresonanz \(f_0\) antwortet der Strom gegenphasig — der Ring baut ein Magnetfeld auf, das dem antreibenden Feld entgegensteht, und zwar stärker, als das Vakuum mitmacht. Die gemittelte Antwort ist ein Lorentz-Term im \(\mu\):

\[ \mu(f) = 1 + \frac{\sigma_R\, f_0^2}{f_0^2 - f^2 - \mathrm{i}\gamma f}, \]

der in einem Band knapp über \(f_0\) negativ wird. Anders als beim Draht (Drude: träge Elektronen, \(\varepsilon < 0\) für alles unter \(f_p\)) ist \(\mu < 0\) also ein Resonanz-Phänomen — es lebt nur in einem schmalen Band, und es wohnt zwangsläufig dort, wo auch die Verluste wohnen (Kapitel 15: an der Resonanz wird absorbiert). Das wird uns noch beschäftigen.

Wir setzen beide Modelle in Meep an. Neu ist nur eines: Dispersion im \(\mu\) schreibt man als H_susceptibilities — das exakte Gegenstück zu den E_susceptibilities aus Kapitel 16.

Zur Wahl der Zahlen: Wir arbeiten bei der Betriebsfrequenz \(f_{\text{op}} = 1\) (also \(\lambda = 1\); Skaleninvarianz wie immer seit Kapitel 11 — dieselbe Rechnung gilt für 10 GHz wie für 300 THz). Die Drude-Stärke wählen wir so, dass \(\varepsilon(1) = -1\) wird, die Ringresonanz legen wir auf \(f_0 = 0{,}8\) mit einer Stärke, die \(\mu(1) = -1\) liefert. Kleine Verluste \(\gamma = 0{,}01\) sind Absicht: Sie machen das Modell ehrlich (echte Metamaterialien sind verlustig) — und sie halten die Simulation stabil.

import numpy as np
import matplotlib.pyplot as plt
import meep as mp

mp.verbosity(0)

F_OP = 1.0                      # Betriebsfrequenz (lambda = 1)
GAMMA = 1e-2                    # Verluste der Bausteine
F0_RING, S_RING = 0.8, 1.125    # Ringresonanz und -stärke

# Drähte: Drude-epsilon (Kapitel 16). Meeps Konvention:
# eps(f) = 1 − sigma·f_ref²/(f² + i·gamma·f) — mit f_ref = F_OP und
# sigma = 2 wird eps(F_OP) = 1 − 2 = −1, exakt.
drude_e = mp.DrudeSusceptibility(frequency=F_OP, gamma=GAMMA, sigma=2.0)

# Ringe: Lorentz-mu (Kapitel-19-Schwingkreis als Materialantwort):
# mu(f) = 1 + sigma·f0²/(f0² − f² − i·gamma·f)
lorentz_h = mp.LorentzianSusceptibility(frequency=F0_RING, gamma=GAMMA,
                                        sigma=S_RING)

MAT_DRAEHTE = mp.Medium(epsilon=1, E_susceptibilities=[drude_e])
MAT_RINGE = mp.Medium(epsilon=1, mu=1, H_susceptibilities=[lorentz_h])
MAT_NIM = mp.Medium(epsilon=1, E_susceptibilities=[drude_e],
                    mu=1, H_susceptibilities=[lorentz_h])


def eps_von(f):
    """Drude-Antwort der Drähte (komplex)."""
    return 1.0 - 2.0 * F_OP**2 / (f**2 + 1j * GAMMA * f)


def mu_von(f):
    """Lorentz-Antwort der Ringe (komplex)."""
    return 1.0 + S_RING * F0_RING**2 / (F0_RING**2 - f**2 - 1j * GAMMA * f)


print(f"eps({F_OP}) = {eps_von(F_OP):.4f}")
print(f"mu({F_OP})  = {mu_von(F_OP):.4f}")
Using MPI version 4.1, 1 processes
eps(1.0) = -0.9998+0.0200j
mu(1.0)  = -0.9985+0.0555j

Beide Antworten stehen bei \(f = 1\) auf \(-1\) (plus kleine Imaginärteile — die Verluste). Bevor wir messen, lohnt der Blick auf die ganzen Kurven, denn sie legen die Bühne fest:

fs = np.linspace(0.3, 1.7, 600)
F_P = np.sqrt(2) * F_OP                          # eps-Nulldurchgang
F_MU_OBEN = F0_RING * np.sqrt(1 + S_RING)        # mu-Nulldurchgang

fig, (a1, a2) = plt.subplots(1, 2, figsize=(9.6, 3.6), sharex=True)
a1.plot(fs, eps_von(fs).real, "C0-", lw=1.6)
a1.axhline(0, color="k", lw=0.6, ls="--")
a1.axvline(F_P, color="C0", lw=0.8, ls=":")
a1.annotate(r"$f_p$", (F_P, 1.5), color="C0")
a1.set_ylim(-6, 3)
a1.set_xlabel("Frequenz f")
a1.set_ylabel(r"$\mathrm{Re}\,\varepsilon$")
a1.set_title("Drähte: Drude")
a2.plot(fs, mu_von(fs).real, "C3-", lw=1.6)
a2.axhline(0, color="k", lw=0.6, ls="--")
a2.axvspan(F0_RING, F_MU_OBEN, color="gold", alpha=0.25)
a2.annotate(r"$\mu < 0$", (0.93, 2.0), color="C3")
a2.set_ylim(-6, 6)
a2.set_xlabel("Frequenz f")
a2.set_ylabel(r"$\mathrm{Re}\,\mu$")
a2.set_title("Ringe: Lorentz (Resonanz bei 0,8)")
for a in (a1, a2):
    a.grid(alpha=0.3)
fig.tight_layout()
plt.show()
Abbildung 28.1: Die beiden Zutaten über der Frequenz. Links die Drude-Antwort der Drähte: ε < 0 für alles unterhalb von f_p = 1,41. Rechts die Lorentz-Antwort der Ringe: µ < 0 nur im Resonanzband 0,8 bis 1,17 (gold). Nur im Überlapp (beide Kurven unter der Nulllinie, gestrichelt) sind beide Vorzeichen negativ — dort, und nur dort, kann die Welle wieder laufen.

Lies die Karte: Die Drähte sperren (allein) das ganze Band bis \(f_p = 1{,}41\). Die Ringe sperren (allein) ihr Resonanzband \(0{,}8\) bis \(1{,}17\). Das Ringband liegt vollständig im Drahtband — und genau dort sind beide Vorzeichen gleichzeitig negativ.

28.3 Das Smith-Experiment: zwei Sperren öffnen ein Fenster

Versuchsaufbau. Frage: Was lässt ein Slab aus Drahtgitter, Ringgitter oder beidem zusammen durch? Bühne: eine 1D-Zelle (dimensions=1, wie der λ/4-Stapel in Kapitel 26), in der Mitte ein Slab der Dicke \(d = 2\) aus dem jeweiligen Material, beidseitig 3 Einheiten Luft und 2 Einheiten PML. Anregung: ein breitbandiger Gauß-Puls (Mitte \(f = 1\), Breite \(1{,}4\) — deckt \(0{,}3\) bis \(1{,}7\) ab). Messgröße: die Transmission \(T(f)\) als Fluss-Verhältnis mit Referenzlauf, exakt das Rezept aus Kapitel 12. Erfolgskriterium: Die drei Kurven müssen auf der parameterfreien Transfer-Matrix-Linie liegen, die wir gleich danach aus Kapitel 26 übernehmen.

Und jetzt du — der Vorhersagepunkt dieses Kapitels:

HinweisVorhersage (PRIMM)

Die Drähte allein sperren bei \(f = 1\) praktisch alles (\(\varepsilon = -1\), \(\mu = +1\) → evaneszent, Kapitel 16). Die Ringe allein sperren bei \(f = 1\) ebenfalls praktisch alles (\(\varepsilon = +1\), \(\mu = -1\) → dieselbe Mathematik, nur mit getauschten Rollen). Was misst du, wenn beide Gitter zusammen im Slab stehen — zwei Sperren hintereinander? (a) Noch dichter — Sperren addieren sich. (b) Genau so dicht. (c) Der Slab wird durchlässig. Entscheide dich, bevor du weiterliest — und notiere in einem Satz, warum.

D_SLAB = 2.0                # Slabdicke
AUFL_1D = 64
DPML_1D, PAD_1D = 2.0, 3.0
SZ_1D = 2 * (DPML_1D + PAD_1D) + D_SLAB


def sweep_1d(material, fcen=1.0, df=1.4, nfreq=281):
    """T(f) durch den Slab: Lauf mit Material / Leerlauf (Kap. 12)."""
    quelle = [mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ex,
                        center=mp.Vector3(0, 0, -SZ_1D / 2 + DPML_1D + 1.0))]
    z_mess = SZ_1D / 2 - DPML_1D - 1.0
    leistungen = []
    for geometrie in ([], [mp.Block(size=mp.Vector3(mp.inf, mp.inf, D_SLAB),
                                    center=mp.Vector3(), material=material)]):
        sim = mp.Simulation(cell_size=mp.Vector3(0, 0, SZ_1D),
                            resolution=AUFL_1D, dimensions=1, Courant=0.25,
                            boundary_layers=[mp.PML(DPML_1D)],
                            geometry=geometrie, sources=quelle)
        tr = sim.add_flux(fcen, df, nfreq,
                          mp.FluxRegion(center=mp.Vector3(0, 0, z_mess)))
        sim.run(until_after_sources=mp.stop_when_fields_decayed(
            50, mp.Ex, mp.Vector3(0, 0, z_mess), 1e-9))
        leistungen.append(np.array(mp.get_fluxes(tr)))
    fs = np.array(mp.get_flux_freqs(tr))
    return fs, leistungen[1] / leistungen[0]


fs_1d, T_draehte = sweep_1d(MAT_DRAEHTE)
_, T_ringe = sweep_1d(MAT_RINGE)
_, T_beide = sweep_1d(MAT_NIM)

i_op = np.argmin(np.abs(fs_1d - F_OP))
print(f"T(f = 1): Drähte allein {T_draehte[i_op]:.5f}, "
      f"Ringe allein {T_ringe[i_op]:.5f}, beide {T_beide[i_op]:.4f}")
T(f = 1): Drähte allein 0.00000, Ringe allein 0.00000, beide 0.3829

Antwort (c): Zwei Sperren hintereinander sind durchlässig. Jede allein blockt auf \(10^{-4}\), beide zusammen lassen 38 % durch. Wer hier (a) getippt hat, hat in Alltagslogik gedacht — zwei Mauern sind dicker als eine. Aber die Mauern sperren aus entgegengesetzten Gründen (\(\varepsilon < 0\) bzw. \(\mu < 0\)), und zusammen heben sich die Gründe auf: \(k^2 = \omega^2\varepsilon\mu > 0\), die Welle läuft. Genau mit dieser Messung — Fenster nur, wenn beide Gitter montiert sind — hat Smith 2000 den negativen Index nachgewiesen.

Jetzt die Theorielinie. Die Transfer-Matrix aus Kapitel 26 kann das alles parameterfrei, wir müssen ihr nur beibringen, mit komplexem \(n\) und komplexem Wellenwiderstand \(Z = \sqrt{\mu/\varepsilon}\) zu rechnen (bisher war stillschweigend \(\mu = 1\)). Dabei fällt die Vorzeichenwahl, über die wir oben geredet haben — und sie folgt zwei physikalischen Geboten:

def schicht_T(fs, eps_funktion, mu_funktion, d):
    """Transmission einer Schicht — Kap.-26-Matrix mit komplexem n, Z.

    Die Wurzel-Zweige wählt die Physik:
      Im(n) >= 0  — die Welle klingt in Laufrichtung ab (Kausalität:
                    ein passives Medium verstärkt nicht);
      Re(Z) >= 0  — die Energie fließt vorwärts.
    Sind eps und mu beide negativ, landet Re(n) damit von selbst
    im Negativen: n = −1 ist keine Setzung, sondern Konsequenz.
    """
    Ts = []
    for f in fs:
        eps, mu = eps_funktion(f), mu_funktion(f)
        n = np.sqrt(eps * mu)
        if n.imag < 0:
            n = -n                       # Zweig mit Im(n) >= 0
        Z = np.sqrt(mu / eps)
        if Z.real < 0:
            Z = -Z                       # Zweig mit Re(Z) >= 0
        delta = 2 * np.pi * f * n * d    # Laufphase (e^{-i delta}, Kap. 26)
        M = np.array([[np.cos(delta), -1j * Z * np.sin(delta)],
                      [-1j * np.sin(delta) / Z, np.cos(delta)]])
        t = 2.0 / ((M[0, 0] + M[0, 1]) + (M[1, 0] + M[1, 1]))
        Ts.append(abs(t)**2)
    return np.array(Ts)


eins = lambda f: 1.0 + 0j
T_tmm_draehte = schicht_T(fs_1d, eps_von, eins, D_SLAB)
T_tmm_ringe = schicht_T(fs_1d, eins, mu_von, D_SLAB)
T_tmm_beide = schicht_T(fs_1d, eps_von, mu_von, D_SLAB)

maske = (fs_1d > 0.4) & (fs_1d < 1.6)
for name, T, Tt in [("Drähte", T_draehte, T_tmm_draehte),
                    ("Ringe", T_ringe, T_tmm_ringe),
                    ("beide", T_beide, T_tmm_beide)]:
    rms = np.sqrt(np.mean((T[maske] - Tt[maske])**2))
    print(f"rms(Meep − Matrix), {name}: {rms:.4f}")
rms(Meep − Matrix), Drähte: 0.0040
rms(Meep − Matrix), Ringe: 0.0008
rms(Meep − Matrix), beide: 0.0054
fig, ax = plt.subplots(figsize=(7.2, 4.2))
for T, Tt, farbe, name in [(T_draehte, T_tmm_draehte, "C0", "Drähte allein"),
                           (T_ringe, T_tmm_ringe, "C2", "Ringe allein"),
                           (T_beide, T_tmm_beide, "C3", "beide zusammen")]:
    ax.plot(fs_1d, T, ".", ms=3.5, color=farbe, label=f"Meep: {name}")
    ax.plot(fs_1d, Tt, "-", lw=1.1, color=farbe, alpha=0.55)
ax.axvspan(F0_RING, F_MU_OBEN, color="gold", alpha=0.2,
           label=r"$\mu < 0$ (Ring-Band)")
ax.axvline(F_P, color="C0", lw=0.8, ls=":")
ax.set_xlim(0.3, 1.7)
ax.set_ylim(-0.03, 1.05)
ax.set_xlabel("Frequenz f")
ax.set_ylabel("Transmission T")
ax.set_title("Zwei Sperren öffnen ein Fenster")
ax.legend(fontsize=8, loc="upper left")
ax.grid(alpha=0.3)
plt.show()
Abbildung 28.2: Das Smith-Experiment in 1D: Transmission durch den Slab für Drähte allein (blau), Ringe allein (grün) und beide zusammen (rot); Punkte sind Meep, Linien die Transfer-Matrix mit komplexem n und Z. Das goldene Band markiert µ < 0. Jedes Gitter allein sperrt dort — beide zusammen öffnen genau dort das Durchlassfenster des negativen Index.

Meep liegt auf der Matrix-Linie (rms unter 0,006 über den ganzen Sweep). Zwei Dinge verdient die Kurve noch:

Erstens, die weichen Flanken. Unterhalb von \(f_0 = 0{,}8\) ist nur \(\varepsilon\) negativ — verboten — und doch sickert bei \(0{,}5\) sichtbar etwas durch. Das ist kein Fehler, sondern Kapitel 17: Der Slab ist nur \(2\lambda\) dick, und evaneszente Wellen tunneln durch dünne verbotene Zonen. Verbotene Bänder haben bei endlicher Dicke immer weiche Ränder.

Zweitens, das Fenster ist nicht klar, sondern milchig: \(T = 0{,}38\) statt \(1{,}0\), obwohl die Anpassung bei \(f = 1\) perfekt ist (\(Z = \sqrt{\mu/\varepsilon} = 1\) — wie Vakuum, nichts wird reflektiert!). Der Fehlbetrag ist reine Absorption, und wir können das quantitativ festnageln: Aus dem komplexen \(n\) folgt die Lambert-Beer-Dämpfung (Kapitel 6) über die Slabdicke:

n_op = np.sqrt(eps_von(F_OP) * mu_von(F_OP))
if n_op.imag < 0:
    n_op = -n_op
# Leistung faellt wie e^(-2 Im(k) d) mit k = 2 pi f n
T_beer = np.exp(-2 * (2 * np.pi * F_OP * n_op.imag) * D_SLAB)
print(f"n(1) = {n_op:.4f}")
print(f"Lambert-Beer aus Im(n): T = {T_beer:.3f}  "
      f"(gemessen: {T_beide[i_op]:.3f})")
n(1) = -0.9993+0.0378j
Lambert-Beer aus Im(n): T = 0.387  (gemessen: 0.383)

\(\mathrm{Im}(n) = 0{,}038\) erklärt das Fenster vollständig. Und dieses \(\mathrm{Im}(n)\) ist kein Schönheitsfehler unserer Parameter, sondern Systemzwang: \(\mu < 0\) gibt es nur nahe der Ringresonanz (Lorentz!), und an der Resonanz wohnen die Verluste (Kapitel 15). Reale Metamaterialien im Optischen verlieren pro Wellenlänge oft mehr, als unsere \(\gamma = 0{,}01\) es tun — das ist bis heute ihr wundester Punkt.

In der Spielwiese kannst du das Fenster selbst verschieben — alles reine NumPy-Transfer-Matrix von oben. Drei Aufträge, erst vorhersagen: (1) Schiebe die Ringresonanz F0_RING auf \(0{,}5\) — wo sitzt das Fenster jetzt? (2) Verdopple die Verluste GAMMA_R — was passiert mit der Fenster-Höhe, was mit den verbotenen Flanken? (3) Mache die Ringe stärker (S_RING = 3.0) — warum wird das Fenster breiter?

28.4 Die Phase läuft rückwärts

Bis hierher könnte ein Skeptiker sagen: „Schönes Fenster — aber woher weißt du, dass \(n\) darin negativ ist und nicht einfach irgendein gewöhnlicher Wert?” Gute Frage. Der Fingerabdruck von \(n < 0\) ist die rückwärts laufende Phase, und die können wir direkt messen.

Versuchsaufbau. Frage: In welche Richtung wächst die Phase der Welle vor, im und hinter dem Slab? Bühne: die 1D-Zelle von oben, etwas verlängert, Slab \(d = 2\) aus dem Doppel-Material in der Mitte. Anregung: Dauerstrich bei \(f = 1\) mit force_complex_fields=True — dann trägt jedes \(E_x(z)\) seine Phase mit sich (Kapitel 15). Messgröße: das Phasenprofil \(\varphi(z) = \arg E_x(z)\), mit np.unwrap zu einer stetigen Kurve zusammengesetzt (das Werkzeug aus Kapitel 27, das dort die Umlaufmode gezählt hat); seine Steigung ist die lokale Wellenzahl. Erfolgskriterium: Steigung \(+k_0\) in der Luft — und wenn \(n = -1\) stimmt, exakt \(-k_0\) im Slab.

SZ_CW = 2 * (DPML_1D + 5.5) + D_SLAB
sim = mp.Simulation(cell_size=mp.Vector3(0, 0, SZ_CW), resolution=AUFL_1D,
                    dimensions=1, Courant=0.25, force_complex_fields=True,
                    boundary_layers=[mp.PML(DPML_1D)],
                    geometry=[mp.Block(size=mp.Vector3(mp.inf, mp.inf, D_SLAB),
                                       center=mp.Vector3(), material=MAT_NIM)],
                    sources=[mp.Source(mp.ContinuousSource(F_OP),
                                       component=mp.Ex,
                                       center=mp.Vector3(0, 0, -SZ_CW / 2
                                                         + DPML_1D + 1.0))])
sim.run(until=400)

ex_cw = sim.get_array(component=mp.Ex, center=mp.Vector3(),
                      size=mp.Vector3(0, 0, SZ_CW - 2 * DPML_1D))
zs_cw = np.linspace(-(SZ_CW / 2 - DPML_1D), SZ_CW / 2 - DPML_1D, len(ex_cw))
phase_cw = np.unwrap(np.angle(ex_cw))
k0 = 2 * np.pi * F_OP

steigungen = {}
for name, z1, z2 in [("davor", -4.5, -1.5), ("im Slab", -0.8, 0.8),
                     ("dahinter", 1.5, 4.5)]:
    m = (zs_cw > z1) & (zs_cw < z2)
    steigungen[name] = np.polyfit(zs_cw[m], phase_cw[m], 1)[0] / k0
    print(f"Phasensteigung {name:9s}: {steigungen[name]:+.4f} · k0")
Phasensteigung davor    : +1.0015 · k0
Phasensteigung im Slab  : -1.0007 · k0
Phasensteigung dahinter : +1.0016 · k0
fig, ax = plt.subplots(figsize=(7.2, 4.0))
ax.axvspan(-D_SLAB / 2, D_SLAB / 2, color="0.85")
ax.plot(zs_cw, phase_cw / k0, "C0-", lw=1.6, label="gemessen")
for z_start, ph_start, stg in [(-4.5, None, 1), (1.5, None, 1)]:
    m0 = np.argmin(np.abs(zs_cw - z_start))
    zz = np.linspace(z_start, z_start + 3, 10)
    ax.plot(zz, (phase_cw[m0] + (zz - z_start) * k0) / k0, "k--", lw=0.9,
            alpha=0.6)
ax.plot([], [], "k--", lw=0.9, alpha=0.6, label="Steigung +k0 (Luft)")
ax.set_xlabel("z (grau: der Slab)")
ax.set_ylabel(r"Phase $\varphi(z)\,/\,k_0$")
ax.set_title(f"Steigung im Slab: {steigungen['im Slab']:+.3f} · k0")
ax.legend(fontsize=9)
ax.grid(alpha=0.3)
plt.show()
Abbildung 28.3: Das Phasenprofil der Dauerstrich-Welle (f = 1) über den Ort. In der Luft steigt die Phase mit +k₀ (gestrichelte Referenzgeraden) — im Slab (grau) fällt sie mit exakt −k₀. Die Welle läuft als Ganzes nach rechts, aber ihre Phase läuft im Material nach links: n = −1, wörtlich gemessen. (Der kleine Knick ganz links liegt vor der Quelle bei z = −5,5: Dort läuft die Welle nach links, ihre Phase fällt also ebenfalls.)

Da steht es: \(+1{,}00\,k_0\) davor, \(-1{,}00\,k_0\) drin, \(+1{,}00\,k_0\) dahinter. Im Slab sind die Wellenberge rückwärts unterwegs.

Und die Energie? Die läuft vorwärts — sonst käme hinten nichts an. Wie schnell, sagt uns ein Puls. Seine Gruppengeschwindigkeit können wir sogar von Hand vorhersagen, mit der Gruppenindex-Formel aus Kapitel 16, \(n_g = n + f\,\mathrm{d}n/\mathrm{d}f\). Die Ableitung entzerren wir wie üblich Schritt für Schritt: Aus \(n^2 = \varepsilon\mu\) folgt mit der Produktregel

\[ 2n\,\frac{\mathrm{d}n}{\mathrm{d}f} = \varepsilon\,\frac{\mathrm{d}\mu}{\mathrm{d}f} + \mu\,\frac{\mathrm{d}\varepsilon}{\mathrm{d}f}. \]

Bei \(f = 1\) ist \(\varepsilon = \mu = n = -1\), und die beiden Steigungen liest man aus den Modellformeln ab: \(\mathrm{d}\varepsilon/\mathrm{d}f = 4\) (Drude) und \(\mathrm{d}\mu/\mathrm{d}f = 11{,}1\) (Lorentz). Also

\[ \frac{\mathrm{d}n}{\mathrm{d}f} = \frac{(-1)(11{,}1) + (-1)(4)}{2 \cdot (-1)} = +7{,}56, \qquad n_g = -1 + 7{,}56 = 6{,}56 . \]

Zwei Dinge stehen in dieser kleinen Rechnung: Erstens ist \(n_g\) positiv — die Energie läuft vorwärts, obwohl \(n\) negativ ist. Zweitens ist \(n_g\) groß: Der Puls kriecht mit \(c/6{,}56 \approx 0{,}15\,c\). Beides prüft die Messung — wir schicken einen schmalbandigen Puls durch einen dickeren Slab (\(d = 3\)) und stoppen die Zeit zwischen zwei Sonden:

D_DICK = 3.0
sz_p = 2 * (DPML_1D + 5.0) + D_DICK
sim = mp.Simulation(cell_size=mp.Vector3(0, 0, sz_p), resolution=AUFL_1D,
                    dimensions=1, Courant=0.25,
                    boundary_layers=[mp.PML(DPML_1D)],
                    geometry=[mp.Block(size=mp.Vector3(mp.inf, mp.inf, D_DICK),
                                       center=mp.Vector3(), material=MAT_NIM)],
                    sources=[mp.Source(mp.GaussianSource(F_OP, fwidth=0.08),
                                       component=mp.Ex,
                                       center=mp.Vector3(0, 0, -sz_p / 2
                                                         + DPML_1D + 0.5))])
vorn, hinten, zeiten = [], [], []


def messe_laufzeit(s):
    zeiten.append(s.meep_time())
    vorn.append(abs(s.get_field_point(
        mp.Ex, mp.Vector3(0, 0, -D_DICK / 2 - 1.0))))
    hinten.append(abs(s.get_field_point(
        mp.Ex, mp.Vector3(0, 0, D_DICK / 2 + 1.0))))


sim.run(mp.at_every(0.5, messe_laufzeit), until=120)
zeiten = np.array(zeiten)
dt_puls = zeiten[np.argmax(hinten)] - zeiten[np.argmax(vorn)]
v_g = D_DICK / (dt_puls - 2.0)     # 2,0 = Vakuumweg zwischen den Sonden
print(f"Laufzeit Sonde→Sonde: {dt_puls:.1f}; "
      f"v_g im Slab = {v_g:.3f} c (Vorhersage 1/6,56 = {1 / 6.56:.3f} c)")
Laufzeit Sonde→Sonde: 21.5; v_g im Slab = 0.154 c (Vorhersage 1/6,56 = 0.152 c)

\(v_g = 0{,}154\,c\) gemessen, \(0{,}152\,c\) vorhergesagt. Energie vorwärts (langsam), Phase rückwärts — beides gleichzeitig, beides quantitativ.

WarnungNaheliegende Vermutung: „Phase rückwärts — das ist doch Zeitumkehr und verletzt die Kausalität“

Warum sie naheliegt: Rückwärtslaufende Wellenberge sehen aus wie ein rückwärts abgespielter Film — und der verstößt gegen den Zeitpfeil aus Kapitel 20.

Was stattdessen stimmt: Transportiert wird nur, was Energie und Information trägt — und beides reist mit der Gruppe, vorwärts und (hier) gemächlich mit \(0{,}15\,c\). Die Phasenberge sind keine Boten; das haben wir schon in Kapitel 16 gemessen, als \(v_p = 1{,}67\,c\) über der Plasmafrequenz auch niemandem Überlicht beschert hat. Ein negatives \(v_p\) ist genauso harmlos wie ein überlichtschnelles. Der Film läuft vorwärts — nur das Riffelmuster darauf bewegt sich gegen den Strom, wie das Profil einer Bohrschnecke, die sich dreht: Das Muster wandert nach oben, das Bohrmehl nach unten.

Der Film macht aus dem gemessenen komplexen CW-Feld eine laufende Welle (der Zeiger-Trick \(\mathrm{Re}(E\, \mathrm{e}^{-\mathrm{i}\varphi})\) aus Kapitel 15). Bevor du startest: Außen laufen die Berge nach rechts — was tun sie im grauen Slab? Verfolge einen einzelnen Berg mit dem Finger.

Code der Animation (nur in der HTML-Fassung)
# von oben: zs_cw, ex_cw, D_SLAB (komplexes CW-Feld bei f = 1)
from matplotlib import animation
from IPython.display import HTML

phasen = np.linspace(0, 2 * np.pi, 16, endpoint=False)
fig, ax = plt.subplots(figsize=(7.2, 3.8))
y_max = 1.15 * np.abs(ex_cw).max()
ax.axvspan(-D_SLAB / 2, D_SLAB / 2, color="0.85")
ax.plot(zs_cw, np.abs(ex_cw), "k-", lw=0.8, alpha=0.45)
ax.plot(zs_cw, -np.abs(ex_cw), "k-", lw=0.8, alpha=0.45)
linie, = ax.plot(zs_cw, np.real(ex_cw), "C0-", lw=1.6)
ax.set_ylim(-y_max, y_max)
ax.set_xlabel("z (grau: der Slab mit n = −1)")
ax.set_ylabel("$E_x$")
titel = ax.set_title("", fontsize=10)
fig.subplots_adjust(top=0.86)


def bild(i):
    zeiger = np.exp(-1j * phasen[i])
    linie.set_ydata(np.real(ex_cw * zeiger))
    titel.set_text("Energie nach rechts (gemessen: 0,15 c) — "
                   "die Berge im Slab laufen nach links\n"
                   f"t = {phasen[i] / (2 * np.pi):.2f} Perioden")
    return [linie, titel]


anim = animation.FuncAnimation(fig, bild, frames=len(phasen), interval=150)
plt.close(fig)
HTML(anim.to_jshtml(default_mode="loop"))

28.5 Negative Brechung: der Knick auf die falsche Seite

Jetzt in die Fläche. Snellius (Kapitel 12) mit \(n = -1\) sagt: \(\sin\theta_2 = \sin\theta_1 / n = -\sin\theta_1\) — der gebrochene Strahl liegt auf der gleichen Seite der Normalen wie der einfallende. Das klingt harmlos und sieht verstörend aus. Wir messen es.

Versuchsaufbau. Frage: Auf welche Seite knickt ein schräger Strahl in einer Platte aus dem Doppel-Material — und wie groß ist der Versatz nach dem Durchgang? Bühne: eine 2D-Zelle (\(20 \times 12\), PML ringsum), darin horizontal eine Platte der Dicke \(t = 3\); einmal aus gewöhnlichem Glas (\(n = 1{,}5\)) als Kontrolle, einmal aus dem Material des Smith-Fensters. Anregung: ein Gauß-Bündel (GaussianBeamSource, Taille \(1{,}5\), Dauerstrich \(f = 1\)), das unter \(30°\) von oben einfällt. Messgröße: der Schwerpunkt des Strahls direkt über dem Eintritt und direkt unter dem Austritt; die Differenz ist der seitliche Versatz über die Platte. Erwartung: Glas bricht zur Normalen hin, \(\theta_2 = 19{,}5°\), Versatz \(+t\tan 19{,}5° = +1{,}06\); bei \(n = -1\) gilt \(\theta_2 = -30°\), Versatz \(-t\tan 30° = -1{,}73\) — anderes Vorzeichen, und betragsmäßig sogar größer als der Einfallswinkel-Versatz von Glas.

Ein Bau-Detail ist diesmal Pflichtlektüre — es hat diese Simulation drei Anläufe gekostet:

WichtigWerkstattkasten: das µ-dispersive Material darf die PML nicht berühren

Der erste Versuch dieser Simulation — Platte quer durch die ganze Zelle, wie wir es seit Kapitel 12 immer bauen — starb nach wenigen hundert Schritten mit Meeps Fehlermeldung „fields are NaN or Inf”. Die Detektivarbeit (Bausteine einzeln testen: nur-\(\varepsilon\) stabil, nur-\(\mu\) explodiert, homogene Zelle mit periodischen Rändern stabil) zeigte: Schuld war der Überlapp von Material und PML. Meeps FAQ nennt den Fall ausdrücklich — in Medien, die Rückwärtswellen tragen, dreht die PML-Streckung ihre Dämpfung in Verstärkung um; auch der einfachere mp.Absorber hilft hier nicht. In 1D war der Slab nie in der PML, darum fiel es erst jetzt auf. Die Regel für alle NIM-Läufe in diesem Kapitel: Blöcke endlich bauen, mit Luftpuffer zur PML — was am Plattenrand an Beugung entsteht, ist der ehrliche Preis (echte Proben sind auch endlich).

AUFL_2D = 32
MAT_GLAS = mp.Medium(index=1.5)


def brechung_slab(material, until=150):
    """Gauß-Bündel, 30° Einfall, Platte t = 3 — komplexes CW-Feld."""
    t, dpml = 3.0, 1.0
    sx, sy = 20.0, 12.0
    theta = np.radians(30)
    quelle = [mp.GaussianBeamSource(
        mp.ContinuousSource(F_OP),
        center=mp.Vector3(-2.5, sy / 2 - dpml - 0.3),
        size=mp.Vector3(8.0, 0), beam_x0=mp.Vector3(0, 0),
        beam_kdir=mp.Vector3(np.sin(theta), -np.cos(theta)),
        beam_w0=1.5, beam_E0=mp.Vector3(0, 0, 1))]
    block = mp.Block(size=mp.Vector3(sx - 2 * dpml - 2.0, t, mp.inf),
                     center=mp.Vector3(), material=material)
    sim = mp.Simulation(cell_size=mp.Vector3(sx, sy), resolution=AUFL_2D,
                        Courant=0.25, force_complex_fields=True,
                        boundary_layers=[mp.PML(dpml)],
                        geometry=[block], sources=quelle)
    sim.run(until=until)
    ez = sim.get_array(component=mp.Ez, center=mp.Vector3(),
                       size=mp.Vector3(sx - 2 * dpml, sy - 2 * dpml))
    return ez, sx - 2 * dpml, sy - 2 * dpml


def strahl_versatz(feld, bx, by, t=3.0):
    """Strahlschwerpunkt unter dem Austritt minus über dem Eintritt."""
    nx, ny = feld.shape
    xs = np.linspace(-bx / 2, bx / 2, nx)
    ys = np.linspace(-by / 2, by / 2, ny)
    intensitaet = np.abs(feld)**2

    def schwerpunkt(y0):
        iy = np.argmin(np.abs(ys - y0))
        zeile = intensitaet[:, iy]
        return np.sum(xs * zeile) / np.sum(zeile)

    return schwerpunkt(-t / 2 - 0.2) - schwerpunkt(+t / 2 + 0.2)


ez_glas, bx_b, by_b = brechung_slab(MAT_GLAS)
ez_nim, _, _ = brechung_slab(MAT_NIM)
v_glas = strahl_versatz(ez_glas, bx_b, by_b)
v_nim = strahl_versatz(ez_nim, bx_b, by_b)
print(f"Versatz Glas: {v_glas:+.3f}  (Theorie +1,06)")
print(f"Versatz NIM:  {v_nim:+.3f}  (Theorie −1,73)")
Versatz Glas: +1.092  (Theorie +1,06)
Versatz NIM:  -1.558  (Theorie −1,73)
fig, achsen = plt.subplots(1, 2, figsize=(10.5, 4.4), sharey=True)
for ax, feld, name, versatz in [(achsen[0], ez_glas, "Glas n = +1,5", v_glas),
                                (achsen[1], ez_nim, "Metamaterial n = −1",
                                 v_nim)]:
    ax.imshow(np.abs(feld).T, origin="lower", cmap="inferno",
              extent=[-bx_b / 2, bx_b / 2, -by_b / 2, by_b / 2],
              vmax=0.85 * np.abs(ez_glas).max())
    ax.axhline(-1.5, color="w", lw=0.7)
    ax.axhline(1.5, color="w", lw=0.7)
    ax.annotate("", xy=(-1.0, 1.5), xytext=(-2.5, 4.2),
                arrowprops=dict(color="w", arrowstyle="->", lw=1.2))
    ax.set_title(f"{name}:  Versatz {versatz:+.2f}")
    ax.set_xlabel("x")
achsen[0].set_ylabel("y")
fig.tight_layout()
plt.show()
Abbildung 28.4: Dasselbe Bündel, zwei Platten: links Glas (n = 1,5), rechts das Metamaterial (n = −1). Gezeigt ist |Ez| (Momentaufnahme des Dauerstrichs), die weißen Linien markieren die Platte, der Pfeil den einfallenden Strahl. Im Glas knickt der Strahl zur Normalen, im Metamaterial auf die andere Seite — der Austrittsversatz wechselt das Vorzeichen (gemessen +1,09 bzw. −1,56; die Theoriewerte +1,06 und −1,73 gelten für den unendlich breiten Strahl).

Im Glasbild läuft der Strahl wie seit Kapitel 12 gewohnt: hinein, zur Normalen geknickt, hinaus, Versatz nach rechts. Im Metamaterial-Bild knickt derselbe Strahl am Eintritt nach links zurück — er kreuzt die eigene Einfallsrichtung — und tritt um \(-1{,}56\) versetzt aus (der Theoriewert \(-1{,}73\) gilt für unendlich breite Strahlen; unser Bündel ist nur \(1{,}5\lambda\) schmal und verliert etwas an die Beugung an den Plattenkanten). Ein Vorzeichen als Messwert — genau wie es 2001 im Mikrowellenlabor gemessen wurde, am Prisma statt an der Platte.

28.6 Die Veselago-Flachlinse

Veselagos hübscheste Folgerung kommt ganz ohne Krümmung aus. Eine Sammellinse funktioniert, weil gekrümmtes Glas schräge Strahlen zur Achse zurückbiegt. Bei \(n = -1\) biegt schon die ebene Grenzfläche jeden Strahl auf die andere Seite der Normalen — eine flache Platte sammelt von selbst: Strahlen, die von einem Punkt vor der Platte auseinanderlaufen, werden am Eintritt zurückgefaltet, kreuzen sich in der Platte und nach dem Austritt noch einmal. Eine Platte der Dicke \(t\) bildet einen Punkt im Abstand \(d_1\) davor in einen Punkt im Abstand \(t - d_1\) dahinter ab. Keine Brennweite, keine Achse, kein Schleifen — jede Punktquelle vor der Platte bekommt ihr Bild dahinter.

HinweisVorhersage (PRIMM)

Unsere Platte ist \(t = 2\) dick, die Punktquelle sitzt \(d_1 = 1\) vor ihrer Vorderseite. Die Vorderseite liegt bei \(y = -1\), die Rückseite bei \(y = +1\). Bei welcher \(y\)-Koordinate erwartest du das Bild? Und: Was erwartest du für eine Quelle, die weiter als \(t\) vor der Platte sitzt (\(d_1 > t\))?

Versuchsaufbau. Frage: Wo bündelt die flache Platte das Licht der Punktquelle? Bühne: 2D-Zelle \(14 \times 10\), Platte \(t = 2\) (endlich breit, Luftpuffer zur PML — der Werkstattkasten gilt weiter). Anregung: Punktquelle (Dauerstrich \(f = 1\)) bei \(y = -2\), also \(d_1 = 1\) vor der Platte. Messgrößen: die Karte \(|E_z|\), das Intensitätsprofil auf der Achse hinter der Platte (Lage des Maximums) und die Halbwertsbreite des Fokus quer zur Achse. Erfolgskriterium: Bild bei \(y = +1 + (t - d_1) = +2\), Fokusbreite nahe der Beugungsgrenze \(\sim \lambda/2\).

T_LINSE, D1 = 2.0, 1.0


def flachlinse(material=None, until=200):
    dpml = 1.0
    sx, sy = 14.0, 10.0
    block = mp.Block(size=mp.Vector3(sx - 2 * dpml - 2.0, T_LINSE, mp.inf),
                     center=mp.Vector3(),
                     material=material if material else MAT_NIM)
    sim = mp.Simulation(cell_size=mp.Vector3(sx, sy), resolution=AUFL_2D,
                        Courant=0.25, force_complex_fields=True,
                        boundary_layers=[mp.PML(dpml)],
                        geometry=[block],
                        sources=[mp.Source(mp.ContinuousSource(F_OP),
                                           component=mp.Ez,
                                           center=mp.Vector3(0, -T_LINSE / 2
                                                             - D1))])
    sim.run(until=until)
    ez = sim.get_array(component=mp.Ez, center=mp.Vector3(),
                       size=mp.Vector3(sx - 2 * dpml, sy - 2 * dpml))
    return ez, sx - 2 * dpml, sy - 2 * dpml


ez_linse, bx_l, by_l = flachlinse()
nx_l, ny_l = ez_linse.shape
xs_l = np.linspace(-bx_l / 2, bx_l / 2, nx_l)
ys_l = np.linspace(-by_l / 2, by_l / 2, ny_l)
int_linse = np.abs(ez_linse)**2

# Fokus: Intensitätsmaximum auf der Achse, oberhalb der Platte
ix0 = np.argmin(np.abs(xs_l))
profil_achse = np.where(ys_l > 1.2, int_linse[ix0, :], 0)
iy_fokus = np.argmax(profil_achse)
y_fokus = ys_l[iy_fokus]

# Quer-Halbwertsbreite am Fokus (Randschutz wie in Kap. 27 gelernt)
quer = int_linse[:, iy_fokus] / int_linse[:, iy_fokus].max()
i_max = int(np.argmax(quer))
lo, hi = i_max, i_max
while lo > 0 and quer[lo - 1] >= 0.5:
    lo -= 1
while hi < len(quer) - 1 and quer[hi + 1] >= 0.5:
    hi += 1
fwhm_fokus = xs_l[hi] - xs_l[lo]
print(f"Fokus auf der Achse bei y = {y_fokus:+.3f}  (Theorie +2,0)")
print(f"Quer-FWHM des Fokus: {fwhm_fokus:.3f}  (lambda/2 = 0,5)")
Fokus auf der Achse bei y = +1.914  (Theorie +2,0)
Quer-FWHM des Fokus: 0.405  (lambda/2 = 0,5)
fig, ax = plt.subplots(figsize=(7.0, 5.0))
ax.imshow(np.abs(ez_linse).T, origin="lower", cmap="inferno",
          extent=[-bx_l / 2, bx_l / 2, -by_l / 2, by_l / 2],
          vmax=0.6 * np.abs(ez_linse).max())
ax.axhline(-T_LINSE / 2, color="w", lw=0.7)
ax.axhline(T_LINSE / 2, color="w", lw=0.7)
ax.plot(0, -T_LINSE / 2 - D1, "wo", ms=7, mfc="none", mew=1.4)
ax.plot(0, y_fokus, "w+", ms=12, mew=1.6)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("Eine flache Platte als Linse")
plt.show()
Abbildung 28.5: Die Flachlinse bei der Arbeit: |Ez|-Karte der Punktquelle (weißer Kreis bei y = −2) vor der Metamaterial-Platte (weiße Linien, t = 2). Das Licht wird in der Platte einmal und dahinter ein zweites Mal gebündelt; das weiße Kreuz markiert das gemessene Bild bei y = +1,91 ≈ t − d₁ hinter der Austrittsfläche — ohne jede Krümmung.

Das Bild sitzt bei \(y = +1{,}91\) — neun Hundertstel vor dem Ideal, der Tribut an die Verluste — und sein Fokus ist mit \(\mathrm{FWHM} = 0{,}41\) etwa beugungsbegrenzt scharf. In der Platte selbst siehst du das innere Bild: die erste Kreuzung der zurückgefalteten Strahlen, bei Tiefe \(d_1\) unter der Eintrittsfläche.

Eine Randnotiz mit großer Geschichte: Im Jahr 2000 zeigte John Pendry, dass die ideale, verlustfreie Veselago-Platte sogar die evaneszenten Felder der Quelle (Kapitel 12) wieder aufrichten würde — sie wäre eine „perfekte Linse” mit Auflösung unterhalb der Beugungsgrenze. Unser gemessenes FWHM von \(0{,}41 < \lambda/2\) deutet diesen Effekt bereits an. Warum daraus trotzdem kein Wunder-Mikroskop wurde, steht im Kleingedruckten — die Kurzfassung hast du in diesem Kapitel schon zweimal gelesen: Verluste wohnen an der Resonanz. Was aus den Verlusten folgt, misst du selbst in Übung 28.4.

28.7 Doppelbrechung: ein Material, zwei Indizes

Bisher haben wir an den Vorzeichen von \(\varepsilon\) gedreht. Die zweite Hälfte des Kapitels dreht an etwas anderem: an der Richtung. Alle unsere Materialien waren bisher isotrop — derselbe \(\varepsilon\)-Wert, egal wie das Feld orientiert ist. Aber viele Kristalle sind anisotrop: Ihre Atome sitzen in Reihen und Schichten, und ein \(\vec{E}\)-Feld längs der Schichten polarisiert sie leichter als eines quer dazu. Dann ist \(\varepsilon\) keine Zahl mehr, sondern ein Tensor — für uns ganz praktisch: drei verschiedene Werte für die drei Feldrichtungen,

\[ \varepsilon = \begin{pmatrix} \varepsilon_x & 0 & 0\\ 0 & \varepsilon_y & 0\\ 0 & 0 & \varepsilon_z \end{pmatrix}, \]

solange wir die Koordinatenachsen auf die Kristallachsen legen. Die Konsequenz steht sofort da: Eine Welle, deren \(\vec{E}\) in \(z\) zeigt, spürt \(\varepsilon_z\) und läuft mit \(n_z = \sqrt{\varepsilon_z}\); eine mit \(\vec{E}\) in \(x\) spürt \(\varepsilon_x\). Ein Material, zwei Brechungsindizes — sortiert nach Polarisation. Das heißt Doppelbrechung, und unser Modellkristall ist der Klassiker: Kalkspat (Calcit), mit \(n_o = 1{,}658\) („ordentlich”, für \(\vec{E}\) senkrecht zur besonderen Kristallachse) und \(n_e = 1{,}486\) („außerordentlich”, für \(\vec{E}\) parallel dazu) — eine der größten Doppelbrechungen unter den häufigen Mineralen.

Unsere 2D-Welt ist dafür wie geschaffen, denn sie hat von Natur aus zwei getrennte Polarisationen (Kapitel 13): die \(E_z\)-Welle (das Feld zeigt aus der Ebene) und die \(H_z\)-Welle (das Feld liegt in der Ebene). Legen wir die besondere Kristallachse in \(z\)-Richtung, dann gilt \(\varepsilon_x = \varepsilon_y = n_o^2\) und \(\varepsilon_z = n_e^2\) — die \(E_z\)-Welle ist der außerordentliche Strahl, die \(H_z\)-Welle der ordentliche.

Versuchsaufbau. Frage: Brechen sich die beiden Polarisationen am selben Kristall verschieden? Bühne: 2D-Zelle \(22 \times 14\), ein Kalkspat-Block der Dicke \(t = 6\) quer durch die Zelle (epsilon_diag — ein gewöhnliches, verlustfreies Material, das auch in der PML brav bleibt). Anregung: dasselbe Gauß-Bündel wie oben, Einfall \(45°\), einmal als \(E_z\)-, einmal als \(H_z\)-Welle. Messgröße: diesmal nicht der Strahlschwerpunkt, sondern der Wellenvektor selbst — wir schneiden die komplexe Feldkarte im Kristall aus und holen uns per 2D-FFT den dominanten \(\vec{k}\): seine Richtung ist der Brechungswinkel, sein Betrag (geteilt durch \(f\)) der Index. Eine Messung, zwei Ergebnisse. Erwartung (Snellius): \(\theta_e = \arcsin(\sin 45° / 1{,}486) = 28{,}4°\) und \(\theta_o = \arcsin(\sin 45° / 1{,}658) = 25{,}2°\).

Warum nicht wieder der Schwerpunkt? Werkstatt-Erfahrung aus dem Aufbau dieses Kapitels: Beim \(45°\)-Einfall ist die interne Reflexion der \(E_z\)-Welle (s-Polarisation, Kapitel 13) kräftig genug, dass ihr Interferenzmuster den Schwerpunkt verfälscht — die \(H_z\)-Welle läuft dagegen nahe am Brewster-Winkel und ist sauber. Die FFT trennt die Laufrichtungen von selbst und ist gegen solche Überlagerungen immun; das Hanning-Fenster und das Zero-Padding kennen wir sinngemäß aus den Spektren von Kapitel 10, nur eben in zwei Dimensionen.

N_O, N_E = 1.658, 1.486
MAT_KALKSPAT = mp.Medium(epsilon_diag=mp.Vector3(N_O**2, N_O**2, N_E**2))
TH_EIN = np.radians(45)
AUFL_DB = 48                       # feiner: wir messen |k| auf Prozente


def brechung_aniso(polarisation, until=150):
    """45°-Bündel auf den Kalkspat-Block; liefert die Feldkarte IM Kristall."""
    t, dpml = 6.0, 1.0
    sx, sy = 22.0, 14.0
    e0 = (mp.Vector3(0, 0, 1) if polarisation == "Ez"
          else mp.Vector3(np.cos(TH_EIN), np.sin(TH_EIN), 0))
    quelle = [mp.GaussianBeamSource(
        mp.ContinuousSource(F_OP),
        center=mp.Vector3(-5.0, sy / 2 - dpml - 0.3),
        size=mp.Vector3(8.0, 0), beam_x0=mp.Vector3(0, 0),
        beam_kdir=mp.Vector3(np.sin(TH_EIN), -np.cos(TH_EIN)),
        beam_w0=1.5, beam_E0=e0)]
    sim = mp.Simulation(cell_size=mp.Vector3(sx, sy), resolution=AUFL_DB,
                        Courant=0.25, force_complex_fields=True,
                        boundary_layers=[mp.PML(dpml)],
                        geometry=[mp.Block(size=mp.Vector3(mp.inf, t, mp.inf),
                                           center=mp.Vector3(),
                                           material=MAT_KALKSPAT)],
                        sources=quelle)
    sim.run(until=until)
    komponente = mp.Ez if polarisation == "Ez" else mp.Hz
    return sim.get_array(component=komponente, center=mp.Vector3(0, 0),
                         size=mp.Vector3(16.0, 5.0))


def k_dominant(feld, bx=16.0, by=5.0):
    """Dominanter k-Vektor einer komplexen Feldkarte per 2D-FFT.

    Zutaten: Hanning-Fenster gegen die harten Ausschnittsränder
    (Kap. 10), Zero-Padding ×4 für ein feines k-Raster, danach
    Schwerpunkt der 3×3-Umgebung des Spektral-Maximums.
    """
    nx, ny = feld.shape
    # np.hanning: das Glockenfenster aus der Kap.-10-Werkstatt, hier in
    # beiden Richtungen als äußeres Produkt über die Karte gelegt
    gefenstert = feld * np.hanning(nx)[:, None] * np.hanning(ny)[None, :]
    # np.fft.fft2 = DFT in beiden Achsen zugleich (Karte -> k-Raum);
    # np.fft.fftshift rückt die Nullfrequenz vom Rand in die Mitte,
    # damit kx/ky wie gewohnt von negativ nach positiv laufen
    spektrum = np.fft.fftshift(np.fft.fft2(gefenstert, s=(4 * nx, 4 * ny)))
    kx = np.fft.fftshift(np.fft.fftfreq(4 * nx, d=bx / nx))
    ky = np.fft.fftshift(np.fft.fftfreq(4 * ny, d=by / ny))
    # argmax liefert den Index im flachgeklopften Array —
    # np.unravel_index übersetzt ihn zurück in das (i, j)-Paar
    i, j = np.unravel_index(np.argmax(np.abs(spektrum)), spektrum.shape)
    umgebung = np.abs(spektrum[i - 1:i + 2, j - 1:j + 2])
    di = (umgebung * np.array([-1, 0, 1])[:, None]).sum() / umgebung.sum()
    dj = (umgebung * np.array([-1, 0, 1])[None, :]).sum() / umgebung.sum()
    return (kx[i] + di * (kx[1] - kx[0]), ky[j] + dj * (ky[1] - ky[0]))


felder_db = {}
for pol, n_soll in [("Ez", N_E), ("Hz", N_O)]:
    feld = brechung_aniso(pol)
    felder_db[pol] = feld
    kxp, kyp = k_dominant(feld)
    winkel = np.degrees(np.arctan2(abs(kxp), abs(kyp)))
    n_mess = np.hypot(kxp, kyp) / F_OP
    soll = np.degrees(np.arcsin(np.sin(TH_EIN) / n_soll))
    print(f"{pol}-Welle: Winkel {winkel:5.2f}° (Snellius {soll:5.2f}°),  "
          f"n = {n_mess:.3f} (Soll {n_soll})")
Ez-Welle: Winkel 28.39° (Snellius 28.41°),  n = 1.479 (Soll 1.486)
Hz-Welle: Winkel 25.59° (Snellius 25.24°),  n = 1.664 (Soll 1.658)
fig, achsen = plt.subplots(1, 2, figsize=(10.5, 3.6), sharey=True)
for ax, pol, name, n_soll in [(achsen[0], "Ez", "außerordentlich (n_e)", N_E),
                              (achsen[1], "Hz", "ordentlich (n_o)", N_O)]:
    feld = felder_db[pol]
    ax.imshow(np.real(feld).T, origin="lower", cmap="RdBu",
              extent=[-8, 8, -2.5, 2.5], aspect="auto")
    th = np.arcsin(np.sin(TH_EIN) / n_soll)
    zz = np.linspace(-2.5, 2.5, 10)
    versatz_strahl = -2.0          # Eintrittspunkt des Bündels (x bei y=2,5)
    ax.plot(versatz_strahl + (2.5 - zz) * np.tan(th), zz, "k--", lw=1.2)
    ax.set_title(f"{pol}-Welle: {name}")
    ax.set_xlabel("x")
achsen[0].set_ylabel("y")
fig.tight_layout()
plt.show()
Abbildung 28.6: Ein Kristall, zwei Brechungswinkel: die Feldkarte im Kalkspat (Ausschnitt im Block) für die Ez-Welle (links, außerordentlich — sie spürt n_e = 1,486) und die Hz-Welle (rechts, ordentlich — sie spürt n_o = 1,658). Die gestrichelten Linien zeigen die jeweilige Snellius-Richtung; die per 2D-FFT gemessenen Wellenvektoren treffen Winkel und Index auf unter ein Prozent.

Zwei Polarisationen, zwei Winkel, zwei Indizes — gemessen auf unter ein Prozent. Schickst du unpolarisiertes Licht (eine Mischung aus beiden, Kapitel 13) durch den Kristall, trennt es sich in zwei Strahlen. Das ist die halbe Erklärung des berühmten Kalkspat-Doppelbilds. Die andere Hälfte ist verblüffender — und kommt jetzt.

28.8 Walk-off: das Doppelbild bei senkrechtem Einfall

Wer einen Kalkspat-Kristall auf eine Buchseite legt, sieht die Schrift doppelt — auch wenn er senkrecht hindurchschaut. Senkrechter Einfall heißt aber \(\theta_1 = 0\), und Snellius liefert dann \(\theta_2 = 0\) für jeden Index. Woher kommt das zweite Bild?

Der Punkt: Auf der Buchseite liegt die Kristallachse im Allgemeinen schräg zur Blickrichtung. Dann ist der \(\varepsilon\)-Tensor in unseren Koordinaten nicht mehr diagonal — er bekommt Außerdiagonal-Einträge, und die verkoppeln die Feldrichtungen: \(D_x\) hängt dann auch von \(E_y\) ab. Für die außerordentliche Welle bedeutet das, dass Wellenfronten und Energiefluss nicht mehr parallel laufen. Die Fronten laufen weiter geradeaus (der \(\vec{k}\)-Vektor gehorcht der Grenzflächen-Logik aus Kapitel 12), aber der Poynting-Vektor (Kapitel 18!) kippt um den „Walk-off”-Winkel \(\rho\) zur Seite. Der ordentliche Strahl kennt das Problem nicht — sein \(\vec{E}\) zeigt aus der Ebene, für ihn ist der Kristall isotrop.

Für die Achse unter \(45°\) liefert die Kristalloptik (wir zitieren sie hier, statt sie herzuleiten)

\[ \tan\rho = \frac{n(45°)^2}{2} \left(\frac{1}{n_e^2} - \frac{1}{n_o^2}\right) \sin(2\cdot 45°), \qquad \frac{1}{n(45°)^2} = \frac{1/2}{n_o^2} + \frac{1/2}{n_e^2}, \]

was für Kalkspat \(\rho = 6{,}23°\) ergibt — über eine Kristalldicke von \(t = 6\) also \(6 \cdot \tan\rho = 0{,}65\) seitlicher Versatz.

Versuchsaufbau. Frage: Wandert der außerordentliche Strahl seitlich, obwohl er senkrecht einfällt — und der ordentliche nicht? Bühne: 2D-Zelle \(14 \times 14\), Kalkspat-Block \(t = 6\), dessen Achse um \(45°\) in der Ebene gedreht ist. Den gedrehten Tensor bauen wir von Hand: Die Drehung mischt \(\varepsilon_o\) und \(\varepsilon_e\) zu \(\varepsilon_{xx} = \varepsilon_{yy} = (\varepsilon_o + \varepsilon_e)/2\) und erzeugt den Koppelterm \(\varepsilon_{xy} = (\varepsilon_o - \varepsilon_e)/2\) — in Meep ist das ein epsilon_offdiag. Anregung: Gauß-Bündel senkrecht von oben, einmal \(E_z\) (ordentlich: spürt nur das ungedrehte \(\varepsilon_{zz} = \varepsilon_o\)), einmal \(H_z\) (außerordentlich: spürt den gedrehten Tensor in der Ebene). Messgröße: der seitliche Strahlversatz über den Block — jetzt darf es wieder der Schwerpunkt sein, denn bei senkrechtem Einfall ist die Reflexion winzig. Erfolgskriterium: \(E_z\)-Versatz \(= 0\), \(H_z\)-Versatz \(= 0{,}65\).

EPS_O, EPS_E = N_O**2, N_E**2
eps_mittel = (EPS_O + EPS_E) / 2
eps_koppel = (EPS_O - EPS_E) / 2
MAT_GEDREHT = mp.Medium(epsilon_diag=mp.Vector3(eps_mittel, eps_mittel,
                                                EPS_O),
                        epsilon_offdiag=mp.Vector3(eps_koppel, 0, 0))
T_WALK = 6.0
n_45 = 1 / np.sqrt(0.5 / N_O**2 + 0.5 / N_E**2)
tan_rho = (n_45**2 / 2) * (1 / N_E**2 - 1 / N_O**2)
print(f"Theorie: rho = {np.degrees(np.arctan(tan_rho)):.2f}°, "
      f"Versatz über t = {T_WALK}: {T_WALK * tan_rho:.3f}")


def walkoff(polarisation, until=150):
    dpml = 1.0
    sx, sy = 14.0, 14.0
    e0 = mp.Vector3(0, 0, 1) if polarisation == "Ez" else mp.Vector3(1, 0, 0)
    quelle = [mp.GaussianBeamSource(
        mp.ContinuousSource(F_OP),
        center=mp.Vector3(0, sy / 2 - dpml - 0.3),
        size=mp.Vector3(10.0, 0), beam_x0=mp.Vector3(0, 0),
        beam_kdir=mp.Vector3(0, -1), beam_w0=1.5, beam_E0=e0)]
    sim = mp.Simulation(cell_size=mp.Vector3(sx, sy), resolution=AUFL_2D,
                        Courant=0.25, force_complex_fields=True,
                        boundary_layers=[mp.PML(dpml)],
                        geometry=[mp.Block(size=mp.Vector3(mp.inf, T_WALK,
                                                           mp.inf),
                                           center=mp.Vector3(),
                                           material=MAT_GEDREHT)],
                        sources=quelle)
    sim.run(until=until)
    komponente = mp.Ez if polarisation == "Ez" else mp.Hz
    feld = sim.get_array(component=komponente, center=mp.Vector3(),
                         size=mp.Vector3(sx - 2 * dpml, sy - 2 * dpml))
    return feld, sx - 2 * dpml, sy - 2 * dpml


feld_o, bx_w, by_w = walkoff("Ez")
feld_e, _, _ = walkoff("Hz")
v_ord = strahl_versatz(feld_o, bx_w, by_w, t=T_WALK)
v_ausser = strahl_versatz(feld_e, bx_w, by_w, t=T_WALK)
print(f"Versatz ordentlich (Ez):        {v_ord:+.3f}")
print(f"Versatz außerordentlich (Hz):   {v_ausser:+.3f}")
Theorie: rho = 6.23°, Versatz über t = 6.0: 0.655
Versatz ordentlich (Ez):        -0.000
Versatz außerordentlich (Hz):   -0.653
fig, achsen = plt.subplots(1, 2, figsize=(10.0, 4.6), sharey=True)
for ax, feld, name, versatz in [(achsen[0], feld_o, "ordentlich (Ez)", v_ord),
                                (achsen[1], feld_e, "außerordentlich (Hz)",
                                 v_ausser)]:
    ax.imshow(np.abs(feld).T, origin="lower", cmap="inferno",
              extent=[-bx_w / 2, bx_w / 2, -by_w / 2, by_w / 2],
              vmax=0.8 * np.abs(feld_o).max())
    ax.axhline(-T_WALK / 2, color="w", lw=0.7)
    ax.axhline(T_WALK / 2, color="w", lw=0.7)
    ax.annotate("", xy=(1.6, 1.6), xytext=(-1.6, -1.6),
                arrowprops=dict(color="w", arrowstyle="<->", lw=1.0,
                                alpha=0.7))
    ax.set_title(f"{name}: Versatz {versatz:+.2f}")
    ax.set_xlabel("x")
achsen[0].set_ylabel("y")
fig.tight_layout()
plt.show()
Abbildung 28.7: Das Doppelbild im Experiment: dasselbe Bündel fällt senkrecht auf den Kalkspat mit 45°-Achse (Doppelpfeil). Links die ordentliche Welle (Ez): sie läuft exakt gerade durch (Versatz 0,000). Rechts die außerordentliche (Hz): ihre Wellenfronten bleiben waagerecht, aber die Energie wandert um 6,2° zur Seite — unten tritt der Strahl 0,65 versetzt aus. Unpolarisiertes Licht enthält beide Anteile: zwei Bilder.

Der ordentliche Strahl: \(0{,}000\). Der außerordentliche: \(-0{,}65\) — Betrag exakt auf der Theorie (das Vorzeichen sagt nur, zu welcher Seite die Achse kippt). Senkrechter Einfall, seitlicher Auszug: Die Wellenfronten laufen stur geradeaus, aber die Energie schert aus — Poynting und \(\vec{k}\) sind im anisotropen Kristall keine Parallelen mehr. Beim NIM zeigten sie in entgegengesetzte Richtungen, hier in verschiedene — zwei Kapitel-18-Pointen in einem Kapitel. Dreht man den Kristall auf der Buchseite, rotiert das zweite Bild ums erste: Die Walk-off-Richtung folgt der Achse.

28.9 Die λ/4-Platte: aus linear wird zirkular

Zum Schluss lösen wir das älteste offene Versprechen dieses Buchs ein. In Kapitel 4 stand: Überlagert man zwei senkrecht zueinander polarisierte Wellen mit 90° Phasenverzug, rotiert der Summenvektor — zirkulare Polarisation. Nur: Wie macht man die 90°? Jetzt haben wir das Werkzeug. In einem doppelbrechenden Plättchen laufen die beiden Feldkomponenten verschieden schnell — die längs der „schnellen” Achse (\(n_s\)) eilt der längs der „langsamen” (\(n_l\)) davon. Der Phasenverzug wächst mit der Dicke:

\[ \Delta\varphi = 2\pi f\, d\,(n_l - n_s) \stackrel{!}{=} \frac{\pi}{2} \quad\Longrightarrow\quad d_{\lambda/4} = \frac{\lambda}{4\,(n_l - n_s)} . \]

Mit \(n_s = 1{,}5\) und \(n_l = 1{,}6\) (ein künstlicher Kristall mit ordentlich Kontrast — echter Quarz hat nur \(\Delta n = 0{,}009\) und braucht entsprechend dickere oder trickreichere Platten) wird \(d = \lambda / 0{,}4 / 1 = 2{,}5\).

Versuchsaufbau. Frage: Macht das Plättchen aus 45°-linearem Licht zirkulares? Bühne: Eigentlich ist das ein 1D-Problem — aber Meeps dimensions=1 kennt nur eine einzige E-Komponente, und wir brauchen zwei (Werkstatt-Notiz). Der Ausweg: eine schmale 2D-Zelle (Breite \(0{,}25\)) mit periodischem Rand (k_point = 0), die Welle läuft in \(y\). Die beiden Polarisationen sind dann \(E_z\) und \(E_x\) — und der Kristalltensor \(\mathrm{diag}(n_s^2,\ \cdot\ , n_l^2)\) gibt jeder ihren Index (\(\varepsilon_y\) spielt bei senkrechtem Durchlauf keine Rolle). Anregung: beide Komponenten gleich stark (das ist die 45°-Linearpolarisation), Dauerstrich \(f = 1\), force_complex_fields — wir wollen die Phasen. Messgröße: die komplexen Amplituden \(E_x\) und \(E_z\) an einem Punkt hinter dem Plättchen: ihr Phasenunterschied und ihr Betragsverhältnis. Erfolgskriterium: \(90°\) (modulo Vorzeichen — welcher von beiden eilt, hängt davon ab, welche Achse die schnelle ist) und Verhältnis nahe 1.

N_SCHNELL, N_LANGSAM = 1.5, 1.6
D_VIERTEL = 1 / (4 * F_OP * (N_LANGSAM - N_SCHNELL))
MAT_PLATTE = mp.Medium(epsilon_diag=mp.Vector3(N_SCHNELL**2, 1.0,
                                               N_LANGSAM**2))


def wellenplatte(dicke, until=300):
    """Quasi-1D: schmale 2D-Zelle, k_point=0, Welle läuft in y."""
    dpml, pad = 2.0, 4.0
    sy = 2 * (dpml + pad) + dicke
    sx = 0.25
    quellen = [mp.Source(mp.ContinuousSource(F_OP), component=c,
                         center=mp.Vector3(0, -sy / 2 + dpml + 1.0),
                         size=mp.Vector3(sx, 0))
               for c in (mp.Ez, mp.Ex)]      # 45°: beide gleich stark
    sim = mp.Simulation(cell_size=mp.Vector3(sx, sy), resolution=64,
                        Courant=0.25, force_complex_fields=True,
                        k_point=mp.Vector3(),
                        boundary_layers=[mp.PML(dpml, direction=mp.Y)],
                        geometry=[mp.Block(size=mp.Vector3(mp.inf, dicke,
                                                           mp.inf),
                                           center=mp.Vector3(),
                                           material=MAT_PLATTE)],
                        sources=quellen)
    sim.run(until=until)
    messpunkt = mp.Vector3(0, sy / 2 - dpml - 1.0)
    vor_punkt = mp.Vector3(0, -dicke / 2 - 0.5)
    return ((complex(sim.get_field_point(mp.Ex, messpunkt)),
             complex(sim.get_field_point(mp.Ez, messpunkt))),
            (complex(sim.get_field_point(mp.Ex, vor_punkt)),
             complex(sim.get_field_point(mp.Ez, vor_punkt))))


(ex_nach, ez_nach), (ex_vor, ez_vor) = wellenplatte(D_VIERTEL)
verzug = np.degrees(np.angle(ex_nach) - np.angle(ez_nach)) % 360
print(f"d = {D_VIERTEL}:")
print(f"  Phasenverzug Ex gegen Ez: {verzug:.1f}°  (Soll 270° ≡ −90°)")
print(f"  Betragsverhältnis |Ex|/|Ez|: {abs(ex_nach) / abs(ez_nach):.3f}")
d = 2.499999999999998:
  Phasenverzug Ex gegen Ez: 269.5°  (Soll 270° ≡ −90°)
  Betragsverhältnis |Ex|/|Ez|: 0.924

Der Verzug trifft die Viertelwelle (es sind \(-90°\), weil bei uns die \(E_x\)-Achse die schnelle ist und ihre Phase weniger nachläuft). Das Betragsverhältnis ist \(0{,}92\) statt \(1{,}00\) — die beiden Achsen sehen leicht verschiedene Fresnel-Verluste und Mehrfachechos im Plättchen (Kapitel 27 lässt grüßen: Jedes Plättchen ist auch ein schwaches Etalon; reale Wellenplatten werden deshalb entspiegelt). Eine \(0{,}92\)-Ellipse ist von einem Kreis mit bloßem Auge kaum zu unterscheiden.

Was haben wir gebaut? Vor dem Plättchen schwingen \(E_x\) und \(E_z\) im Takt — der Summenvektor pendelt auf einer 45°-Linie. Dahinter ist \(E_x\) um eine Viertelperiode verschoben: Wenn \(E_z\) maximal ist, geht \(E_x\) gerade durch null und umgekehrt — der Summenvektor rotiert. Das ist zirkulare Polarisation, das Versprechen aus Kapitel 4, eingelöst mit einem Stück anisotropem Material.

Der Film zeichnet den E-Vektor als Uhrzeiger: links die Eingangswelle (45°-linear — per Konstruktion, beide Quellen gleich stark), rechts die gemessenen komplexen Amplituden hinter dem Plättchen, über eine volle Periode abgerollt. Sieh links das Pendeln auf der Linie — und rechts den Kreisverkehr. (Würden wir direkt vor dem Plättchen messen, sähe die Linie leicht verbogen aus — das schwache Etalon-Echo von eben mischt dort mit.)

Code der Animation (nur in der HTML-Fassung)
# von oben: (ex_nach, ez_nach) — gemessene komplexe Amplituden
from matplotlib import animation
from IPython.display import HTML

phasen_p = np.linspace(0, 2 * np.pi, 32, endpoint=False)
fein = np.linspace(0, 2 * np.pi, 200)
fig, (a1, a2) = plt.subplots(1, 2, figsize=(8.4, 4.4))
pfeile, spuren = [], []
for a, (ex, ez), name in [(a1, (1 + 0j, 1 + 0j),
                           "Eingangswelle (45°-linear)"),
                          (a2, (ex_nach, ez_nach),
                           "hinter dem Plättchen (gemessen)")]:
    norm = max(abs(ex), abs(ez))
    bahn_x = np.real(ex / norm * np.exp(-1j * fein))
    bahn_z = np.real(ez / norm * np.exp(-1j * fein))
    a.plot(bahn_x, bahn_z, "-", color="0.75", lw=1.0)   # die volle Bahn
    pfeil, = a.plot([0, bahn_x[0]], [0, bahn_z[0]], "C0-", lw=2.2)
    punkt, = a.plot([bahn_x[0]], [bahn_z[0]], "C0o", ms=7)
    pfeile.append((pfeil, punkt, ex / norm, ez / norm))
    a.set_xlim(-1.25, 1.25)
    a.set_ylim(-1.25, 1.25)
    a.set_aspect("equal")
    a.set_xlabel("$E_x$")
    a.set_title(name, fontsize=10)
a1.set_ylabel("$E_z$")
titel = fig.suptitle("", fontsize=10)
fig.subplots_adjust(top=0.82)


def bild(i):
    zeiger = np.exp(-1j * phasen_p[i])
    artists = []
    for pfeil, punkt, exn, ezn in pfeile:
        x, z = np.real(exn * zeiger), np.real(ezn * zeiger)
        pfeil.set_data([0, x], [0, z])
        punkt.set_data([x], [z])
        artists += [pfeil, punkt]
    titel.set_text("Die λ/4-Platte als Polarisations-Dreher: "
                   "Linie hinein, Kreis heraus\n"
                   f"t = {phasen_p[i] / (2 * np.pi):.2f} Perioden")
    return artists + [titel]


anim = animation.FuncAnimation(fig, bild, frames=len(phasen_p), interval=90)
plt.close(fig)
HTML(anim.to_jshtml(default_mode="loop"))

Die Spielwiese rechnet beliebige Platten mit dem Jones-Kalkül (die Zwei-Zeilen-Buchhaltung der Polarisationsoptik: ein Vektor \((E_x, E_z)\), eine Matrix pro Bauteil). Drei Aufträge, erst vorhersagen: (1) GAMMA_GRAD = 180 — was wird aus der 45°-Linie? (2) Verzug 90°, aber ACHSE_GRAD = 0 (Licht schwingt längs der schnellen Achse) — warum passiert dann gar nichts? (3) Bei welchem Verzug entsteht die „fetteste” Ellipse zwischen Linie und Kreis?

Wozu das alles? Zirkular polarisiertes Licht steckt in 3D-Kinobrillen (links- und rechtsdrehend für die beiden Augen — und es bleibt zirkular, wenn man den Kopf neigt, anders als bei linearen Filtern), in CD/DVD-Laufwerken (der Isolator aus Polarisator + λ/4-Platte schützt die Laserdiode vor ihrem eigenen Reflex) und in jedem LCD-Bildschirm — dessen Flüssigkristalle sind nichts anderes als schaltbare Wellenplatten, deren Verzug eine Spannung einstellt (mehr im Kleingedruckten).

28.10 Das Kapitel-Programm

programme/kap28/kap28_metamaterial_doppelbrechung.py bündelt alle Befunde eigenständig und mit assert-Schranken: die Eichung \(\varepsilon(1) = \mu(1) = -1\), die Smith-Dreierreihe samt Transfer-Matrix-Vergleich und Lambert-Beer-Gegenprobe, das Phasenprofil (\(\pm 1{,}00\,k_0\)), die Gruppenlaufzeit gegen das handgerechnete \(n_g\), den Vorzeichenwechsel des Brechungsversatzes, die Flachlinse (Fokuslage und FWHM), die Doppelbrechungs-Winkel und -Indizes per 2D-FFT, den Walk-off (ordentlich exakt null), beide Wellenplatten und die Drahtgitter-Klippen. Läuft in wenigen Minuten.

TippMerkkasten
  • Metamaterial ≠ Kristall: Bausteine \(a \ll \lambda\) wirken als homogenes Medium (effektives \(\varepsilon\), \(\mu\)); Bausteine \(a \sim \lambda\) machen Interferenz (Kapitel 26).
  • Ein negatives Vorzeichen sperrt, zwei öffnen: \(\varepsilon < 0\) allein → evaneszent (Kapitel 16/17); \(\varepsilon, \mu < 0\) → Welle läuft mit \(n = -\sqrt{\varepsilon\mu}\). Energie vorwärts, Phase rückwärts (\(\vec{S}\) antiparallel \(\vec{k}\)); Brechung knickt auf die falsche Seite, eine flache Platte wird zur Linse (Bild bei \(t - d_1\)).
  • Negative Konstanten gibt es nur mit Dispersion — und \(\mu < 0\) nur nahe einer Resonanz (Split-Ring = Kapitel-19-Schwingkreis). Darum sind Verluste der Systempreis der Metamaterialien (Kapitel 15: die Absorption wohnt an der Resonanz).
  • Anisotropie = ein Index pro Feldrichtung: \(\varepsilon\) als Tensor (epsilon_diag/offdiag). Doppelbrechung sortiert nach Polarisation (\(n_o\)/\(n_e\)); bei schräger Achse scheren Energiefluss und Wellenvektor auseinander (Walk-off → Doppelbild).
  • Wellenplatten sind Phasen-Werkzeuge: \(d_{\lambda/4} = \lambda/(4\,\Delta n)\) macht aus 45°-linear zirkular (Kapitel-4-Zusage); die λ/2-Platte spiegelt die Polarisationsrichtung.
  • Werkstatt: µ-dispersive Materialien dürfen weder PML noch Absorber berühren (Rückwärtswellen drehen die Dämpfung um) — Blöcke endlich bauen, Luftpuffer lassen.

Roter Faden

Dieses Kapitel erntet quer durchs Buch: Die Drude-Klippe aus Kapitel 16 wird zur Stellschraube (Drahtgitter), der Schwingkreis aus Kapitel 19 wird zum Materialparameter (Split-Ring macht \(\mu(f)\) zur Lorentz-Linie aus Kapitel 15 — mitsamt ihrer Verluste), Snellius aus Kapitel 12 bekommt ein Vorzeichen, der Poynting-Vektor aus Kapitel 18 trennt sich vom Wellenvektor (antiparallel im NIM, schräg im Kristall), das Tunneln aus Kapitel 17 erklärt die weichen Flanken des Smith-Fensters, und die Polarisation aus Kapitel 4/13 wird endlich zum Werkzeug: Die λ/4-Platte löst die älteste offene Zusage des Buchs ein. Nach vorn: Zweimal haben wir uns in diesem Kapitel auf „das gibt es nur mit Dispersion” berufen — Kapitel 29 zeigt, dass dahinter ein Satz steckt: Kausalität erzwingt Dispersion (Kramers-Kronig), und Real- und Imaginärteil von \(\varepsilon(f)\) sind keine unabhängigen Größen.

Übungen

Ü 28.1 — Die Flachlinse vermessen (Verstehen). Eine Veselago-Platte (\(n = -1\)) ist \(t = 4\) cm dick. (a) Eine Punktquelle sitzt \(d_1 = 1{,}5\) cm vor der Platte — wo liegen das innere und das äußere Bild? Konstruiere beide mit dem Snellius-Knick \(\theta_2 = -\theta_1\) (eine Skizze mit zwei, drei Strahlen genügt). (b) Was ändert sich bei \(d_1 = 3\) cm? (c) Was passiert für \(d_1 = 5\) cm — und was sagt das über den Unterschied zwischen dieser „Linse” und einer Glaslinse mit Brennweite?

Bei \(n = -1\) wird jeder Strahl an beiden Flächen exakt gespiegelt gebrochen (\(\theta_2 = -\theta_1\)): Ein Strahl, der unter dem Winkel \(\theta\) von der Quelle wegläuft, läuft in der Platte unter \(-\theta\) weiter und hinter ihr wieder unter \(\theta\). Die Geometrie ist dadurch ein reines Dreisatz-Problem: Alle Strahlen kreuzen die Achse erst bei der Tiefe \(d_1\) in der Platte (inneres Bild), dann im Abstand \(t - d_1\) hinter ihr.

import numpy as np
import matplotlib.pyplot as plt

T_P, D1_U = 4.0, 1.5
fig, ax = plt.subplots(figsize=(7.0, 3.6))
ax.axvspan(0, T_P, color="0.88")
for steigung in [0.5, 0.25, -0.25, -0.5]:
    # Quelle bei x = −d1; Knick an x = 0, Gegenknick in der Platte
    x_pfad = [-D1_U, 0, T_P, T_P + (T_P - D1_U)]
    y_pfad = [0, steigung * D1_U,
              steigung * D1_U - steigung * T_P,
              0]
    ax.plot(x_pfad, y_pfad, "C0-", lw=1.2)
ax.plot(-D1_U, 0, "ko", ms=6)
ax.plot(D1_U, 0, "C3x", ms=8)
ax.plot(T_P + (T_P - D1_U), 0, "C3x", ms=8)
ax.annotate("Quelle", (-D1_U, 0.15), ha="center")
ax.annotate("inneres Bild", (D1_U, -0.45), ha="center", color="C3")
ax.annotate("äußeres Bild", (T_P + (T_P - D1_U), 0.3), ha="center",
            color="C3")
ax.set_xlabel("Achse (grau: die Platte, t = 4)")
ax.set_ylabel("Querrichtung")
ax.set_ylim(-1.4, 1.4)
ax.grid(alpha=0.3)
plt.show()
print(f"(a) d1 = 1,5: inneres Bild {D1_U} cm tief, "
      f"äußeres {T_P - D1_U} cm hinter der Platte")
print(f"(b) d1 = 3,0: inneres Bild 3 cm tief, äußeres {T_P - 3.0} cm dahinter")

Strahlkonstruktion zu Ü 28.1 (a): Jeder Strahl knickt an beiden Flächen auf die andere Seite der Normalen. Inneres Bild bei Tiefe d₁ = 1,5 in der Platte, äußeres bei t − d₁ = 2,5 dahinter.
(a) d1 = 1,5: inneres Bild 1.5 cm tief, äußeres 2.5 cm hinter der Platte
(b) d1 = 3,0: inneres Bild 3 cm tief, äußeres 1.0 cm dahinter

(c) Für \(d_1 = 5 > t\) gibt es kein Bild mehr: Die Strahlen knicken zwar, kreuzen die Achse aber schon nicht mehr innerhalb der Platte (\(d_1\) passt nicht hinein) — und hinter der Platte laufen sie auseinander wie zuvor. Die Veselago-Platte hat also keine Brennweite: Sie kann nur abbilden, was näher als \(t\) vor ihr liegt, und sie kann parallele Strahlen (Quelle im Unendlichen) überhaupt nicht bündeln. Eine Glaslinse bildet jede Objektweite ab — die Flachlinse verschiebt jedes Nahfeld-Bild nur um die feste Strecke \(2t\) nach hinten (\(d_1 + \text{Bildabstand} = d_1 + t - d_1 + t\)… besser gemerkt: Quelle→Bild ist immer \(2t\)).

Ü 28.2 — Die λ/2-Platte (Verändern). Verdopple im Wellenplatten-Aufbau die Dicke auf \(d = 5\). (a) Sage den Phasenverzug vorher. (b) Miss ihn nach. (c) Die Eingangswelle war 45°-linear (\(E_x = E_z\)). Was für eine Polarisation kommt heraus? Zeichne die E-Vektor-Bahn aus den gemessenen Amplituden und bestimme ihren Winkel — welche Operation führt die λ/2-Platte an der Polarisationsrichtung aus?

(a) Doppelte Dicke, doppelter Verzug: \(180°\).

(ex_h, ez_h), _ = wellenplatte(2 * D_VIERTEL)
verzug_h = np.degrees(np.angle(ex_h) - np.angle(ez_h)) % 360
print(f"(b) Verzug = {verzug_h:.1f}°  (Soll 180°), "
      f"|Ex|/|Ez| = {abs(ex_h) / abs(ez_h):.4f}")

# (c) E-Vektor-Bahn und ihr Hauptachsen-Winkel
fein = np.linspace(0, 2 * np.pi, 400)
bahn_x = np.real(ex_h * np.exp(-1j * fein))
bahn_z = np.real(ez_h * np.exp(-1j * fein))
i_spitze = np.argmax(np.hypot(bahn_x, bahn_z))
winkel_bahn = np.degrees(np.arctan2(bahn_x[i_spitze], bahn_z[i_spitze]))
print(f"(c) Bahn-Winkel = {winkel_bahn:+.1f}° gegen die z-Achse "
      f"(Eingang: +45°)")
(b) Verzug = 179.4°  (Soll 180°), |Ex|/|Ez| = 1.0000
(c) Bahn-Winkel = -45.0° gegen die z-Achse (Eingang: +45°)

(c) Heraus kommt wieder lineares Licht (das Betragsverhältnis ist 1,000 und die Bahn eine Linie), aber die Schwingungsrichtung liegt bei \(-45°\): Die λ/2-Platte spiegelt die Polarisationsrichtung an ihren Achsen — aus \(+45°\) wird \(-45°\), allgemein aus dem Winkel \(\alpha\) der Winkel \(-\alpha\). Deshalb ist sie das Standard-Werkzeug, um Polarisationen zu drehen (Platte um \(\beta/2\) drehen → Polarisation dreht um \(\beta\)), etwa vor polarisierenden Strahlteilern.

Ü 28.3 — Das Drahtgitter in echt (Übertragen). Im Kapitel haben wir den Drähten einfach eine Drude-Formel zugeschrieben. Baue das Metamaterial-Atom jetzt wirklich: eine einzelne Reihe-um-Reihe-Zelle (\(s_y = a = 0{,}2\), periodischer Rand mit k_point = 0) mit zehn PEC-Stäbchen (mp.metal) im Abstand \(a = 0{,}2\) quer zum Strahl, Radius \(r = 0{,}01\) — das Gitter ist also viel feiner als die Wellenlänge (\(a/\lambda = 0{,}2\) bei \(f = 1\)). (a) Miss \(T(f)\) von \(0{,}1\) bis \(2{,}5\) und finde die Klippe. (b) Wiederhole mit \(r = 0{,}02\): Wandert die Klippe nach oben oder unten? Begründe mit dem Induktivitäts-Argument aus Kapitel 22. (c) Die Pendry-Faustformel \(f_p \approx \sqrt{2\pi/\ln(a/r)}\,/\,(2\pi a)\) liefert \(1{,}15\) bzw. \(1{,}3\) — du wirst deutlich höhere Klippen messen. Sieh dir \(a/\lambda\) an der Klippe an: In welchem Regime sind wir dort, und welches Kapitel ist dann zuständig?

A_G, N_REIHEN = 0.2, 10


def drahtgitter_sweep(radius, fcen=1.3, df=2.4, nfreq=361):
    dpml, pad = 1.0, 2.0
    breite = N_REIHEN * A_G
    sx = 2 * (dpml + pad) + breite
    sy = A_G
    leistungen = []
    for mit_gitter in (False, True):
        geometrie = [mp.Cylinder(radius=radius,
                                 center=mp.Vector3(-breite / 2 + A_G / 2
                                                   + i * A_G, 0),
                                 material=mp.metal)
                     for i in range(N_REIHEN)] if mit_gitter else []
        quelle = [mp.Source(mp.GaussianSource(fcen, fwidth=df),
                            component=mp.Ez,
                            center=mp.Vector3(-sx / 2 + dpml + 0.5, 0),
                            size=mp.Vector3(0, sy))]
        sim = mp.Simulation(cell_size=mp.Vector3(sx, sy), resolution=200,
                            k_point=mp.Vector3(),
                            boundary_layers=[mp.PML(dpml, direction=mp.X)],
                            geometry=geometrie, sources=quelle)
        x_mess = sx / 2 - dpml - 0.5
        tr = sim.add_flux(fcen, df, nfreq,
                          mp.FluxRegion(center=mp.Vector3(x_mess, 0),
                                        size=mp.Vector3(0, sy)))
        sim.run(until_after_sources=mp.stop_when_fields_decayed(
            50, mp.Ez, mp.Vector3(x_mess, 0), 1e-9))
        leistungen.append(np.array(mp.get_fluxes(tr)))
    fs = np.array(mp.get_flux_freqs(tr))
    return fs, leistungen[1] / leistungen[0]


fig, ax = plt.subplots(figsize=(7.0, 3.8))
for radius, farbe in [(0.01, "C0"), (0.02, "C3")]:
    fs_d, T_d = drahtgitter_sweep(radius)
    ueber = T_d > 0.5
    klippe = fs_d[np.argmax(ueber)] if ueber.any() else float("nan")
    pendry = np.sqrt(2 * np.pi / np.log(A_G / radius)) / (2 * np.pi * A_G)
    ax.plot(fs_d, T_d, "-", color=farbe,
            label=f"r = {radius}: Klippe {klippe:.2f} (Pendry {pendry:.2f})")
    ax.axvline(klippe, color=farbe, lw=0.8, ls=":")
    print(f"r = {radius}: Klippe (T = 0,5) bei f = {klippe:.3f}, "
          f"Pendry-Formel {pendry:.2f}, a/lambda an der Klippe: "
          f"{A_G * klippe:.2f}")
ax.set_xlabel("Frequenz f")
ax.set_ylabel("Transmission T")
ax.set_title("Zehn Reihen dünner Stäbchen = ein Plasma")
ax.legend(fontsize=9)
ax.grid(alpha=0.3)
plt.show()
r = 0.01: Klippe (T = 0,5) bei f = 1.447, Pendry-Formel 1.15, a/lambda an der Klippe: 0.29
r = 0.02: Klippe (T = 0,5) bei f = 1.813, Pendry-Formel 1.31, a/lambda an der Klippe: 0.36

Das Drahtgitter als Plasma: Transmission durch zehn Stäbchen-Reihen für zwei Drahtradien. Unterhalb der Klippe ist das Gitter dicht wie die Ionosphäre für Kurzwelle — obwohl es zu über 99 % aus Luft besteht. Dickere Drähte verschieben die Klippe nach oben.

(a), (b) Die Klippe liegt für \(r = 0{,}01\) bei \(f \approx 1{,}45\) und wandert für \(r = 0{,}02\) auf \(\approx 1{,}8\)dickere Drähte sperren bis höher hinauf. Das Induktivitäts-Argument: \(f_p\) ist die Frequenz, bei der die träge Strom-Antwort des Gitters die Feld-Antwort gerade noch einholt. Die Trägheit kommt hier nicht von der Elektronenmasse, sondern von der Induktivität der dünnen Drähte (Kapitel 22: \(L' \propto \ln(a/r)\)) — dünnere Drähte sind induktiver, also träger, also sperrt ihr Gitter nur bis zu einer tieferen Frequenz.

(c) An der Klippe ist \(a/\lambda = 0{,}29\) bzw. \(0{,}36\) — das ist kein sauberes Metamaterial-Regime mehr (\(a \ll \lambda\) verlangt eher \(a/\lambda \lesssim 0{,}1\)), sondern der Vorhof von Kapitel 26: Bragg-Interferenz zwischen den Reihen mischt mit, und die hübsche Pendry-Formel (hergeleitet für \(a/\lambda \to 0\)) wird zur Größenordnungs-Schätzung. Merke: Ob etwas „Medium” oder „Kristall” ist, entscheidet nicht das Bauteil, sondern das Verhältnis \(a/\lambda\) — ein und dasselbe Drahtgitter ist bei tiefen Frequenzen ein Plasma und bei hohen ein photonischer Kristall.

Ü 28.4 — Die Achillesferse (Übertragen). Die Flachlinse aus dem Kapitel lief mit \(\gamma = 0{,}01\). Reale Metamaterialien — besonders optische — sind verlustreicher. Wiederhole den Linsenlauf mit \(\gamma = 0{,}05\) und \(\gamma = 0{,}15\) (in beiden Bausteinen). (a) Sage qualitativ vorher, was mit dem Fokus passiert. (b) Miss Fokuslage und Spitzenintensität. (c) Erkläre mit Kapitel 15, warum man dieses Problem nicht einfach „wegoptimieren” kann.

def nim_mit(gamma):
    return mp.Medium(
        epsilon=1,
        E_susceptibilities=[mp.DrudeSusceptibility(frequency=F_OP,
                                                   gamma=gamma, sigma=2.0)],
        mu=1,
        H_susceptibilities=[mp.LorentzianSusceptibility(frequency=F0_RING,
                                                        gamma=gamma,
                                                        sigma=S_RING)])


print("gamma   Fokuslage   Spitzenintensität")
for gamma_u in [0.01, 0.05, 0.15]:
    ez_u, bx_u, by_u = flachlinse(material=nim_mit(gamma_u))
    nx_u, ny_u = ez_u.shape
    xs_u = np.linspace(-bx_u / 2, bx_u / 2, nx_u)
    ys_u = np.linspace(-by_u / 2, by_u / 2, ny_u)
    int_u = np.abs(ez_u)**2
    profil_u = np.where(ys_u > 1.2, int_u[np.argmin(np.abs(xs_u)), :], 0)
    iy_u = np.argmax(profil_u)
    print(f"{gamma_u:5.2f}   y = {ys_u[iy_u]:+.2f}    {profil_u[iy_u]:.4f}")
gamma   Fokuslage   Spitzenintensität
 0.01   y = +1.91    0.4449
 0.05   y = +1.76    0.0041
 0.15   y = +1.29    0.0000

(a), (b) Der Fokus stirbt dramatisch: Schon \(\gamma = 0{,}05\) kostet einen Faktor 100 an Intensität (von \(0{,}44\) auf \(0{,}004\)), bei \(0{,}15\) bleibt praktisch nichts — und was bleibt, rückt zur Platte hin (das „Bild” ist dann nur noch das durchsickernde Nahfeld der Plattenrückseite). Zum Vergleich: \(\gamma = 0{,}05\) heißt im Lorentz-Modell immer noch ein Resonator-\(Q\) von \(f_0/\gamma = 16\) — kein schlechter Schwingkreis.

(c) Das Kapitel-15-Argument: \(\mu < 0\) existiert nur nahe der Ringresonanz, und ein getriebener Resonator setzt genau dort am meisten Leistung um — Verlust und Funktion wohnen an derselben Adresse, man kann nicht das eine behalten und das andere räumen. Bei Mikrowellen (Kupferringe, hohes \(Q\)) geht es leidlich; im Sichtbaren müssen die „Ringe” aus Plasmonik-Metallen gebaut werden (Kapitel 16: Gold mit seinen Drude-Verlusten), und dort hat genau diese Rechnung der Pendry-Superlinse die Flügel gestutzt. Merksatz fürs Berufsleben: Wenn ein Wundermaterial eine Resonanz braucht, frage zuerst nach dem Imaginärteil.

Das Kleingedruckte

Die perfekte Linse und ihr Ende. Pendrys Rechnung von 2000 zeigte: Die verlustfreie \(n = -1\)-Platte verstärkt evaneszente Wellen exponentiell wieder auf — die Information über Sub-Wellenlängen-Details, die nach Kapitel 12 im Abstand \(\sim\lambda\) verloren geht, käme im Bild vollständig an („perfect lens”, Auflösung unbegrenzt). Der Haken steckt im selben Exponenten: Schon winzige Verluste kappen die Verstärkung der feinsten Details zuerst, und mit realistischen \(\mathrm{Im}(\varepsilon)\) bleibt von „perfekt” ein Faktor zwei, drei unter der Beugungsgrenze — im Nahfeld, als „Superlinse” aus einem schlichten Silberfilm sogar experimentell bestätigt (2005). Unser FWHM von \(0{,}41\lambda\) statt \(0{,}5\lambda\) ist ein matter Abglanz dieses Effekts.

Transformationsoptik und Tarnkappen. Die Metamaterial-Idee in voller Größe: Weil Maxwell unter Koordinatentransformationen forminvariant ist, kann man eine gedachte Verbiegung des Raums (etwa: „alle Strahlen fließen um diese Kugel herum”) in ein Rezept für ortsabhängige \(\varepsilon\)- und \(\mu\)-Tensoren übersetzen. Die Mikrowellen-„Tarnkappe” von 2006 (Schurig/Smith) war genau das — ein Ring aus Split-Ring-Zellen mit radial verlaufenden Materialprofilen. Die Grenzen sind dieselben wie immer in diesem Kapitel: schmalbandig (Dispersion!) und verlustig.

Flüssigkristalle: die schaltbare Wellenplatte. Die stäbchenförmigen Moleküle eines nematischen Flüssigkristalls sind anisotrop polarisierbar (\(\Delta n\) bis \(0{,}3\) — mehr als Kalkspat) und lassen sich mit ein paar Volt ausrichten. Ein LCD-Pixel ist eine Zelle, deren Verzug die Spannung zwischen „λ/2-Platte” (Licht passiert den gekreuzten Ausgangspolarisator) und „neutral” (Licht wird geschluckt) durchstimmt. Dass dein Monitor unter polarisierter Sonnenbrille seltsam aussieht, ist angewandte Doppelbrechung.

Chiralität und Bianisotropie. Unsere Tensoren waren symmetrisch und reell. Materialien mit Schraubenstruktur (Zuckerlösung, manche Metamaterial-Zellen) koppeln \(\vec{E}\) an \(\vec{B}\) — sie drehen die Polarisationsebene kontinuierlich (optische Aktivität) oder unterscheiden links- und rechtszirkular in der Absorption. Der Saccharimeter des Zuckersieders ist Polarisationsoptik von 1840, die bis heute in jeder Zuckerfabrik steht.

Negativer Index ohne Ringe. Photonische Kristalle (Kapitel 26) können nahe der Bandkante effektiv \(n < 0\) zeigen — über Interferenz statt über Resonanzen, mit geringeren Verlusten, aber nur in schmalen Bändern und Winkelbereichen. Auch hier gilt: Die Natur gibt das Vorzeichen nicht her, man muss es ihr abringen.