Numerical modeling of solar wind-magnetosphere interactions: MHD and beyond-MHD models

Yuxi Chen

University of Michigan

Outline

• Physics processes control SW-M interactions

​

• Global numerical models for SW-M coupling
• Fluid models (MHD, Extended MHD, high order moment methods)
• Kinetic models (Fully kinetic, Hybrid)
• Scaling kinetic scales

​

• Capabilities of global kinetic models: using MHD-AEPIC as examples
• Magnetopause reconnection
• Shock

2

Physical processes control SW-M coupling

3

Dorelli+2007

Hasegawa+2004

• Southward IMF
• Reconnection at low latitude magnetopause
• Northward IMF
• Reconnection at high latitude magnetopause
• Kelvin–Helmholtz instability (KHI) at flank magnetopause
• During the non-linear phase of KHI, solar wind plasma enters magnetosphere through either reconnection or diffusion

SW-M coupling: the role of shock

4

Balogh+2013

• Turbulent flow, waves, and jets originating from the shock may propagate downstream, reach magnetopause and impact SW-M coupling, for example, trigger reconnection

Tools for investigating SW-M interactions

• Satellite data
• In-situ observations: Cluster, MMS…
• Remote sensing: IMAGE, SMILE…
• Ground-based observations
• SuperDARN radar
• Ground-based imager
• Magnetometer stations
• Numerical simulations
• Local simulations with simplified geometries and boundary conditions
• Global simulations

5

SW-M interactions are governed by Vlasov-Maxwell equations

6

High order moment fluid equations

MHD equations

Moment methods

Direct solvers

Assumptions

Represent velocity space with spectrum: spectral methods

Derive ideal MHD from Vlasov equation

7

• Fluid equations for each species can be obtained from Vlasov equation by multiplying moments and integrating over the velocity space
• Equation of the n-th moment depends on the the (n+1)-th. Equation chain should be truncated somewhere.
• The equation chain above is valid for collision plasma of any distributions. Not limited to Maxwellian
• Obtain ideal MHD with some assumptions:
• quasi charge neutrality (eliminate electron continuity equation)
• ignoring electron inertial (obtain Ohm’s law from electron momentum equation)
• isotropic pressure (pressure tensor -> scalar)
• truncate the moment equations by assuming the heat flux Q = 0 (requires symmetries for velocity space distributions, but do NOT require to be Maxwellian)

Simplification with assumptions

Physical capabilities and limitations of ideal MHD models

• Physical capabilities for SW-M coupling simulations
• Global structures of magnetosphere. Location and shape of bow shock and magnetopause
• Generation and propagation of MHD waves (Alfven wave, slow/fast waves) and instabilities (KH instability)
• Reasonably good at simulating the global response of magnetopause reconnection (although the local reconnection mechanism is unphysical)
• Physical limitations for SW-M coupling simulations
• Local reconnection physics is not physical
• Kinetic shock structures are missing: reflection and acceleration of particles, generation of kinetic waves.

8

Numerical accuracy of a model

• Equations are usually solved on discretized mesh

​

​

​

• Discretization introduces discretization error, which is the difference between the original equation and the equation solved numerically.
• How to reduce the error?
• Remove the low order discretization errors --> high order scheme
• Reduce cell size Δx. Number of cells ~ (Δx)³
• How to improve simulation efficiency?
• Adaptive Mesh Refinement (AMR): reduce the total number of cells
• …

​

9

Reducing numerical errors: High-order schemes

10

• High-order scheme reduces numerical diffusion/dissipation, enables development of KHI with coarse grid

11

Merkin+2013

Guo+2010

Reducing numerical errors: High-order schemes

LFM/GAMERA 

CAS group

Reducing numerical errors: Adaptive Mesh Refinement (AMR)

• Refine the regions of interest and/or the regions with large spatial gradients (boundary layers)
• Dynamic refinement
• For SW-M coupling, usually refined magnetopause

12

Gombosi+2003

Chen+2017

Scales in numerical models

• Scales in an ideal MHD model
• Ideal MHD equations do not any intrinsic length scale
• The location and shape of the boundary layers, i.e., bow shock and magnetopause are determined by the strength of the dipole field (inner boundary condition) and solar wind parameters (outer boundary condition)
• In ideal MHD, the thickness of the boundary layers are proportional to the simulation grid size
• Hall-MHD
• Introduce ion inertial length. Should be resolved to incorporate Hall physics
• Hall-MHD is much more computationally expensive than ideal MHD
• Models with kinetic ions
• Ion gyro-radius ~ ion inertial length
• Models with kinetic ions and electrons
• Electron gyro-radius ~ electron inertial length
• Debye length

​

​

13

Derive ideal MHD from Vlasov equation

14

• Obtain ideal MHD with some assumptions:
• quasi charge neutrality (eliminate electron continuity equation)
• ignoring electron inertial (obtain Ohm’s law from electron momentum equation)
• isotropic pressure (pressure tensor -> scalar)
• truncate the moment equations by assuming the heat flux Q = 0 (requires symmetries for velocity space distributions, but do NOT require to be Maxwellian)
• …

Simplification with assumptions

Derive ideal MHD from Vlasov equation

15

• Obtain ideal MHD with some assumptions:
• quasi charge neutrality (eliminate electron continuity equation)
• ignoring electron inertial (obtain Ohm’s law from electron momentum equation)
• isotropic pressure (pressure tensor -> scalar)
• truncate the moment equations by assuming the heat flux Q = 0 (requires symmetries for velocity space distributions, but do NOT require to be Maxwellian)
• …

Simplification with assumptions

Extended MHD: anisotropic pressure

• Solve for both parallel and perpendicular pressure
• Global models: BATS-R-US

​

16

Meng+2012

Derive ideal MHD from Vlasov equation

17

• Obtain ideal MHD with some assumptions:
• quasi charge neutrality (eliminate electron continuity equation)
• ignoring electron inertial (obtain Ohm’s law from electron momentum equation)
• isotropic pressure (pressure tensor -> scalar)
• truncate the moment equations by assuming the heat flux Q = 0 (requires symmetries for velocity space distributions, but do NOT require to be Maxwellian)
• …

Simplification with assumptions

Moment methods

18

• Equations
• Solve continuity, momentum and pressure equations for both electrons and ions
• Obtain E from Ampere’s law instead of Ohm’s law
• Solve for isotropic pressure: 5-moment method
• Solve for P_par and P_perp: 6-moment method
• Solver for full pressure tensor: 10-moment method
• Light wave and Langmuir oscillations are included
• Global models:
• GkeyII (5- and 10-moment)
• BATS-R-US (5- and 6-moment)
• GkeyII has been applied to simulated Ganymede’s and Mercury’s magnetosphere
• Have not been applied to model Earth’s magnetosphere yet.

Wang+2015

SW-M interactions are governed by Vlasov-Maxwell equations

19

High order moment fluid equations

MHD equations

Moment methods

Direct solvers

Assumptions

Represent velocity space with spectrum: spectral methods

Resolve velocity space with different methods

20

Palmroth+2018

Vencels, 2016

• Solve for the change of the phase space density in each cell
• Noise free
• Slower than PIC

• Update the velocity and location of each macro-particle
• Statistical noise
• Faster than Vlasov solvers

• Represent the velocity space distribution with the sum of a serial of Hermite basis functions
• Solve for the change of the basis function coefficients

# of velocity space cells per spatial cell >> # of macro-particles per cell

Memory requirement for simulating Earth’s magnetosphere with PIC

• Size of Earth’s magnetosphere: 100 R_E x 100 R_E x 100 R_E
• Electron gyro-radius in the magnetosheath is ~6km ~1/1000 R_E
• Cell size ~ electron scale ~ 1/1000 R_E
• Cell number in each direction ~ 100,000
• Cell number = 10⁵x10⁵x10⁵ =10¹⁵
• 200 macro-particles per cell
• Store 6 real numbers (velocity, location) for each particle
• A double precision real number is represented with 8 bytes

​

​

Memory requirement: 10¹⁵ x 200 x 6 x 8 ~ 10¹⁹ bytes ~ 10¹⁰ GB

A personal laptop has 16/32 GB memory

10¹⁰ GB ~ 1 billion x personal laptop

​

21

SW-M interactions are governed by Vlasov-Maxwell equations

22

High order moment fluid equations

MHD equations

Moment methods

Direct solvers

Assumptions

Represent velocity space with spectrum: spectral methods

Hybrid methods: kinetic ions and fluid electrons

23

ANGIE3D (Auburn University): PIC hybrid

Palmroth+2018

Vlasiator: Vlasov hybrid

• Electron is treated as a massless fluid
• Electron kinetic processes are eliminated and no need to resolve the electron scales
• Electric field is obtained from Ohm’s law instead of Ampere’s law

Embed local PIC code into global fluid model: MHD-EPIC

24

• MHD with embedded PIC model (MHD-EPIC): combine the efficiency of the global fluid code with the physics capabilities of the local PIC code
• Two-way coupled
• MHD provides the initial state and boundary conditions for PIC
• PIC overwrites the overlapped MHD cells

From MHD-EPIC to MHD-AEPIC

25

Chen+2023

MHD with embedded particle-in-cell (MHD-EPIC)

• PIC domains are static boxes

MHD with adaptively embedded particle-in-cell (MHD-AEPIC)

• Arbitrarily shaped PIC region
• Allows dynamic adaptation

SW-M interactions are governed by Vlasov-Maxwell equations

26

High order moment fluid equations

MHD equations

Moment methods

Direct solvers

Assumptions

Represent velocity space with spectrum: spectral methods

Reduce separation between electron scale and ion scale

• The real ratio between the ion scale and the electron scale is sqrt(1836) ~ 40. Difficult to resolve electron scale in 3D since the simulation domain is usually 1 or 2 orders larger than the ion scale
• Artificially increase electron scale by reducing the ion-electron mass ratio m_i/m_e
• Reconnection process weakly depends on the mass ratio
• A reduced mass ratio is widely used in almost all the PIC simulations

27

Shay+1998

28

Reduce separation between global scale and ion scale

Toth+2017

Outline

• Physics processes control SW-M interactions

​

• Global numerical models for SW-M coupling
• Fluid models (MHD, Extended MHD, high order moment methods)
• Kinetic models (Fully kinetic, Hybrid)
• Scaling kinetic scales

​

• Capabilities of global kinetic models: using MHD-AEPIC as examples
• Magnetopause reconnection
• Shock

29

MHD-EPIC simulation of magnetopause reconnection

30

• Show the kinetic electron feature and the evolution of FTEs in one simulation

Chen+2017

Reconnection onset and X-line expansion

31

• Beautiful results, but are they physical?

Compare simulations with SuperDARN data�(Wei Zhang, Toshi Nishimura, Yuxi Chen)

• MHD-AEPIC vs. observations: Similar flow structures and expansion speeds

​

32

Poleward flow in ionosphere

Kinetic shock: turbulent magnetosheath flow

33

• Quasi-perpendicular shock
• Extreme case with high-Mach number (M_A ~ 50)
• Turbulent magnetosheath

Kinetic shock: Reflection of both ions and electrons

34

Summary

• Physics required for descripting SW-M interactions
• Ideal MHD can capture global structures
• Beyond-MHD effects are required for modeling shock and reconnection
• Improve numerical accuracy
• High-order accurate methods
• Adaptive mesh refinement (AMR)
• Beyond ideal MHD models
• Extended MHD (Hall term, anisotropic pressure…)
• High moments methods (5-, 6- or 10-moment)
• Kinetic models.
• Resolve velocity space with macro-particle (PIC), mesh (Vlasov solver) or spectrum (spectral methods)
• Hybrid methods
• Traditional hybrid models: kinetic ions and fluid electrons
• Embed fully kinetic code into a global fluid model (MHD-AEPIC)

35

• Speed up simulations by artificially reducing the separation between scales
• Reduce the separation between electron scale and ion scale by reducing ion-electron mass ratio (m_i/m_e)
• Reduce the separation between global scale and ion scale by artificially increasing ion inertial length
• MHD-AEPIC simulations
• Evolution of FTEs
• Reconnection onset and X-line expansion. Compare with SuperDARN observations (ongoing)
• Include shock in the PIC region. Investigate how the kinetic shock impact SW-M coupling (ongoing)

​