FDTD en OghmaNano

1. Introducción

Finite-Difference Time-Domain (FDTD) es un método numérico de onda completa para resolver las ecuaciones de Maxwell directamente en el dominio temporal sobre una rejilla espacial discreta. A diferencia del trazado de rayos, los solucionadores de modos o los enfoques de matriz de transferencia, FDTD no realiza suposiciones geométricas ni modales. En su lugar, los campos eléctricos y magnéticos se avanzan explícitamente en el tiempo mediante la discretización de las ecuaciones de rotacional —la ley de Ampère y la ley de Faraday— de modo que las ondas electromagnéticas se propaguen, dispersen, interfieran y difracten de forma natural a través del dominio de simulación.

\[ \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}. \]

En el método de diferencias finitas en el dominio temporal (FDTD), las ecuaciones de Maxwell se resuelven directamente en el dominio temporal mediante el reemplazo de derivadas espaciales por diferencias finitas sobre una rejilla escalonada (normalmente una red de Yee) y avanzando los campos en el tiempo. Debido a que este enfoque prácticamente no realiza suposiciones sobre la dirección de propagación, la estructura modal o el comportamiento del campo, una única simulación FDTD puede capturar respuesta de banda ancha, difracción e interferencia, efectos de campo cercano, geometría sub-longitud de onda y fenómenos electromagnéticos transitorios dentro de un mismo marco. Esta ausencia de suposiciones simplificadoras es lo que hace que FDTD sea tan general, pero también lo que lo hace computacionalmente costoso: la rejilla espacial debe resolver la escala de longitud física más pequeña, el paso temporal está fuertemente restringido por la estabilidad numérica, y las simulaciones tridimensionales se vuelven rápidamente intensivas en memoria y tiempo incluso para geometrías que parecen modestas en tamaño.

2. Cuándo usar FDTD y cuándo no usar FDTD

En la mayoría de los problemas ópticos prácticos, FDTD no debería ser el primer método al que recurrir. Para sistemas ópticos de gran escala, largas distancias de propagación o situaciones dominadas por óptica geométrica, FDTD es una elección extremadamente ineficiente. Sistemas de imagen, lentes macroscópicas, aperturas y ensamblajes ópticos con dimensiones del orden de milímetros a centímetros se tratan mucho mejor utilizando trazado de rayos, mientras que estructuras multicapa planas como células solares de película delgada, pilas OLED y recubrimientos se modelan con mayor precisión y eficiencia mediante métodos de matriz de transferencia. En estos regímenes, FDTD actúa como un martillo computacional: aumenta drásticamente el tiempo de ejecución y el uso de memoria mientras aporta poca o ninguna información física adicional.

En la práctica, FDTD debe considerarse como un método de último recurso, reservado para casos en los que modelos más simples realmente dejan de ser válidos. Su fortaleza reside en situaciones donde el campo electromagnético completo debe resolverse explícitamente, como la difracción de redes sub-longitud de onda, cristales fotónicos, metasuperficies, interfaces fuertemente desordenadas o rugosas y efectos de acoplamiento de campo cercano que no pueden capturarse mediante óptica de rayos o modelos de medios en capas. Incluso en estos casos, FDTD debe aplicarse de forma limitada y deliberada, con una comprensión clara de su coste numérico, restricciones de estabilidad y requisitos de memoria. Las secciones siguientes derivan por tanto las ecuaciones discretas de actualización utilizadas en OghmaNano directamente a partir de las leyes de Ampère y Faraday, no para fomentar su uso indiscriminado, sino para hacer explícito tanto lo que FDTD puede hacer como por qué debe utilizarse con moderación.

3. Fundamentos teóricos

Esta sección del manual tiene como objetivo describir el código FDTD en su totalidad con derivaciones detalladas para ayudar a la comprensión/detección de errores.

La ley de Ampère se expresa como \[\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}\]

que puede expandirse como

\[\sigma E_{x} + \epsilon \frac{\partial E_{x}}{\partial t} = \frac{\partial H_{z}}{\partial y}-\frac{\partial H_{y}}{\partial z}\]

\[\sigma E_{y} + \epsilon \frac{\partial E_{y}}{\partial t} = -\frac{\partial H_{z}}{\partial x}+\frac{\partial H_{x}}{\partial z}\]

\[\sigma E_{z} + \epsilon \frac{\partial E_{z}}{\partial t} = \frac{\partial H_{y}}{\partial x}-\frac{\partial H_{x}}{\partial y}\]

Para el caso \(\frac{\partial}{\partial y}=0\)

\[\begin{split} &\sigma E_{x} + \epsilon \frac{\partial E_{x}}{\partial t} =-\frac{\partial H_{y}}{\partial z}\\ &\sigma E_{y} + \epsilon \frac{\partial E_{y}}{\partial t} = -\frac{\partial H_{z}}{\partial x}+\frac{\partial H_{x}}{\partial z}\\ &\sigma E_{z} + \epsilon \frac{\partial E_{z}}{\partial t} = \frac{\partial H_{y}}{\partial x} \end{split}\]

para \(E_{x}\) \[\begin{split} &\sigma E_{x} + \epsilon \frac{\partial E_{x}}{\partial t} =-\frac{\partial H_{y}}{\partial z}\\ &\sigma \frac{E_{x}^{t+1}[]+E_{x}^{t}[]}{2} + \epsilon \frac{E_{x}^{t+1}[]-E_{x}^{t}[]}{\Delta t} = -\frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}\\ &\sigma \frac{E_{x}^{t+1}[]}{2} + \epsilon \frac{E_{x}^{t+1}[]}{\Delta t} = -\frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{x}^{t}[]}{2}+\epsilon \frac{E_{x}^{t}[]}{\Delta t}\\ &\sigma \frac{E_{x}^{t+1}[]}{2} + \epsilon \frac{E_{x}^{t+1}[]}{\Delta t} = -\frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{x}^{t}[]}{2}+\epsilon \frac{E_{x}^{t}[]}{\Delta t}\\ & \frac{\sigma \Delta t + 2 \epsilon }{ 2 \Delta t}E_{x}^{t+1}[] = -\frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{x}^{t}[]}{2}+\epsilon \frac{E_{x}^{t}[]}{\Delta t}\\ & E_{x}^{t+1}[] = \left ( -\frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{x}^{t}[]}{2}+\epsilon \frac{E_{x}^{t}[]}{\Delta t} \right ) \frac{2 \Delta t}{\sigma \Delta t + 2 \epsilon} \end{split}\]

para \(E_{y}\) \[\begin{split} &\sigma E_{y} + \epsilon \frac{\partial E_{y}}{\partial t} = -\frac{\partial H_{z}}{\partial x}+\frac{\partial H_{x}}{\partial z}\\ &\sigma \frac{E_{y}^{t+1}[]+E_{y}^{t}[]}{2} + \epsilon \frac{E_{y}^{t+1}[]-E_{y}^{t}[]}{\Delta t} = -\frac{H_{z}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{z}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}+\frac{H_{x}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{x}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}\\ &\sigma \frac{E_{y}^{t+1}[]}{2} + \epsilon \frac{E_{y}^{t+1}[]}{\Delta t} = -\frac{H_{z}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{z}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}+\frac{H_{x}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{x}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{y}^{t}[]}{2} + \epsilon \frac{E_{y}^{t}[]}{\Delta t}\\ &E_{y}^{t+1}[] = \left ( -\frac{H_{z}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{z}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}+\frac{H_{x}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{x}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{y}^{t}[]}{2} + \epsilon \frac{E_{y}^{t}[]}{\Delta t} \right ) \frac{2 \Delta t}{\sigma \Delta t + 2 \epsilon} \end{split}\]

para \(E_{z}\) \[\begin{split} &\sigma E_{z} + \epsilon \frac{\partial E_{z}}{\partial t} = \frac{\partial H_{y}}{\partial x}\\ &\sigma \frac{E_{z}^{t+1}[]+E_{z}^{t}[]}{2} + \epsilon \frac{E_{z}^{t+1}[]-E_{z}^{t}[]}{\Delta t} = \frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}\\ &\sigma \frac{E_{z}^{t+1}[]}{2} + \epsilon \frac{E_{z}^{t+1}[]}{\Delta t} = \frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}-\sigma \frac{E_{z}^{t}[]}{2} + \epsilon \frac{E_{z}^{t}[]}{\Delta t}\\ &E_{z}^{t+1}[]= \left ( \frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}-\sigma \frac{E_{z}^{t}[]}{2} + \epsilon \frac{E_{z}^{t}[]}{\Delta t} \right ) \frac{2 \Delta t}{\sigma \Delta t + 2 \epsilon}\\ \end{split}\]

La ley de Faraday se expresa como \[-\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}\]

que puede expandirse para dar:

\[-\sigma_{m} H_{x} - \mu \frac{\partial H_{x}}{\partial t} = \frac{\partial E_{z}}{\partial y}-\frac{\partial E_{y}}{\partial z}\]

\[-\sigma_{m} H_{y} - \mu \frac{\partial H_{y}}{\partial t} = -\frac{\partial E_{z}}{\partial x}+\frac{\partial E_{x}}{\partial z}\]

\[-\sigma_{m} H_{z} - \mu \frac{\partial H_{z}}{\partial t} = \frac{\partial E_{y}}{\partial x}-\frac{\partial E_{x}}{\partial y}\]

Con \(\sigma_m=0\) y \(\frac{\partial}{\partial y}=0\)

\[\begin{split} &\frac{\partial H_{x}}{\partial t} = \frac{1}{\mu} \left ( \frac{\partial E_{y}}{\partial z} \right )\\ &\frac{\partial H_{y}}{\partial t} = \frac{1}{\mu} \left ( \frac{\partial E_{z}}{\partial x}-\frac{\partial E_{x}}{\partial z} \right )\\ &\frac{\partial H_{z}}{\partial t} = - \frac{1}{\mu} \left ( \frac{\partial E_{y}}{\partial x} \right ) \end{split}\]

que al discretizarse da

\[\begin{split} & H_{x}^{t+1} = \frac{1}{\mu} \left ( \frac{E_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-E_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z} \right ) \Delta t + H_{x}^{t}[]\\ & H_{y}^{t+1} = \frac{1}{\mu} \left ( \frac{E_{z}^{t+\frac{1}{2}}[\frac{1}{2}]-E_{z}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}-\frac{E_{x}^{t+\frac{1}{2}}[\frac{1}{2}]-E_{x}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z} \right ) \Delta t+ H_{y}^{t}[]\\ & H_{z}^{t+1} = \frac{1}{\mu} \left ( - \frac{E_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-E_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x} \right ) \Delta t + H_{x}^{z}[] \end{split}\]

4. Solución numérica en un ordenador

Las ecuaciones derivadas anteriormente describen la evolución en tiempo continuo y espacio continuo de los campos eléctricos y magnéticos. Para resolverlas en un ordenador, tanto el espacio como el tiempo se discretizan. El dominio de simulación se divide en una rejilla regular de celdas, y los componentes del campo se almacenan en posiciones espaciales discretas y pasos temporales discretos. En la formulación FDTD estándar, los campos eléctricos y magnéticos están escalonados tanto en espacio como en tiempo, de modo que cada componente del campo se evalúa donde es naturalmente requerido por las ecuaciones de rotacional.

En la práctica, esto significa que los componentes del campo eléctrico se almacenan en pasos temporales enteros (\(t, t+\Delta t, t+2\Delta t\)), mientras que los componentes del campo magnético se almacenan en pasos temporales semi-enteros (\(t+\tfrac{1}{2}\Delta t, t+\tfrac{3}{2}\Delta t\)). Derivadas espaciales como \(\partial / \partial x\) y \(\partial / \partial z\) se reemplazan por diferencias finitas centradas entre puntos vecinos de la rejilla. Este escalonamiento conduce a un esquema de actualización simple y explícito en el que los campos se adelantan entre sí en el tiempo.

En cada paso temporal, el campo magnético se actualiza primero utilizando la forma discretizada de la ley de Faraday, con el rotacional del campo eléctrico evaluado a partir de los valores actuales del campo eléctrico. Una vez que el campo magnético ha avanzado medio paso temporal, el campo eléctrico se actualiza utilizando la forma discretizada de la ley de Ampère, incluyendo parámetros del material como permitividad, permeabilidad y conductividad. No se requiere resolver una matriz global: cada componente del campo en cada celda de la rejilla se actualiza utilizando únicamente sus vecinos locales del paso temporal anterior.

Este procedimiento explícito de avance temporal se repite hasta que se alcanza el tiempo de simulación deseado. Debido a que el método está en el dominio temporal, una única simulación contiene de forma natural información sobre un amplio rango de frecuencias; cantidades en el dominio de la frecuencia como espectros o distribuciones de campo en estado estacionario se obtienen normalmente registrando los valores del campo en función del tiempo y procesándolos posteriormente mediante transformadas de Fourier. Fuentes, fronteras y capas absorbentes (como las capas perfectamente adaptadas) se incorporan directamente en las ecuaciones de actualización en las posiciones apropiadas de la rejilla.

La simplicidad de las ecuaciones de actualización hace que FDTD sea conceptualmente sencillo y altamente paralelizable, pero también impone estrictas restricciones numéricas. El paso temporal debe satisfacer una condición de estabilidad que depende del espaciamiento de la rejilla espacial y de las propiedades del material, y el uso de memoria escala directamente con el número de celdas de la rejilla y los componentes de campo almacenados. Estas restricciones limitan en última instancia el tamaño y la duración de las simulaciones que pueden realizarse en la práctica y motivan una elección cuidadosa de la resolución de la rejilla, el tamaño del dominio y la complejidad del modelo.