Maxwell en una rejilla: FDTD

FDTD — Finite-Difference Time-Domain — es el método más directo de todos: toma las ecuaciones de Maxwell, sustituye las derivadas por diferencias finitas, y avanza paso a paso en el tiempo. No hay matrices que invertir, no hay funciones base que elegir. Solo la rejilla, las ecuaciones, y el tiempo. FDTD fue propuesto por Kane Yee en 1966 y sigue siendo, sesenta años después, el método más usado en electromagnetismo computacional.

La rejilla de Yee

La clave de FDTD es la rejilla de Yee (1966). En lugar de poner todos los campos en los mismos puntos, Yee intercala los campos eléctricos y magnéticos: $E$ vive en los nodos de la rejilla, $H$ en los centros de las celdas (y viceversa según la componente). Esta disposición garantiza que las ecuaciones de Maxwell se satisfacen automáticamente a segundo orden.

Para el modo TE en 2D (campo $E_z$ perpendicular al plano), las ecuaciones de actualización son:

H_x^{n+\frac{1}{2}} = H_x^{n-\frac{1}{2}} - \frac{\Delta t}{\mu_0} \frac{\partial E_z}{\partial y}

H_y^{n+\frac{1}{2}} = H_y^{n-\frac{1}{2}} + \frac{\Delta t}{\mu_0} \frac{\partial E_z}{\partial x}

E_z^{n+1} = E_z^{n} + \frac{\Delta t}{\varepsilon} \left(\frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y}\right)

Derivación: de Maxwell a las ecuaciones de actualización

Partimos de las ecuaciones de Maxwell rotacionales:

\nabla \times \mathbf{E} = -\mu_0 \frac{\partial \mathbf{H}}{\partial t} \qquad \nabla \times \mathbf{H} = \varepsilon \frac{\partial \mathbf{E}}{\partial t}

Para el modo TE en 2D ( $E_z, H_x, H_y$ , campos independientes de $z$ ), las componentes del rotacional dan:

\frac{\partial E_z}{\partial y} = -\mu_0 \frac{\partial H_x}{\partial t} \qquad -\frac{\partial E_z}{\partial x} = -\mu_0 \frac{\partial H_y}{\partial t} \qquad \frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y} = \varepsilon \frac{\partial E_z}{\partial t}

Ahora discretizamos. Cada derivada temporal se reemplaza por una diferencia centrada entre pasos de tiempo: $\partial H_x / \partial t \approx (H_x^{n+\frac{1}{2}} - H_x^{n-\frac{1}{2}}) / \Delta t$ . Cada derivada espacial, por la diferencia entre celdas vecinas: $\partial E_z / \partial y \approx (E_z^{i,j+1} - E_z^{i,j}) / \Delta y$ . Despejando el campo al paso siguiente, se obtienen directamente las ecuaciones de actualización. La rejilla de Yee garantiza que todos los campos están evaluados en los puntos correctos — sin interpolar.

Cada derivada parcial se aproxima con la diferencia entre celdas vecinas. El algoritmo es un leapfrog: actualiza $H$ a medio paso, luego $E$ a paso completo. Repite.

La condición CFL: estabilidad

FDTD es un método explícito — no resuelve sistemas lineales, solo aplica las ecuaciones de actualización. Pero eso tiene un precio: la estabilidad no está garantizada. Si el paso temporal es demasiado grande, los errores crecen exponencialmente y la simulación explota.

La condición de estabilidad es la condición CFL (Courant-Friedrichs-Lewy):

c \, \Delta t \leq \frac{\Delta x}{\sqrt{d}}

donde $d$ es la dimensión (1, 2 o 3) y $c$ la velocidad de la luz en el medio. La interpretación es geométrica: la onda no puede viajar más de una celda por paso de tiempo. Si lo hace, la información «escapa» de la rejilla numérica y la simulación diverge.

¿De dónde sale la condición CFL?

Sustituimos una onda plana $e^{i(kx - \omega t)}$ en las ecuaciones discretizadas. Después de algo de álgebra, la relación de dispersión numérica del esquema FDTD en 1D es:

\sin^2\!\left(\frac{\omega \Delta t}{2}\right) = \left(\frac{c \Delta t}{\Delta x}\right)^2 \sin^2\!\left(\frac{k \Delta x}{2}\right)

Para que $\omega$ sea real (no hay crecimiento exponencial), el lado derecho no puede exceder 1. El máximo de $\sin^2(k\Delta x/2)$ es 1 (para $k = \pi/\Delta x$ , la mayor frecuencia que la rejilla puede representar). Por tanto:

\left(\frac{c \Delta t}{\Delta x}\right)^2 \leq 1 \quad \Rightarrow \quad c \Delta t \leq \Delta x

En $d$ dimensiones, las contribuciones de cada eje se suman, dando $c \Delta t \leq \Delta x / \sqrt{d}$ . El número $S = c\Delta t / \Delta x$ se llama número de Courant: debe ser ≤ 1 en 1D, ≤ 0.707 en 2D, ≤ 0.577 en 3D.

Dispersión numérica

Incluso cuando la simulación es estable, la rejilla introduce un artefacto: la velocidad de la onda en la rejilla no es exactamente $c$ . Depende de la dirección de propagación y de la frecuencia — es la dispersión numérica.

La regla práctica: necesitas al menos 10–20 celdas por longitud de onda ( $\Delta x \leq \lambda/10$ a $\lambda/20$ ) para mantener el error de fase por debajo del 1%. Para propagaciones largas (muchos $\lambda$ de distancia), los errores se acumulan y puedes necesitar 30 o más celdas por $\lambda$ .

La simulación en vivo

Aquí tienes una simulación FDTD real corriendo en tu navegador. Una onda plana entra por la izquierda y encuentra un obstáculo. Observa la difracción, la reflexión, y cómo la onda se ralentiza dentro del dieléctrico:

Explorar

Obstáculo

Permitividad 4.0

Campo E_z (tiempo real)

La onda se dispersa al encontrar el obstáculo. Mayor permitividad = más reflexión y menor velocidad dentro.

Prueba la doble rendija: verás el patrón de difracción formarse en tiempo real — exactamente como en el artículo 02 del módulo anterior, pero ahora calculado directamente desde Maxwell, sin aproximaciones.

PML: simulando el infinito

La rejilla FDTD es finita, pero el espacio no. Si la onda llega al borde, se refleja — un artefacto que arruina la simulación. La solución estándar son las PML (Perfectly Matched Layers), inventadas por Bérenger en 1994.

La idea es rodear el dominio con una capa de material artificial que absorbe la onda sin reflejarla. El truco matemático: en la PML, la coordenada espacial se estira al plano complejo ( $x \to x + i\sigma(x)/\omega$ ). Esto convierte la onda propagante en una onda que decae exponencialmente, sin cambiar la impedancia del medio — sin reflexión.

En la práctica, una PML de 8–16 celdas con conductividad creciente (perfil polinómico o geométrico) da reflexiones menores que $10^{-6}$ . La simulación de arriba usa exactamente esta técnica: el borde de la rejilla tiene una rampa de absorción que «traga» la onda.

Ventajas y limitaciones

FDTD tiene ventajas potentes:

Intuitivo: «ves» la onda propagarse. El resultado es una película, no un número. No hay método más visual.
Broadband: una sola simulación con un pulso corto te da la respuesta en todo el espectro vía FFT temporal.
Geometría arbitraria: cualquier distribución de $\varepsilon(x,y)$ cabe en la rejilla. Materiales dispersivos (Drude, Lorentz) se incluyen con ecuaciones auxiliares.
Paralelizable: cada celda solo interactúa con sus vecinas inmediatas. Esto hace que FDTD escale muy bien en GPUs y clusters — es un método «embarazosamente paralelo».

Pero también tiene precio:

Rejilla cartesiana: las curvas se aproximan como «escaleras» ( staircasing). Para geometrías curvas, necesitas rejillas muy finas — o pasar a DGTD (Artículo 08).
Coste en 3D: escala como $O(N^4)$ — $N^3$ celdas × $N$ pasos de tiempo. Una simulación 3D de una nanopartícula de 100 nm con $\Delta x = 1$ nm necesita ~10⁶ celdas y puede tardar horas.
Dispersión numérica: la velocidad de fase depende de la dirección y de la resolución. Se acumula en propagaciones largas.

FDTD en la práctica

Los dos softwares FDTD más usados en fotónica son Lumerical FDTD Solutions (comercial, interfaz gráfica, estándar industrial) y Meep (MIT, open-source, scriptable en Python). Ambos incorporan PML, materiales dispersivos, fuentes broadband, y monitores de campo cercano/lejano. Una simulación típica de nanopartícula plasmónica en 3D con Meep: ~30 minutos en un portátil, ~3 minutos en GPU.

Ejercicios

Ejercicio 1

Quieres simular una nanopartícula de oro de 80 nm de diámetro iluminada a $\lambda = 520$ nm. Si usas $\Delta x = 2$ nm, ¿cuántas celdas por longitud de onda tienes? ¿Cuál es el paso temporal máximo permitido por la condición CFL en 3D? ¿Cuántos pasos necesitas para simular 50 fs?

Solución

Celdas por λ: $520/2 = 260$ — más que suficiente (el mínimo práctico es ~10–20).

CFL en 3D: $\Delta t \leq \Delta x / (c\sqrt{3}) = 2 \times 10^{-9} / (3 \times 10^8 \times 1.732) \approx 3.85 \times 10^{-18}$ s = 3.85 as (attosegundos).

Pasos para 50 fs: $50 \times 10^{-15} / 3.85 \times 10^{-18} \approx 13\,000$ pasos. Cada paso actualiza todas las celdas: si el dominio es $200 \times 200 \times 200 = 8 \times 10^6$ celdas, son ~10¹¹ operaciones en total. Factible en minutos en GPU.

Ejercicio 2

La simulación FDTD de arriba usa una rejilla de 150 × 150 en 2D. Si el dominio físico es de 3 μm × 3 μm, ¿cuál es $\Delta x$ ? ¿Cuántas celdas por longitud de onda tienes a $\lambda = 500$ nm? ¿Es suficiente? ¿Y a $\lambda = 300$ nm?

Solución

$\Delta x = 3000/150 = 20$ nm.

A 500 nm: $500/20 = 25$ celdas/λ — bien. A 300 nm: $300/20 = 15$ celdas/λ — aceptable pero al límite. Para UV profundo o materiales de alto índice (donde $\lambda_{\text{eff}} = \lambda/n$ puede ser mucho menor), habría que refinar.

Ejercicio 3

¿Cuánta memoria necesita una simulación FDTD 3D con dominio de $5 \times 5 \times 5$ μm³ y $\Delta x = 5$ nm? Cada celda almacena 6 componentes de campo (3 de E, 3 de H) en doble precisión (8 bytes cada una). ¿Cuántos GB? ¿Cabe en la RAM de un portátil?

Solución

Celdas por eje: $5000/5 = 1000$ . Total: $1000^3 = 10^9$ celdas.

Memoria: $10^9 \times 6 \times 8 = 48 \times 10^9$ bytes = 48 GB. No cabe en un portátil típico (16–32 GB). Habría que reducir el dominio, aumentar $\Delta x$ , o usar un cluster. Esta es la principal limitación de FDTD en 3D: el coste de memoria escala como $N^3$ .