Lo mejor de dos mundos: DGTD

FDTD es intuitivo y broadband, pero su rejilla cartesiana escalonea las curvas. FEM tiene mallas elegantes que siguen la geometría, pero trabaja en dominio frecuencial — una simulación por frecuencia. ¿Y si pudieras tener ambas cosas? Mallas triangulares adaptativas Y evolución temporal. Eso es DGTD.

El problema del escalonado

En el artículo 04 vimos FDTD corriendo en tiempo real. Funciona, pero tiene un problema fundamental: la rejilla de Yee es cartesiana. Cuando el objeto tiene curvas — un cilindro, una esfera, un resonador de anillo — la rejilla lo aproxima con escalones. Eso introduce errores de fase que limitan la precisión a segundo orden, por mucho que refines la rejilla.

Para un resonador óptico de alto $Q$ , esos errores de fase son fatales. La frecuencia de resonancia converge lentamente, y el factor de calidad $Q$ puede estar mal por órdenes de magnitud incluso con rejillas finas (Niegemann et al., J. Opt. A, 2009).

Compara las dos aproximaciones del mismo cilindro:

Explorar

Resolución 10

Mismo cilindro, misma resolución. FDTD escalonea las curvas (rojo). DGTD las sigue fielmente (verde).

A la izquierda, FDTD: la línea roja punteada es el contorno real. Los bloques azules son las celdas que la rejilla clasifica como «dentro del cilindro». Las escaleras son inevitables. A la derecha, DGTD: los triángulos se adaptan al contorno y lo siguen sin escalones. Aumenta la resolución y verás cómo el FDTD sigue escalonando mientras el DGTD refina suavemente.

¿Cómo funciona?

DGTD toma prestado de tres tradiciones:

De FEM: mallas no estructuradas (triángulos en 2D, tetraedros en 3D) con funciones base polinómicas de orden $p$ en cada elemento.
De FDTD: evolución temporal explícita (típicamente Runge-Kutta de orden 4). No hay sistemas lineales grandes que resolver en cada paso de tiempo.
De volúmenes finitos: el acoplamiento entre elementos vecinos se maneja mediante flujos numéricos en las interfaces. Como las funciones base son locales a cada elemento (de ahí «discontinuas»), los campos pueden ser diferentes a cada lado de una interfaz — la información se intercambia solo a través del flujo.

Matemáticamente, en cada elemento $\Delta$ se proyectan las ecuaciones de Maxwell sobre funciones test locales $L_i$ :

\int_\Delta \left(\mathscr{Q} \cdot \partial_t \mathbf{q}_N + \nabla \cdot \mathbf{F}(\mathbf{q}_N)\right) \cdot L_i \, d^3r = 0

donde $\mathbf{q} = (E, H)^T$ es el vector de estado y $\mathbf{F}$ es el flujo de Maxwell. La clave: esta ecuación es completamente local — solo involucra campos dentro del elemento. El acoplamiento con los vecinos entra como un término de flujo en la frontera del elemento, que reemplaza $\nabla \cdot \mathbf{F}$ por un flujo numérico (típicamente upwind) evaluado en las caras compartidas.

¿Por qué «discontinuo»?

En FEM convencional, las funciones base son continuas entre elementos — los campos se comparten en los nodos. Esto crea un sistema de ecuaciones global (matriz dispersa grande) que hay que resolver.

En DG, las funciones base son locales a cada elemento. Los campos pueden saltar en la frontera entre elementos — de ahí «discontinuo». Esto hace que la matriz de masas sea bloque-diagonal (un bloque por elemento), trivialmente invertible. El resultado: stepping temporal explícito sin resolver sistemas lineales grandes.

El precio: necesitas un mecanismo de acoplamiento (el flujo numérico) para que la información pase de un elemento a otro. Pero ese flujo es local — solo involucra los dos elementos que comparten la cara.

Convergencia exponencial

La ventaja más espectacular de DGTD sobre FDTD es la convergencia. FDTD converge como $O(h^2)$ — segundo orden. Duplicar la resolución (y multiplicar el coste por 16 en 3D) solo gana un factor 4 en precisión.

DGTD con polinomios de orden $p$ converge como $O(h^{p+1})$ — y para problemas suaves, la convergencia es exponencial en $p$ . Esto significa que subir el orden polinómico de 3 a 4 puede mejorar la precisión más que duplicar la malla.

El paper de Niegemann et al. (2009) lo demuestra con un resonador cilíndrico: FDTD necesita ~128 celdas/radio para bajar el error en $Q$ al 1%. DGTD con orden $p = 3$ lo consigue con ~8 elementos/radio — usando una fracción de la memoria.

CFL en mallas no estructuradas

DGTD, como FDTD, es un método explícito — y tiene su propia condición CFL. Pero en mallas no estructuradas, la condición es más restrictiva:

\Delta t \leq C \, \frac{h_{\min}}{c \, p^2}

donde $h_{\min}$ es el tamaño del elemento más pequeño, $p$ el orden polinómico, y $C$ una constante que depende de la geometría del elemento. El factor $p^2$ es el precio de la convergencia de alto orden: a mayor $p$ , paso temporal más pequeño.

El problema: si la malla tiene unos pocos elementos muy pequeños (cerca de una punta o una rendija), esos elementos dictan el paso temporal de toda la simulación. La solución es el local time stepping (LTS): cada elemento avanza con su propio $\Delta t$ , proporcional a su tamaño. Los elementos grandes dan pasos más largos; los pequeños, más cortos. Esto puede acelerar la simulación 10–100× en problemas con mallas muy heterogéneas.

¿Cuándo usar DGTD?

Resonadores de alto Q: microdiscos, anillos, cavidades fotónicas. FDTD no puede capturar Q > 10⁴ de forma fiable; DGTD sí.
Geometrías curvas: nanopartículas esféricas o elipsoidales, guías de onda curvas, lentes. DGTD las discretiza sin escaleras.
Problemas multiscala: la malla puede ser gruesa en el vacío y ultra-fina cerca de puntas, rendijas o interfaces. La adaptatividad h-p (refinar malla Y subir orden) es natural.
Broadband: como FDTD, un pulso corto te da todo el espectro. Pero con mejor precisión.

El precio: DGTD es más complejo de implementar que FDTD (generación de malla, funciones base de alto orden, flujos numéricos). Pero las implementaciones open-source (la librería de Hesthaven & Warburton, Nodal DG; el código del grupo de Busch, BEAST) y los solvers comerciales lo han democratizado.

DGTD en la práctica

El grupo de Kurt Busch (Humboldt-Universität zu Berlin) ha sido pionero en aplicar DGTD a nanofotónica. Sus trabajos demuestran ventajas claras en:

Resonadores de cristal fotónico: cavidades con $Q > 10^6$ , donde FDTD no puede converger (Busch et al., Laser & Photonics Reviews, 2011).
Plasmónica: nanotriángulos, bow-ties, nanoantenas con puntas afiladas. La malla refina en la punta sin penalizar el paso temporal global (con LTS).
Metamateriales: celdas unitarias con geometrías complejas (split-ring resonators, fish-nets) donde las escalas van de 10 nm a varios μm.

El paisaje completo

Con DGTD cerramos un panorama que cubre todas las necesidades de la simulación en fotónica:

TMM: 1D, instantáneo — capas planas.
Mie/Rayleigh: analítico — esferas y partículas pequeñas.
FDTD: tiempo, rejilla cartesiana — rápido, intuitivo.
DDA: dipolos en espacio abierto — nanopartículas arbitrarias.
FEM: frecuencia, mallas adaptativas — cavidades y guías.
BEM: superficie, espacio abierto — plasmónica, campos cercanos.
DGTD: tiempo + mallas adaptativas — resonadores de alto Q, multiscala.

No hay un método universal. Pero con estos siete (+1), tienes las herramientas para atacar prácticamente cualquier problema de nanofotónica. El arte — como siempre — está en elegir el correcto.

Ejercicios

Ejercicio 1

Un resonador de microdisco tiene un factor de calidad $Q = 10^5$ . Eso significa que la luz da $Q/2\pi \approx 16\,000$ vueltas dentro del disco antes de disiparse. Si el disco tiene radio $R = 5$ μm y opera a $\lambda = 1.55$ μm, ¿cuántos periodos ópticos dura la simulación temporal? ¿Y cuántos pasos de tiempo si $\Delta t \approx T/20$ (donde $T$ es el periodo óptico)?

Solución

Periodo óptico: $T = \lambda/c = 1.55 \times 10^{-6} / 3 \times 10^8 \approx 5.2$ fs. Tiempo de vida del fotón en el resonador: $\tau = Q \cdot T / 2\pi = 10^5 \times 5.2 / 6.28 \approx 83\,000$ fs = 83 ps.

Periodos ópticos: $\tau/T \approx 16\,000$ . Pasos de tiempo: $16\,000 \times 20 = 320\,000$ . Con FDTD (rejilla fina para reducir staircasing) cada paso es caro. Con DGTD (malla gruesa + alto orden), cada paso cuesta menos y la convergencia es mejor. Pero 320.000 pasos son 320.000 pasos — por eso los resonadores de alto Q son un desafío computacional incluso con DGTD.

Ejercicio 2

Compara los costes de convergencia. FDTD tiene error $O(h^2)$ ; DGTD con orden $p$ tiene error $O(h^{p+1})$ . Si quieres reducir el error por un factor 100 (dos órdenes de magnitud):

Con FDTD, ¿por cuánto necesitas reducir $h$ ? ¿Cuánto sube el coste en 3D?
Con DGTD a $p = 3$ , ¿por cuánto necesitas reducir $h$ ? ¿Cuánto sube el coste?

Solución

FDTD: error $\propto h^2$ . Reducir 100×: $h \to h/10$ . Coste 3D: $(10)^3 \times 10 = 10^4\times$ más ( $N^3$ celdas × $N$ pasos por CFL).

DGTD p=3: error $\propto h^4$ . Reducir 100×: $h \to h/\sqrt[4]{100} \approx h/3.16$ . Coste 3D: $(3.16)^3 \times 3.16 \approx 100\times$ más.

10.000× vs 100× — dos órdenes de magnitud de diferencia. Esto explica por qué DGTD es indispensable para problemas que requieren alta precisión: resonadores de alto Q, propagación larga, o geometrías con features sub-λ.

Resumen en frío · Óptica numérica

Todo lo que el curso deja utilizable para dimensionar una simulación: qué ecuación resuelve cada método, qué criterio fija la discretización y qué coste se paga por ella. Si para decidir cuántas celdas por λ hacen falta o cuánta memoria va a comer el problema hay que volver al texto a buscar un umbral, esta tabla ha fallado.

Qué	Fórmula o valor	Dónde
Parámetro de tamaño	x = 2πa/λ = ka. Rayleigh si x ≪ 1, resonancias de Mie en x ~ 1. Mie (1908) resuelve la esfera homogénea en armónicos esféricos, exacta a cualquier x	art. 01
Primera resonancia de Mie	x ≈ 4 con vidrio (n = 1.50) y x ≈ 1.5 con silicio (n = 3.50): resuena cuando dentro de la partícula cabe una λ/n	art. 01, ej. 1
Paradoja de extinción	Q_ext → 2 cuando x ≫ 1, no 1: la esfera intercepta el doble de luz que su sección geométrica	art. 01
Fresnel y Snell	r_s = (n₁cos θ₁ − n₂cos θ₂)/(n₁cos θ₁ + n₂cos θ₂) · r_p = (n₂cos θ₁ − n₁cos θ₂)/(n₂cos θ₁ + n₁cos θ₂), con n₁ sen θ₁ = n₂ sen θ₂ y R = \|r\|²	art. 02
Incidencia normal y Brewster	R = ((n₁ − n₂)/(n₁ + n₂))²: aire/vidrio (1.00/1.52) da 4.3% por cara, y las 15 superficies de un objetivo dejan pasar 0.957¹⁵ ≈ 52%. tan θ_B = n₂/n₁ anula r_p	art. 02
Antirreflejante de una capa	Espesor d = λ₀/4n_f e índice ideal n_f = √(n₀·n_s); aire/vidrio pide 1.23 y el MgF₂ da 1.38, así que el mínimo se queda en 1.3% en vez de cero	art. 02, ej. 1
Matriz característica de una capa	M_j = [[cos δ_j, −(i/n_j)·sen δ_j], [−i·n_j·sen δ_j, cos δ_j]] con δ_j = 2π n_j d_j/λ; N capas → M = M₁M₂···M_N, coste lineal: 100 capas en microsegundos	art. 02
Reflectancia del apilamiento	r = (n₀M₀₀ + n₀n_sM₀₁ − M₁₀ − n_sM₁₁)/(n₀M₀₀ + n₀n_sM₀₁ + M₁₀ + n_sM₁₁), y R = \|r\|². La misma fórmula para 2 capas que para 200	art. 02
Espejo de Bragg	Pares λ/4 de alto/bajo índice; con TiO₂/SiO₂ (2.30/1.46), 3 pares pasan del 95%, ~5 del 99% y 10 del 99.9%	art. 02, ej. 2
Fabry-Perot y finesse	Dos Bragg con espaciador λ/2: ventana estrecha dentro de la banda prohibida, F = π√R/(1 − R) con R la de cada espejo	art. 02, ej. 3
Rayleigh / cuasiestática	Q_sca = (8/3)·x⁴·\|(m² − 1)/(m² + 2)\|², vía α = 4πa³(ε − ε_m)/(ε + 2ε_m). A tamaño fijo Q_sca ∝ λ⁻⁴: el azul se dispersa ~10 veces más que el rojo	art. 03
Hasta dónde vale Rayleigh	x < 0.3–0.5. Una partícula de a = 40 nm a λ = 500 nm da x ≈ 0.50, justo en el borde: ~5% de desvío con n = 1.50 y mucho más con n = 3.50	art. 03, ej. 1
LSPR y biosensado	Resuena donde Re(ε) = −2ε_m. Oro: ~520 nm en aire, ~560 nm en agua (n = 1.33) y ~10 nm más de 1.33 a 1.40 → sensibilidad local ≈ 10/0.07 = 140 nm/RIU. La plata cae en el UV/azul y es más estrecha	art. 03, ej. 2
Maxwell-Garnett (1904)	ε_eff = ε_h·(ε_i + 2ε_h + 2f(ε_i − ε_h))/(ε_i + 2ε_h − f(ε_i − ε_h)); vale para inclusiones diluidas, f < 0.3	art. 03
Bruggeman (1935)	f·(ε_i − ε_eff)/(ε_i + 2ε_eff) + (1 − f)·(ε_h − ε_eff)/(ε_h + 2ε_eff) = 0, simétrico. TiO₂ en polímero (2.30/1.50) a f = 0.25: n_eff 1.70 frente a 1.67 de Maxwell-Garnett	art. 03, ej. 3
Rejilla de Yee	E en los nodos y H en los centros de celda, actualizados en leapfrog a medio paso: Maxwell a segundo orden sin interpolar (Yee, 1966)	art. 04
CFL y resolución de rejilla	Estabilidad: c·Δt ≤ Δx/√d, o sea S = cΔt/Δx ≤ 1 (1D), ≤ 0.707 (2D), ≤ 0.577 (3D). Precisión: 10–20 celdas por λ para error de fase < 1%, y 30 o más en propagaciones largas	art. 04
PML	Coordenada estirada al plano complejo, x → x + iσ(x)/ω (Bérenger, 1994): 8–16 celdas con conductividad creciente dan reflexiones menores que 10⁻⁶	art. 04
Coste temporal de FDTD en 3D	O(N⁴) = N³ celdas × N pasos. Con Δx = 2 nm la CFL 3D da Δt ≤ 3.85 as, y simular 50 fs son ~13 000 pasos	art. 04, ej. 1
Memoria de FDTD en 3D	6 campos × 8 bytes por celda: (5 μm)³ con Δx = 5 nm son 10⁹ celdas = 48 GB, fuera de un portátil de 16–32 GB	art. 04, ej. 3
DDA: el sistema	(α⁻¹ − G)·p = E_inc, matriz densa de 3N × 3N, con α_j = (3V/4π)·(ε_j − 1)/(ε_j + 2) (Purcell y Pennypacker, 1973)	art. 05
Corrección LDR	α_LDR⁻¹ = α_CM⁻¹ − (2/3)i k³ − (k²/d)·b₁ (Draine y Goodman, 1993); sin el término radiativo, DDA viola el teorema óptico	art. 05
Cuántos dipolos hacen falta	\|m\|·k·d < 1 con k = 2π/λ: vidrio (\|m\| = 1.5) ~10 dipolos por λ, silicio (\|m\| ≈ 4) ~25, metales (\|m\| ~ 10) ~60. Oro a 600 nm (\|m\| = 3.3): d < 29 nm	art. 05, ej. 1
Precisión de DDA frente a Mie	Con la corrección LDR y \|m\|kd ≤ 0.5, el error en Q_ext es típicamente < 1%	art. 05
Memoria de DDA	9N² entradas complejas de 16 bytes: N = 10⁴ → 14.4 GB; N = 10⁵ → 1.44 TB. Por eso es obligatorio FFT (O(N log N)) con solucionador iterativo	art. 05, ej. 2
FEM: el sistema	∇×(μ_r⁻¹ ∇×E) − k₀²·ε_r·E = 0 pasa, vía formulación débil, a (K − k₀²M)·e = f, con K y M dispersas	art. 06
Elementos de arista	Funciones base en las aristas y no en los nodos (Nédélec, 1980): imponen solo la continuidad tangencial y eliminan las soluciones espurias	art. 06
Coste de FEM	Matrices dispersas: los solucionadores directos (MUMPS, PARDISO) escalan O(N^1.5) en 3D. Convergencia O(h^(p+1)), exponencial con refinamiento hp	art. 06
Elementos y DOF de una malla	~6–10 elementos por λ con p = 1. Una malla 2D de N_T triángulos tiene ≈ 1.5·N_T aristas, una por DOF: 10 000 triángulos → 15 000 DOF, 30 000 con p = 2	art. 06, ej. 1 y 2
BEM: solo la frontera	Incógnitas J_s = n̂ × H y M_s = −n̂ × E sobre la superficie, con G(r, r′) = e^(ik\|r − r′\|)/(4π\|r − r′\|); formulación PMCHWT	art. 07
Cuánto ahorra discretizar solo el borde	Cubo con L/h = 100: (L/h)³ = 10⁶ incógnitas de volumen frente a 6(L/h)² = 6×10⁴ de superficie, 17× menos	art. 07, ej. 2
El precio de la matriz densa	Almacenar O(N²), resolver O(N³): esas 6×10⁴ incógnitas son 3.6×10⁹ entradas contra 2×10⁷ del FEM disperso. El FMM lo baja a O(N log N)	art. 07, ej. 2
DGTD: de dónde viene cada pieza	Malla no estructurada de FEM + avance temporal explícito (Runge-Kutta 4) de FDTD + flujos numéricos upwind de volúmenes finitos. La matriz de masas queda bloque-diagonal: no hay sistema global que resolver	art. 08
Convergencia de DGTD	O(h^(p+1)), exponencial en p si el problema es suave. Resonador cilíndrico: FDTD pide ~128 celdas por radio para bajar del 1% de error en Q; DGTD con p = 3, ~8 elementos por radio	art. 08
CFL en malla no estructurada	Δt ≤ C·h_min/(c·p²): el elemento más pequeño dicta el paso de toda la simulación. El local time stepping da a cada elemento su Δt y acelera 10–100× en mallas heterogéneas	art. 08
Lo que cuesta un factor 100 de precisión	En 3D: FDTD, con error O(h²), exige h/10 y multiplica el coste por 10⁴; DGTD con p = 3, con O(h⁴), exige h/3.16 y lo multiplica por ~100	art. 08, ej. 2
Microdisco de alto Q	Q = 10⁵ a λ = 1.55 μm: T = λ/c ≈ 5.2 fs, τ = Q·T/2π ≈ 83 ps y ~16 000 periodos ópticos; con Δt = T/20, 320 000 pasos	art. 08, ej. 1
El paisaje completo	TMM capas 1D · Mie/Rayleigh esferas y partículas pequeñas · FDTD tiempo en rejilla cartesiana · DDA nanopartículas arbitrarias en espacio abierto · FEM frecuencia con malla adaptativa · BEM superficie y campo cercano · DGTD tiempo + malla, alto Q y multiescala	art. 08