Ion-Neutral Fluids

Introduction

The ion–neutral module evolves a partially ionized plasma as two coupled fluids: a magnetized ion fluid (solved with the MHD solver) and a neutral fluid (solved with the hydro solver). The two fluids exchange momentum through a collisional drag term and, optionally, mass through ionization and recombination. Because the drag can be stiff, the coupling is integrated implicitly using an ImEx (implicit–explicit) Runge–Kutta scheme, while the fluxes of each fluid are still handled explicitly.

The physical model for the momentum and mass exchange is

\[\mathbf{F} = -\gamma \rho_i \rho_n (\mathbf{u}_i - \mathbf{u}_n) + \xi \rho_n \mathbf{u}_n - \alpha \rho_i^2 \mathbf{u}_i, \qquad G = \xi \rho_n - \alpha \rho_i^2,\]

where subscripts \(i\) and \(n\) denote the ion and neutral fluids, \(\gamma\) is the drag coefficient, \(\xi\) the ionization rate, and \(\alpha\) the recombination rate.

Inputs

Ion–neutral runs require three blocks in the input file:

  • <hydro> — configures the neutral fluid (see Hydrodynamics);

  • <mhd> — configures the ion fluid (see Magnetohydrodynamics);

  • <ion-neutral> — configures the coupling (note the hyphen in the block name).

The <ion-neutral> block has three parameters:

  • drag_coeff (mandatory): the ion–neutral drag coupling coefficient \(\gamma\).

  • ionization_coeff (default 0.0): the ionization rate \(\xi\).

  • recombination_coeff (default 0.0): the recombination rate \(\alpha\).

Note

  • Whenever both <hydro> and <mhd> are present, the <ion-neutral> block is required (and vice versa) — a single run cannot otherwise contain both fluids.

  • The time integrator must be an ImEx scheme: set <time>/integrator to imex2, imex2+, or imex3. Explicit integrators (rk1–rk4) will abort.

Outputs

There are no ion-neutral-specific output variables. Instead, write the two fluids using their respective standard variables (see Hydrodynamics, Magnetohydrodynamics, and Outputs):

  • the neutral fluid via hydro_w (or hydro_u, individual components, etc.);

  • the ion fluid via mhd_w_bcc (primitives plus the cell-centered magnetic field).

Example: a perpendicular C-shock

A C-shock (continuous shock) is a steady, magnetically moderated shock structure that arises in weakly ionized media when the neutrals stream through the magnetized ions. AthenaK provides a 1D test at inputs/ion-neutral/perp_cshock.athinput using the built-in cshock problem generator.

1. Build (the cshock generator is built in, so the default build suffices):

cmake -B build
cmake --build build -j

2. Note the setup. The input file couples an isothermal ion (MHD) and neutral (hydro) fluid with an ImEx integrator and a unit drag coefficient:

<time>
integrator = imex3      # ImEx integrator is required

<ion-neutral>
drag_coeff = 1.0

<hydro>
eos             = isothermal
iso_sound_speed = 1.0

<mhd>
eos             = isothermal
iso_sound_speed = 1.0

3. Slightly modify it. The shipped file writes .tab line-outs; to make the drag coupling more visible you can raise the drag coefficient, which compresses the shock front:

<ion-neutral>
drag_coeff = 2.0

4. Run and visualize. For 1D runs the .tab text output is convenient. Each line-out file (named <basename>.<variable>.<number>.tab) has columns gid i x1v dens velx vely velz, so the neutral density profile can be plotted directly:

./build/src/athena -i inputs/ion-neutral/perp_cshock.athinput -d run

import numpy as np
import matplotlib.pyplot as plt

d = np.loadtxt("run/tab/cshock.hydro_w.00040.tab")   # columns: gid i x1v dens velx ...
x, dens, velx = d[:, 2], d[:, 3], d[:, 4]

fig, ax = plt.subplots()
ax.plot(x, dens, label="neutral density")
ax.set_xlabel("x"); ax.set_ylabel("density")
fig.savefig("cshock_dens.png", dpi=150)

The neutral density rises smoothly across the C-shock rather than through a sharp discontinuity (shown here together with the neutral velocity):

Neutral density profile across a perpendicular C-shock

Overlaying the ion density (from the mhd_w_bcc output) shows how the two fluids decouple in the shock front, with the separation set by drag_coeff.