Geometric numerical discretization of electromagnetic quasineutral models

Planteamiento del Problema
La simulación numérica de plasmas magnetizados es esencial para la fusión de energía, el clima espacial y la astrofísica, pero enfrenta desafíos debido a la vasta separación de escalas de longitud y tiempo. Aunque los modelos cinéticos basados en el sistema Vlasov-Maxwell son necesarios para capturar distribuciones no maxwellianas e interacciones onda-partícula, resolver la longitud de Debye y las ondas electromagnéticas de alta frecuencia (como las ondas de Langmuir y de luz) es computacionalmente prohibitivo para muchas aplicaciones. Los modelos cuasin neutrales, que asumen la neutralidad de carga ($\rho=0$) y filtran estas escalas pequeñas, ofrecen una alternativa más eficiente. Sin embargo, los enfoques numéricos existentes para sistemas cuasin neutrales a menudo dependen de esquemas de preservación asintótica que requieren amortiguamiento numérico para suprimir los modos de alta frecuencia a medida que la longitud de Debye tiende a cero, o carecen de la preservación de la estructura geométrica (como la conservación de la energía y las restricciones de divergencia exacta) que se encuentra en los métodos modernos de Particle-in-Cell (PIC) que preservan la estructura.

Metodología
Los autores extienden el marco de PIC electromagnético geométrico, GEMPICX, para resolver las ecuaciones de Vlasov-Maxwell cuasin neutrales y totalmente cinéticas. La metodología se construye sobre los siguientes pilares:

• Marco Geométrico: El esquema utiliza diferencias finitas miméticas en mallas duales (mallas primal y dual). Las variables de campo se definen en espacios discretos (nodos, aristas, caras, volúmenes) que forman un complejo de de Rham discreto, asegurando la preservación exacta de las propiedades topológicas como la condición de divergencia nula del campo magnético.
• Principio Variacional: La discretización se deriva de un principio de acción discreto basado en el Lagrangiano por Tronci & Camporeale (2015). Esto asegura que las ecuaciones de movimiento semi-discretas para partículas y campos sean consistentes con la estructura Hamiltoniana subyacente.
• Formulación Cuasin Neutral: A diferencia de los enfoques de preservación asintótica que resuelven el sistema completo de Vlasov-Maxwell y toman el límite $\epsilon_0 \to 0$, este método se deriva directamente en el límite cuasin neutral. En consecuencia, el término de corriente de desplazamiento está ausente y el campo eléctrico no puede evolucionarse mediante una derivada temporal. En su lugar, el campo eléctrico se obtiene implícitamente en cada paso de tiempo resolviendo una ecuación de curl-curl discretizada (Ecuación 2.25), la cual es no singular y definida positiva.
• Corrección de la Divergencia de Corriente: Para satisfacer la condición de cuasin neutralidad $\nabla \cdot \mathbf{J} = 0$ con precisión de máquina (lo cual no se satisface automáticamente por la deposición de partículas debido al ruido estadístico), se introduce un multiplicador de Lagrange. Esto se manifiesta como un campo eléctrico "corrector de divergencia" ($\mathbf{E}^\star = -\nabla \psi^\star$) resuelto mediante una ecuación tipo Poisson (Ecuación 6.8) para ajustar las velocidades de las partículas antes del paso de tiempo principal.
• Paso de Tiempo: El algoritmo emplea un esquema dividido de segundo orden de precisión. Implica corregir la divergencia de la corriente, resolver para el potencial vectorial $\mathbf{A}$ (usando la calada de Coulomb $\nabla \cdot \mathbf{A} = 0$), resolver para el campo eléctrico $\mathbf{E}$, y finalmente actualizar las posiciones y velocidades de las partículas utilizando un empujador de Boris (Boris pusher).

Contribuciones Clave

1. Extensión de GEMPICX: El artículo presenta la primera extensión del marco GEMPIC de diferencias finitas miméticas al límite cuasin neutral para electrones e iones totalmente cinéticos.
2. Derivación Directa Cuasin Neutral: El esquema se deriva directamente en el límite cuasin neutral en lugar de depender de la consistencia asintótica de un resolvedor completo. Esto asegura la ausencia de ondas electromagnéticas y de Langmuir espurias en la solución numérica sin depender del amortiguamiento numérico.
3. Preservación de la Estructura: La formulación semi-discreta demuestra matemáticamente la conservación de la suma de las energías magnética y cinética. También mantiene el requisito de divergencia cero para la densidad de corriente con precisión de máquina mediante el uso de un multiplicador de Lagrange.
4. Resolvedor de Campo Eléctrico Implícito: Los autores derivan un sistema lineal robusto (Ecuación 5.12) para obtener el campo eléctrico implícitamente, evitando la necesidad de potenciales eléctricos en la actualización del campo y evitando los problemas de singularidad asociados con el operador curl-curl estándar en ausencia de corriente de desplazamiento.
5. Escalabilidad: La implementación se integra en el framework AMReX con paralelización MPI, demostrando capacidades de escalado fuerte y débil en recursos de computación de alto rendimiento (HPC).

Resultados
El marco fue validado mediante una serie de pruebas numéricas con condiciones de contorno periódicas:

• Relaciones de Dispersión: Las simulaciones de una sola especie (solo electrones) y de doble especie (electrón-ión) reprodujeron con éxito las relaciones de dispersión teóricas para ondas que se propagan paralela y perpendicularmente al campo magnético de fondo.
  • Propagación Paralela: El esquema capturó correctamente la onda de ciclotrón de electrones (ECW) y la onda de ciclotrón de iones (ICW), así como las estructuras cónicas de los modos de Alfvén-ciclotrón de orden superior amortiguados. Como era de esperar, las ondas de Langmuir y las ondas electromagnéticas L/R estuvieron ausentes.
  • Propagación Perpendicular: Las simulaciones reprodujeron las ondas de Bernstein y las resonancias de armónicos de ciclotrón. En los casos de doble especie, se observaron las ondas de Bernstein de iones y la onda híbrida. Los modos O y X estuvieron correctamente ausentes debido al supuesto de cuasin neutralidad.
• Amortiguamiento de Ondas: Una prueba que involucró una perturbación sinusoidal simple demostró la capacidad del código para capturar la tasa de amortiguamiento predicha teóricamente de una onda de ciclotrón, validando la precisión del esquema en el manejo de procesos cinéticos de orden superior.
• Eficiencia Computacional: Una comparación entre la resolución del sistema curl-curl completo para el campo eléctrico (Ecuación 5.12) y una aproximación diagonal más económica (Ecuación 5.13) mostró que la aproximación redujo el costo computacional en aproximadamente un 25% con una pérdida de precisión insignificante para los casos probados.
• Escalabilidad: Los estudios de escalado fuerte y débil en una arquitectura de memoria distribuida mostraron una aceleración casi ideal hasta 32 tareas MPI y mantuvieron una alta eficiencia hasta 128 tareas, lo que indica que el código es adecuado para simulaciones a gran escala.

Significancia y Reivindicaciones
Los autores afirman que este trabajo logra cerrar la brecha entre los métodos de preservación de estructura geométrica y el modelado de plasma cuasin neutral. Al derivar el esquema directamente en el límite cuasin neutral, el método filtra inherentemente los modos de alta frecuencia, ofreciendo una alternativa más eficiente y físicamente consistente a los esquemas de preservación asintótica que dependen del amortiguamiento. La preservación de la estructura de de Rham discreta y la conservación de la energía (al nivel semi-discreto) aseguran la estabilidad a largo plazo y la fidelidad física.

El documento señala modestamente las limitaciones y áreas para trabajos futuros:

• Condiciones de Contorno: La implementación actual depende de condiciones de contorno periódicas para imponer fácilmente la calada de Coulomb e invertir el operador de Laplace. Extender esto a contornos no periódicos (por ejemplo, para dispositivos de fusión) requiere análisis de compatibilidad y restricciones adicionales.
• Orden de Precisión: El trabajo utiliza operadores de Hodge de segundo orden. Aunque es posible utilizar operadores de orden superior dentro del marco mimético, su beneficio específico para simulaciones cuasin neutrales (que carecen de modos de alta frecuencia) sigue siendo una cuestión abierta.
• Manejo de Choques: El método asume campos suaves. Aunque los métodos PIC generalmente manejan los choques a través de las distribuciones de partículas, el comportamiento de este esquema de discretización geométrica específico en regímenes con gradientes fuertes o choques requiere mayor investigación.
• Discretización Temporal: La conservación de la energía es exacta solo al nivel semi-discreto; tras la discretización temporal, la conservación es precisa solo hasta el orden del integrador de tiempo. Se sugiere como extensión futura un esquema de discretización temporal totalmente variacional para lograr una conservación de energía discreta exacta.