Data-driven multi-species heat flux closures for two-stream-unstable plasmas with nonlinear sparse regression

TL;DR

This work extends data-driven heat-flux closures to multi-species collisionless plasmas by generalizing a six-term closure and predicting the three key coefficients $A_1$, $A_4$, and $A_5$ from box-averaged fluid quantities. It combines sparse regression (SINDy) with a nonlinear rational form and Bernstein polynomial basis to produce interpretable, pole-free closures, and benchmarks them against neural networks using OSIRIS PIC data. The results show that the closures capture 80–90% of heat-flux variation and 85–95% of the pressure-time change across beam, core, and combined electron populations, offering a practical path toward efficient, physics-informed fluid models that incorporate kinetic effects. The study also analyzes the theoretical consistency with multi-species linear theory and discusses avenues for improvement, including local averaging and extensions to relativistic regimes.

Abstract

The dual aims of accuracy and computational efficiency in computational plasma physics lend themselves well to the use of fluid models. The first of these goals, however, is only satisfied for such models insofar as the utilized closure can capture the neglected kinetic physics -- something which has proven challenging for multi-scale collisionless processes. In a recent article [E. R. Ingelsten et al. (2025) J. Plasma Phys. 91 E64], we used the data-driven method of sparse regression to discover a novel heat flux closure for electrostatic phenomena. Here, we generalize the six-term closure model found in that work from single- to multi-species modeling. Using data from OSIRIS particle-in-cell simulations over a range of initial conditions, we then demonstrate how the unknown coefficients in front of the three most important terms in the closure can be estimated from box-averaged fluid quantities. Both neural networks and a newly developed framework for nonlinear sparse regression are showcased. The resulting models predict the heat flux for each species with a typical accuracy of 80-90 % and regularly account for 85-95 % of the rate of change in the pressure. The models are also compared with results from multi-species linear collisionless theory.

Figures (8)

• Figure 1: Log-scale plot of training data FVU vs iteration for the neural network model estimating the combined-species $A_1$ coefficient (blue), along with a 51-iteration moving average (black). The 51-iteration average FVU attained at the end of training, $21%$, is highlighted with a red horizontal line.
• Figure 2: The time evolution of the spatial variance in (a) $n_\sigma$, (b) $V_\sigma$, (c) $p_\sigma$ and (d) $q_\sigma$, with beam quantities in dashed blue, core quantities in dotted orange and combined-species quantities in dash-dotted green, along with $\mathop{\mathrm{var}}\limits_x E$ in solid black, rescaled to fit the plotted ranges in each panel. In the depicted simulation, an initial condition of $\qty{n_b, V_\text{rel}, v_{\text{th},b}, v_{\text{th},c}} =$$\{0.17\,\bar{n},3.4e-2c, 5.0e-3c, 6.4e-3c\}$ was used. The linear growth phase, defined as the time range of linear $E$-field perturbation growth, is highlighted in red.
• Figure 3: Snapshots of (a:i) $p_\sigma$ and (a:ii) $q_\sigma$ at the end of linear growth, with beam in dashed blue, core in dotted orange, combined-species in dash-dotted green and the “extra term” in long dashed red. Additionally, we show the time evolution of the spatial variance in $p_b$ and $q_b$ (dashed blue), $p_c$ and $q_c$ (dotted orange) and the extra $p$ and $q$ terms (long dashed red), normalized to the combined-species variances $\mathop{\mathrm{var}}\limits_x p$ and $\mathop{\mathrm{var}}\limits_x q$. Pressure quantities are shown in (b:i) and heat flux quantities in (b:ii). The vertical dashed red lines in these plots mark the time of the snapshot panels. For reference we also show the time evolution of the logarithmized $E$-field variance in gray in arbitrary units, and highlight the linear growth phase in red. As before, the depicted case uses the initial condition $\qty{n_b, V_\text{rel}, v_{\text{th},b}, v_{\text{th},c}} =$$\{0.17\,\bar{n},3.4e-2c, 5.0e-3c, 6.4e-3c\}$.
• Figure 4: Time evolution of 6-term closure coefficients compared to growth rates $\gamma_E$ and $\gamma_q$ for the (a) beam species, (b) core species and (c) combined species, together with (d) the FVU over time for the discovered models of $q_b$ (dashed blue), $q_c$ (dotted orange) and $q$ (dash-dotted green). In the former three panels, $q_\text{even}$ ($q_\text{odd}$) coefficients are shown in blue (red), while $\gamma_E$ is shown in tightly dotted gray and $\gamma_q$ is shown in faded blue (orange, green) for the beam (core, combined) cases, respectively. In the latter panel, $\log \mathop{\mathrm{var}}\limits_x E$ in arbitrary units is additionally shown in gray for reference. In every subplot, the linear growth phase is highlighted in red.
• Figure 5: The two 3-term closure coefficients $A_1$ (solid blue) and $A_5$ (dash-dottd red) compared to predictions from linear theory at $\gamma = 0$ (in faded blue and red, respectively) shown for the beam (top), core (middle) and combined population (bottom) in the simulation with initial condition $\qty{n_b, V_\text{rel}, v_{\text{th},b}, v_{\text{th},c}} =$$\{0.17\,\bar{n},3.4e-2c, 5.0e-3c, 6.4e-3c\}$. Furthermore, the growth rate of the electric field perturbation is plotted in dotted gray and the linear growth phase is highlighted in red.