Was ist FDTD?

1. Einleitung

Die Methode der endlichen Differenzen im Zeitbereich (FDTD) löst die Maxwell-Gleichungen – die Licht beschreiben – direkt im Zeitbereich auf einem diskreten räumlichen Gitter. Im Gegensatz zu Raytracing, Modenlösungsverfahren oder Transfermatrix-Ansätzen trifft FDTD keine vereinfachenden Annahmen über die Ausbreitung oder die Modenstruktur. Stattdessen werden die elektrischen und magnetischen Felder explizit im Laufe der Zeit durch Diskretisierung der Rotationsgleichungen – des Ampèreschen Gesetzes und des Faradayschen Gesetzes – weiterentwickelt, sodass sich elektromagnetische Wellen auf natürliche Weise im gesamten Simulationsbereich ausbreiten, streuen, interferieren und beugen. Diese Allgemeingültigkeit macht FDTD sehr genau und vielseitig, stellt aber auch hohe Anforderungen an den Speicher und die Rechenzeit.

\[ \nabla \times \boldsymbol{H} = \sigma \boldsymbol{E} + \epsilon \frac{\partial \boldsymbol{E}}{\partial t}, \qquad \nabla \times \boldsymbol{E} = - \mu \frac{\partial \boldsymbol{H}}{\partial t}. \]

Da FDTD im Wesentlichen keine Annahmen über die Ausbreitungsrichtung, die Struktur des Bauelements oder das Feldverhalten trifft, kann eine einzige Simulation ein breites Frequenzspektrum, Beugungs- und Interferenzphänomene, Nahfeld-Effekte, Strukturen mit Abmessungen unterhalb der Wellenlänge und transiente elektromagnetische Phänomene innerhalb eines einzigen Rahmens erfassen. Diese Allgemeingültigkeit macht FDTD so leistungsstark, hat aber auch ihren Preis: Das räumliche Gitter muss die kleinste physikalische Längenskala auflösen, die Zeitschrittweite ist durch die numerische Stabilität stark begrenzt, und dreidimensionale Simulationen werden schnell sowohl in Bezug auf den Speicherbedarf als auch auf die Rechenzeit anspruchsvoll, selbst für Geometrien, die auf den ersten Blick nicht sehr komplex erscheinen.

2. Warum Sie die FDTD-Methode nicht verwenden sollten ...

Bei den meisten praktischen optischen Problemen sollte die FDTD-Methode nicht die erste Wahl sein. Sie ist rechenintensiv, und für viele Systeme gibt es deutlich effizientere Ansätze. Die Entwicklung von Linsensystemen wie Kameras oder Teleskopen, die Modellierung großflächiger Bildgebungssysteme oder die Verfolgung von Lichtstrahlen durch makroskopische optische Anordnungen lassen sich viel besser mit Hilfe der Strahlenverfolgungsmethode bewerkstelligen. Ebenso lassen sich planare Mehrschichtstrukturen wie Dünnschicht-Solarzellen, OLED-Stapel und optische Beschichtungen viel effizienter mit Hilfe von Transfermatrixmethoden behandeln. In diesen Fällen verursacht die FDTD-Methode einen erheblichen zusätzlichen Rechenaufwand, ohne zusätzliche physikalische Erkenntnisse zu liefern.

Die meisten Methoden zur elektromagnetischen Simulation basieren letztendlich auf den Maxwell-Gleichungen. Verfahren wie Modenlösungsverfahren und Transfermatrixmethoden werden durch die Einführung vereinfachender Annahmen über die Geometrie, die Dimensionalität oder die Feldstruktur entwickelt, um die Berechnungen zu beschleunigen. Diese Näherungen sind äußerst wirksam, aber nur, wenn die zugrunde liegenden Annahmen zutreffen.

Die FDTD-Methode steht am anderen Ende dieses Spektrums. Sie ist äußerst leistungsfähig, insbesondere für komplexe dreidimensionale Systeme, unregelmäßige Geometrien und Probleme, bei denen das Feld nicht in einfache Modi oder Schichten zerlegt werden kann – wie beispielsweise bei photonischen Kristallen, Strukturen mit Abmessungen unterhalb der Wellenlänge, stark streuenden Medien und Nahfeldkopplungen. Der Preis für diese Flexibilität ist der hohe Rechenaufwand, sowohl in Bezug auf den Speicherbedarf als auch auf die Rechenzeit.

Entscheidend ist die Wahl des richtigen Werkzeugs. Verwenden Sie vereinfachte Modelle, wenn dies physikalisch möglich ist, und setzen Sie FDTD ein, wenn diese Annahmen nicht mehr gelten. Die Stärke von FDTD liegt in seiner Flexibilität und physikalischen Vollständigkeit; der Nachteil ist der höhere Rechenaufwand.

3. Die Maxwell-Gleichungen als Grundlage für die FDTD-Methode

Im Kern ist die FDTD-Methode sehr einfach. Man nimmt die Maxwell-Gleichungen, zerlegt sie in ihre Komponenten für jede Raumrichtung und löst diese Gleichungen dann direkt auf einem Gitter. Es gibt keine komplizierten Tricks – die Methode ist konzeptionell unkompliziert und in vielerlei Hinsicht einfacher als komplexere numerische Ansätze wie Drift-Diffusions-Solver. Zum Einstieg erweitern wir das Ampèresche Gesetz. In einem linearen Medium mit Leitfähigkeit \(\sigma\), Permittivität \(\epsilon\) und Permeabilität \(\mu\) kann dies wie folgt geschrieben werden:

\[ \sigma \boldsymbol{E} + \epsilon \frac{\partial \boldsymbol{E}}{\partial t} = \nabla \times \boldsymbol{H} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ H_x & H_y & H_z \end{vmatrix} \]

was wie folgt erweitert werden kann:

Das Faradaysche Gesetz lautet wie folgt:

\[ -\sigma_m \boldsymbol{H} - \mu \frac{\partial \boldsymbol{H}}{\partial t} = \nabla \times \boldsymbol{E} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ E_x & E_y & E_z \end{vmatrix} \]

was wie folgt erweitert werden kann:

Die in Blau und Rot hervorgehobenen Gleichungen stellen die Definition von FDTD dar und bilden das Kernstück der Methode. Es gibt nichts Komplizierteres. Alles andere ist eine Frage der Umsetzung: Ersetzen von Ableitungen durch finite Differenzen, Auswahl von Randbedingungen, Definition von Lichtquellen (Anregungen), Auswertung der Ergebnisse mithilfe von Detektoren und Schreiben von effizientem (oft parallelem) Code, um die Entwicklung der Felder im Zeitverlauf zu simulieren. Diese zusätzlichen Komponenten machen einen praktischen FDTD-Solver komplexer, ändern aber nicht die grundlegende Einfachheit der Kernmethode.

4. Ein einfaches Beispiel für eine eindimensionale FDTD-Simulation mit einem Spielzeug

Um zu zeigen, wie einfach die FDTD-Methode ist – sowohl konzeptionell als auch in der Umsetzung –, reduzieren wir die vollständigen Gleichungen nun auf ein eindimensionales Modell, das in wenigen Dutzend Zeilen MATLAB- oder Python-Code implementiert werden kann. Ziel ist es hier nicht, die physikalischen Grundlagen zu verändern, sondern unnötige Komplexität zu eliminieren, damit die Struktur der Methode vollständig nachvollziehbar ist.

Um dies zu erreichen, betrachten wir die Variation nur entlang der \(x\)-Achse und setzen die übrigen räumlichen Ableitungen explizit auf Null: \[ \frac{\partial}{\partial y} = 0, \qquad \frac{\partial}{\partial z} = 0. \]. Dadurch entsteht ein transversales elektromagnetisches System mit einem elektrischen Feld \(E_z(x,t)\) und einem magnetischen Feld \(H_y(x,t)\).

Unter diesen Annahmen lassen sich Maxwells Wirbelgleichungen wie folgt vereinfachen:

\[ \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x} \]

\[ \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x}. \]

Dies sind die grundlegenden kontinuierlichen Gleichungen für den eindimensionalen Fall. In dieser Form können sie jedoch nicht direkt auf einem Computer gelöst werden. Um dies zu ermöglichen, ersetzen wir die Ableitungen durch finite Differenzen auf einem diskreten Gitter. Dies ist das Wesen der Methode: eine allgemeine Übersicht finden Sie unter finite Differenzen.

Durch die Einführung eines räumlichen Gitters \(x_i = i\Delta x\) und von Zeitschritten \(t^n = n\Delta t\) sowie durch die Verwendung der üblichen, versetzten (alternierenden) Anordnung ergeben sich die folgenden Aktualisierungsgleichungen:

\[ E_z^{n+1}(i) = E_z^n(i) + \frac{\Delta t}{\epsilon(i)\,\Delta x} \left[ H_y^{n+\frac{1}{2}}\!\left(i+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i-\frac{1}{2}\right) \right] \]

\[ H_y^{n+\frac{1}{2}}\!\left(i+\frac{1}{2}\right) = H_y^{n-\frac{1}{2}}\!\left(i+\frac{1}{2}\right) + \frac{\Delta t}{\mu\,\Delta x} \left[ E_z^n(i+1) - E_z^n(i) \right]. \]

Dies ist bereits die vollständige FDTD-Methode in ihrer einfachsten Form. Das elektrische Feld wird anhand räumlicher Differenzen im magnetischen Feld aktualisiert, und das magnetische Feld wird anhand räumlicher Differenzen im elektrischen Feld aktualisiert, wobei sich die beiden Felder zeitlich abwechselnd aktualisieren.

Für ein nicht-magnetisches Medium \(\mu = \mu_0\) besteht eine Beziehung zwischen der Permittivität und dem Brechungsindex, die durch folgende Formel ausgedrückt wird:

\[ \epsilon(x) = \epsilon_0 n(x)^2. \]

Das bedeutet, dass eine Änderung des Brechungsindex einfach durch eine Veränderung von \(\epsilon\) über das Gitter hinweg erreicht wird. Im folgenden Beispiel betrachten wir ein Medium, das von \(n=1\) zu \(n=3\) übergeht, sodass ein von links ausgesandter Impuls an der Grenzfläche teilweise reflektiert und teilweise durchgelassen wird.

In der Praxis werden die oben genannten Gleichungen einfach in eine Zeitschleife eingebettet: Das Magnetfeld wird anhand des elektrischen Feldes aktualisiert, die Randbedingungen werden angewendet, dann wird das elektrische Feld anhand des Magnetfeldes aktualisiert, und der Prozess wiederholt sich. Die Felder werden zeitlich versetzt, sodass sie in einem einfachen, expliziten Schema schrittweise weiterentwickelt werden, wobei jede Aktualisierung nur von lokalen, benachbarten Werten abhängt. Dies wird schematisch in ?? dargestellt und stellt im Wesentlichen alles dar, was der Computer bei jedem Zeitschritt ausführt.

5. Ein einfaches Simulationsbeispiel für ein Spielzeug mithilfe der 1D-FDTD-Methode in Matlab/Python.

Das gleiche eindimensionale FDTD-Verfahren lässt sich mit nur wenigen Dutzend Codezeilen implementieren. Im Folgenden sind äquivalente Implementierungen in MATLAB und Python aufgeführt. Beide Programme erzeugen einen Gauß-Impuls und zeigen die Reflexion an einer Sprungstelle des Brechungsindex.

clear;
clc;

% Physical constants
c0   = 299792458;
mu0  = 4*pi*1e-7;
eps0 = 1/(mu0*c0^2);

% Grid
Nx = 400;
dx = 20e-9;
dt = dx/(2*c0);

% Material profile
n = ones(1, Nx);
n(220:end) = 3.0;
eps_r = n.^2;

% Fields
Ez = zeros(1, Nx);
Hy = zeros(1, Nx-1);

% Source
source_pos = 80;
t0 = 40;
spread = 12;

Nt = 600;

figure;

for t = 1:Nt

    % Update H
    for i = 1:Nx-1
        Hy(i) = Hy(i) + (dt/(mu0*dx))*(Ez(i+1) - Ez(i));
    end

    % Update E
    for i = 2:Nx-1
        Ez(i) = Ez(i) + (dt/(eps0*eps_r(i)*dx))*(Hy(i) - Hy(i-1));
    end

    % Source
    Ez(source_pos) = Ez(source_pos) + exp(-((t - t0)/spread)^2);

    % Simple absorbing boundaries
    Ez(1)   = Ez(2);
    Ez(end) = Ez(end-1);

    % Plot every few steps
    if mod(t,5) == 0
        plot(1:Nx, Ez, 'b', 'LineWidth', 1.5);
        hold on;
        plot([220 220], [-1.5 1.5], '--r', 'LineWidth', 1.0);
        hold off;
        xlabel('Grid index');
        ylabel('E_z');
        title(['1D FDTD pulse reflecting from an index step, t = ', num2str(t)]);
        axis([1 Nx -1.5 1.5]);
        grid on;
        drawnow;
    end
end

import numpy as np
import matplotlib.pyplot as plt

c0 = 299792458.0
mu0 = 4.0e-7 * np.pi
eps0 = 1.0 / (mu0 * c0**2)

Nx = 400
dx = 20e-9
dt = dx / (2 * c0)

n = np.ones(Nx)
n[220:] = 3.0
eps_r = n**2

Ez = np.zeros(Nx)
Hy = np.zeros(Nx - 1)

source_pos = 80
t0 = 40
spread = 12

Nt = 600

plt.ion()
fig, ax = plt.subplots()

for t in range(Nt):

    Hy += (dt / (mu0 * dx)) * (Ez[1:] - Ez[:-1])
    Ez[1:-1] += (dt / (eps0 * eps_r[1:-1] * dx)) * (Hy[1:] - Hy[:-1])

    Ez[source_pos] += np.exp(-((t - t0) / spread) ** 2)

    Ez[0] = Ez[1]
    Ez[-1] = Ez[-2]

    if t % 5 == 0:
        ax.clear()
        ax.plot(Ez)
        ax.axvline(220, linestyle='--')
        ax.set_xlim(0, Nx)
        ax.set_ylim(-1.5, 1.5)
        plt.pause(0.01)

plt.ioff()
plt.show()

6. Durchführung der Simulation

Das Verhalten der Simulation ist in ??–?? dargestellt. Zunächst breitet sich der Gauß-Impuls ungehindert durch den Bereich \(n=1\) aus. Wenn er den Bereich mit dem veränderten Brechungsindex erreicht, teilt sich das Feld in zwei Komponenten auf: ein reflektierter Impuls, der sich zurück in den Bereich mit dem niedrigeren Brechungsindex bewegt, und ein transmittierter Impuls, der in den Bereich mit dem höheren Brechungsindex eintritt. Da der Brechungsindex auf \(n=3\) ansteigt, verlangsamt sich die transmittierte Welle und ihre Wellenlänge wird kürzer, während die reflektierte Komponente ihre ursprünglichen Ausbreitungseigenschaften beibehält. Dies entspricht genau dem erwarteten physikalischen Verhalten und zeigt, dass selbst diese minimale Implementierung die wesentlichen physikalischen Aspekte der Wellenausbreitung, Reflexion und Transmission erfasst.

11. Was ist das nächste Ziel?

Die oben genannten Beispiele für MATLAB und Python zeigen, wie einfach das FDTD-Verfahren im Prinzip ist, aber es handelt sich dabei nur um vereinfachte Implementierungen. In der Praxis erfordern realistische optische Simulationen große Berechnungsraster, lange Rechenzeiten und eine sorgfältige Behandlung von Randbedingungen und Quellen. Um das FDTD-Verfahren für nützliche, praxisrelevante Probleme einzusetzen, benötigen Sie einen Solver, der hochoptimiert, parallelisiert und sowohl auf CPUs als auch auf GPUs ausführbar ist.

Genau das bietet OghmaNano: einen umfassenden 3D-FDTD-Solver, der in C geschrieben und mit OpenCL beschleunigt ist. Dadurch kann dasselbe Modell nahtlos auf Mehrkern-CPUs oder GPUs ausgeführt werden. Anstatt einen eigenen Solver von Grund auf neu zu entwickeln, können Sie die hier gezeigten Ideen übernehmen und sie sofort auf realistische photonische Strukturen, Geräte und Probleme der Wellenoptik anwenden.

Um loszulegen, können Sie die Seite FDTD in OghmaNano aufrufen, um einen Überblick über die Funktionen zu erhalten, sich die Informationen zu was FDTD ist oder die vollständige mathematische Ableitung noch einmal ansehen oder direkt zu den praktischen Einstellungen übergehen, indem Sie Lichtquellen, Randbedingungen und Detektoren verwenden.

Es stehen auch zahlreiche sofort einsatzbereite Beispiele zur Verfügung, darunter Ausbreitung in freiem Raum, Doppelspaltexperiment, Photonische Kristallwellenleiter, Bragg-Gitter, Fabry-Pérot-Resonatoren, Geneigte Wellenleiterflächen, Verjüngte Strukturen, Mach-Zehnder-Bauelemente, Ringresonatoren und Antireflexbeschichtungen.

Kurz gesagt: Die grundlegende Methode ist einfach, wie oben gezeigt – aber OghmaNano ist so konzipiert, dass es diese Einfachheit nutzt und sie zu einem praktischen, leistungsstarken Werkzeug für realistische optische Simulationen ausbaut.